Dynamic analysis of a malaria reaction-diffusion model with periodic delays and vector bias

One of the most important vector-borne disease in humans is malaria, caused by Plasmodium parasite. Seasonal temperature elements have a major effect on the life development of mosquitoes and the development of parasites. In this paper, we establish and analyze a reactiondiffusion model, which includes seasonality, vector-bias, temperature-dependent extrinsic incubation period (EIP) and maturation delay in mosquitoes. In order to get the model threshold dynamics, a threshold parameter, the basic reproduction number R0 is introduced, which is the spectral radius of the next generation operator. Quantitative analysis indicates that when R0 < 1, there is a globally attractive disease-free ω-periodic solution; disease is uniformly persistent in humans and mosquitoes if R0 > 1. Numerical simulations verify the results of the theoretical analysis and discuss the effects of diffusion and seasonality. We study the relationship between the parameters in the model and R0. More importantly, how to allocate medical resources to reduce the spread of disease is explored through numerical simulations. Last but not least, we discover that when studying malaria transmission, ignoring vector-bias or assuming that the maturity period is not affected by temperature, the risk of disease transmission will be underestimate.


Introduction
Malaria is a lethal infectious disease caused by Plasmodium parasites, which infects humans through the bite of infected female geneus Anopheles mosquitoes. It is endemic in about 100 countries and an estimated 3.3 billion people are at risk of being infected in the word [1]. Therefore, it is necessary to study the spread of malaria in the population. The pioneering work was given by Ross [2] which showed that malaria is transmitted by mosquitoes and proposed the first mathematical model of this disease. Subsequently, Macdonald promoted this [3], which plays a particularly important role in understanding the disease. However, the Ross-Macdonald model is extremely simplified, ignoring many critical factors of real-world ecology and epidemiology [4]. One omission is spatial heterogeneity, for example, the distribution of urban and rural areas and the spatial movement of mosquitoes and people, which increases the spatial transmission of malaria. Hence, in the malaria epidemiological model, spatial heterogeneity must be considered [5][6][7]. In 1951, Skellam specially examined the reaction-diffusion model of species population density in a bounded habitat [8]. One advantage of reaction-diffusion models is that they can easily incorporate simple rules about edge effects. Use these features of reaction-diffusion model to study how environmental heterogeneity affects malaria [9][10][11]. However, in these studies, it seems that the following important biological aspects are rarely paid, simultaneously.
First, "vector-bias" is a typical feature of malaria transmission, which describes the different probabilities of a mosquito bitting humans. Lacroix et al. [12] proved that mosquitoes are generally more inclined to bite infected individuals. To describe this preference, Chamchod and Britton incorporated the vector-bias effect into malaria model [13]. Subsequently, a wealth of papers take into account the effect of vector-bias when studying malaria transmission [14][15][16]. Due to the different densities of mosquitoes in various regions, the infectivity of malaria is also different. The analysis of vector-bias effects is of great significance for understanding the dynamics of malaria transmission.
Secondly, in recent years, climate has become more fluctuating [17], and the rise in temperature is related to global warming. Climate warming is directly connected with the number of Anopheles, the main vector of malaria. Because temperature affects the immature development rate, immature stage survival rate, adult size, adult lifespan, blood feeding, and reproduction ability of mosquitoes [18]. It is widely accepted that both the life cycle of mosquitoes and the development of parasites are extremely sensitive to temperature. Empirical evidence indicates that the life cycle of mosquitoes experiences four different stages: eggs, larvae, pupae and adults. The first three stages are also called aquatic or immature stages [19]. In addition, a very important process in the development of parasites is the extrinsic incubation period (EIP), i.e., mosquitoes may not be able to transmit the disease to humans for a period of time after taking infected blood meal. The lifespan of female adult mosquitoes is generally 3 to 100 days, and EIP can reach 30 days [20]. These infected mosquitoes that live pass EIP will keep infectious until they die. There is considerable evidence that temperature rise will shorten the length of each immature stage [21] and EIP is extremely sensitive to seasonal changes in temperature [22,23]. In view of climate change, exploring the effect of temperature in the spread of malaria is particularly essential. Therefore, these seasonal forced EIP and mosquito life cycles need to be considered in the malaria model.
Note that mathematical models of malaria transmission rarely consider key features such as vector-bias, climate/seasonality, temperature-dependent mosquito life cycle and EIP in a spatially heterogeneous environment, simultaneously. A reaction-diffusion model with periodic delays is proposed to study the influence of the above factors on the transmission of malaria. This work is inspired by the following biological questions: (I) What is the impact of vector-bias on malaria transmission? (II) What effect does the maturity and incubation period of seasonal temperature changes have on the spread of malaria? (III) Does the movement of mosquitoes and humans have different effects for the transmission of malaria in different areas? These issues that are closely related to the spread and control of malaria need to be explored in more depth. How to add the temperature-dependent maturity period to the model is the first issue we consider, and the detailed model derivation is also the difficulty of this paper. The existence of temperature-dependent periodic delays makes it difficult to understand the asymptotic behavior of the model. We consider the stage structure, because the mosquitoes in the aquatic stage cannot fly, which makes it impossible for us to directly use the theoretical results in references [22,24]. Our theoretical results indicate the model has a threshold parameter, the basic reproduction number R 0 , which determines wether the malaria can persist or become extinct in susceptible populations.
The structure of the rest of this paper is shown below. Section 2 establishes human and mosquito populations models and malaria transmission model. This section also gives the biological significance of these parameters in models. Section 3 is devoted to well-posedness and some properties of the model. In Section 4, we first define the basic reproduction number R 0 , and then analyze the threshold dynamics of models in terms of R 0 . Section 5 gives the basic framework of the system threshold dynamics through numerical simulation, and discusses the influence of parameters on basic reproduction number, mosquitoes and humans in the process of malaria transmission. A brief discussion section concludes the paper.

Mosquito stage structure model with seasonal effect
Motivated by models in references [19,25,26], we introduce a seasonal stage-structured model of malaria to explain the cross-infection between mosquitoes and individuals. In this subsection, we first consider the mosquito population. According to the biological cycle, mosquitoes can go through several stages and have different habits at various stages. We divide them into immature aquatic stage M(t, x) and mature stage V(t, x) at time t ≥ 0 and location x in a spatial bounded domain Ω ∈ R n with smooth boundary ∂Ω. Let η(t, a, x) denote the density of the mosquito population at time t, chronological age a ≥ 0, and location x. Metz and Diekmann [27] gave through the standard argument for structured populations and spatial diffusion. Here, ∆ denotes the usual Laplacian operator.φ(·, x) represents the non-negative initial age and location distribution of mosquitoes, the term Λ v (t, x) is the recruitment rate of the immature stage, which is Hölder continuous and nonnegative function on R ×Ω, D(a) is the age-related diffusion coefficient, and d(t, a, x) indicates the death rate of mosquitoes with time t, age a and location x, which is Hölder continuous and positive on R ×Ω. Λ v (t, ·) and d(t, ·, ·) are periodic in t and have the same period ω. Let us first assume some age thresholds so that the growth of mosquitoes can be divided into discrete age stages. When modeling species dynamics with a stage structure, we divide mosquitoes into immature stage (also known as aquatic stage) and mature stage, the time delay is the mature period. The temperature T (t) is an ω-periodic function of time t. Generally speaking, the maturity age is determined by seasonal changes in temperature. We will use a time-dependent positive function τ 1 (t) to describe the duration from newborn to adult. That is, an individual at time t becomes mature only if its age exceeds τ 1 (t), in other words, the adult mosquitoes that develop at time t enter the aquatic stage at t − τ 1 (t). We assume that the maximum chronological age of aquatic mosquitoes at time t is A(t). Here and in the following, τ 1 (t) and A(t) are C 1 (a section of continuous differentiable functions) periodic functions of time t with the same periodic ω. The relationship between time-dependent threshold age and time-dependent delay can be indicated as follows: Here we introduce another quantity f to describe the maturity level of mosquitoes. Let f = f M = 0, that is, the maturity level of mosquitoes that have just entered the aquatic stage is 0. At time t, the mosquitoes reach the maximum chronological age A(t) of the aquatic stage. At this time, the maturity level is also the maximum, f V . Then one has (r 1 )dr 1 whereρ(r 1 ) is the proportion of maturity at time r 1 . Taking the derivative of the above formula with respect to t, we can get , Recall that τ 1 (t) represents the time it takes for a mosquito to go from the aquatic stage to the adult stage at time t. For the case where δ > 0 is small enough, if the mosquito is an adult mosquito at time t, then it is still an adult mosquito at time t + δ. Thus τ 1 (t + δ) < τ 1 (t) + δ, and hence, τ ′ 1 (t) = lim δ→0 τ 1 (t+δ)−τ 1 (t) δ ≤ 1. This condition also makes sense biologically. 1 − τ ′ 1 (t) > 0, is called maturation rate [26], which excludes the possibility of mosquitoes returning from the mature stage to the immature stage unless they are born. Mathematically speaking, t − τ 1 (t) is strictly increased with respect to t, that is, since the developmental rate depends only on time, it is impossible for juveniles to reach maturity before they are born.
M(t, x) and V(t, x) denote the density of mosquitoes in the aquatic and mature stages at time t and location x, respectively. The accumulative density between the two age thresholds can be recorded as the population size at each stage. In particular, for the aquatic stage M and the mature stage V (for the sake of simplicity, we use M and V instead of M(t, x) and V(t, x), respectively), we have the following mathematical expressions: Next, we put forward natural biological assumptions about birth rate, death rate and diffusion coefficient in Eq (2.1). Since immature mosquitoes cannot fly, we assume that the diffusion coefficient satisfies At each age stage, all mosquitoes undergo the same birth rate and death rate which are age-independent, then the natural death rates in the aquatic and mature stages are µ m (t, x) and µ v (t, x), respectively, i.e., where µ m (t, x) and µ v (t, x) are non-negative non-trivial and Hölder continuous on R ×Ω and ω-periodic in t. Differentiating the equations in system (2.2) with respect to time t on both sides yields In a biological sense, because no mosquito can live forever, η(t, ∞, x) = 0. η(t, 0, x) represents the inflow rate of mosquitoes in the aquatic stage at time t, then η(t, 0, x) = Λ v (t, x) is a function of t.
Mathematically, assume that the delay τ 1 (t) is continuously differentiable in [0, ∞), bounded and far away from zero and infinity. The expression of η(t, A(t), x) in the above system will be obtained by integrating along characteristics. For t ≥ A(t) and t < A(t), η(t, A(t), x) has different expressions. Without loss of generality, we study t ≥ A(t) (see, [28]), which is practicable due to the boundedness of A(t), because we are concerned with long-term behavior of population dynamics. Here, For t ≥ A(t), setting t 0 = t − A(t), h = A(t) and a 0 = 0, we have Since we are concerned about the long-term behavior of population dynamics, the model for t ≥ A(t) is as follows where ν is the outward unit normal vector on ∂Ω and ∂ ∂ν indicates the normal derivative along ν on ∂Ω, e − t t−τ 1 (t) µ m (r 2 ,x)dr 2 represents probability surviving through natural death during aquatic stage. It is easy to see that the equation V(t, x) is decoupled, so we can get the following system (2.5) for any t ≥ A(t). Alternatively, according to biological significance, M and V can be written as integral form. Note that τ 1 (t) is the maturation time of aquatic mosquitoes at time t. Thence, aquatic mosquitoes at time t contain all newly born mosquitoes at a previous time ξ 1 with ξ 1 ∈ (t − τ 1 (t), t) and survive to time t. Hence, one has Similarly, for any x ∈Ω, The first term in Eq (2.7) represents those mosquitoes that were adults at time 0 at position y and are still in V class at time t at position x by diffusion. The second term in Eq (2.7) denotes the total density of mosquitoes at time t. These mosquitoes were newly developed adults at time ξ 1 and survived until time t. These new adults developed from the aquatic stage at previous time ξ 1 − τ 1 (ξ 1 ), which grown through τ 1 (ξ 1 ) periodic of time and developed into adult stage at time ξ 1 with "maturation rate" x)V that meets the Neumann boundary condition.

The human population dynamics model
Since the incubation period of the pathogen in the human is 1-12 days, which is short compared with the human life span, we do not consider the incubation period of the pathogen in the human body. In this subsection, inspired by [13,19,20,24,29], we give a model of malaria transmission, including the mobility of mosquitoes and individuals, vector-bias, temperature-dependent maturation delay and EIP. Because mosquitoes are the vectors of malaria transmission from person to person, the modeling takes into account the dynamics of human and mosquito populations. One supposes that individuals will become susceptible again after recovery. Thus, humans include two epidemiological types: susceptible (S h ) and infected (I h ) compartments. Assume that all populations walk unbiased and random. Since climate change has almost no impact on human birth and death, according to [29], the following reaction-diffusion equation is suitable: D h > 0 is the human diffusion coefficient. Suppose the dispersal pattern is random diffusion, that is, individual walkers use a fixed step to walk randomly on the solid line [30,31]. Λ h (x) and µ h (x) are the recruitment rate and natural death rate of human at location x, respectively. In particular, we should note that the disease dose not have a significant impact on the mobility of individuals and mosquitoes. It is mathematically assumed that the diffusion coefficient D h > 0 of all individuals is the same, and the diffusion coefficient D v > 0 of all mosquitoes is recognized as the same. Assuming that no population flow crosses the boundary ∂Ω, therefore, the Neumann boundary condition is imposed: here, the definition of ν is same as Section 2.1. Following Theorems 3.1.5 and 3.1.6 in reference [32] and Section 2 in reference [29], we realize that the positive steady state N(x) ∈ C Ω , R + \{0} of systems (2.8) and (2.9) is globally attractive, where C Ω , R + is Banach space of continuous functions fromΩ to R + . According to Appendix 8.1 in reference [33], it is assumed that the total human density of at time t and location x is stable at N(x), i.e., N h (t, x) ≡ N(x), for ∀t ≥ 0 and x ∈ Ω. Then we know that S h = N(x) − I h . Here and later, for the sake of simplicity, we use S h and I h instead of S h (t, x) and I h (t, x), respectively. From a biological and epidemiological point of view, immature mosquitoes do not take part in infection, and are primarily in a waiting period, because they cannot fly and bite humans. Since adult male mosquitoes do not eat blood meal, they cannot transmit and acquire the virus, so only adult female mosquitoes are modeled. Mature female mosquitoes are divided into susceptible (S v (t, x)), exposed (E v (t, x)) and infectious (E v (t, x)), t is time and x denotes location. The exposed mosquitoes are already infected, but not infectious. From Section 2.1, we know that V = S v + E v + I v is the total density of female adult mosquitoes and satisfies reaction-diffusion system (2.5). Here and later, for the sake of simplicity, we use V, S v , E v and Vector-bias effect is a situation in which mosquitoes display the weight of preference when choosing a host. In order to add a vector-bias term to the model, parameters p and l are introduced, which are the probability of mosquitoes arriving at humans randomly, and select infectious and susceptible humans, respectively [13]. Since mosquitoes are more attracted to infected humans, we assume that the vector-bias parameter χ = p l ≥ 1. To comprehend the influences of spatial structure and climate change, we allow a space-dependent bite rate β(t, x) of female mosquitoes to be Hölder continuous non-negative non-trivial on R ×Ω with β(t, x) 0, representing the number of mosquitoes biting per unit time at location x, and β(t, ·) is a ω-period function of time t. Let B represent the mosquitoes bite incident, and A be the event of mosquitoes reaching humans infected with malaria. Then P(A) = I h N(x) , P(A c ) = 1 − I h N(x) , P(B|A) = p and P(B|A c ) = l. So the probability that an infectious human is chosen as P(A|B), then one has and which are the probability of infected individuals and susceptible individuals are bitten. Then at time t and location x, the numbers of newly infected individuals and mosquitoes per unit time are respectively, where c is the probability that an infected mosquito will spread to susceptible individuals per bite, and α is probability that an infected individual will transmission to susceptible mosquitoes per bite.
During the malaria transmission cycle, the EIP of parasite in mosquitoes is one of the most pivotal parameters. Mosquitoes can fly around during EIP. The malaria parasites undergo different growth stages before migrating to the salivary glands where can spread to humans. The degree of this development is related to environmental factors, such as temperature. EIP is exceedingly sensitive to temperature. Therefore, it is necessary to incorporate this seasonal forced incubation period when describing malaria transmission. We will use a time-dependent positive function τ 2 (t) to describe the duration from new infection to infectiousness. That is, a mosquito is infectious at time t only if its age of infection exceeds τ 2 (t), where τ 2 (t) is a C 1 periodic function in [0, ∞) with the periodic ω, and t − τ 2 (t) is strictly increasing in t.
Based on the above discussions, we will pay attention to the dynamics of the following malaria transmission reaction-diffusion model 10) for t > 0. Similar to the derivation in Section 2.1, we get Substituting E m (t, x) into system (2.10), one has 11) for t > 0, where α, p, l and c are positive constants, and ϱ(x) denotes the recovery rate of humans in position x, which is Hölder continuous and positive on Ω. Neumann boundary condition indicates that all individuals and mosquitoes still remain in this domain Ω.
Since M(t, x) and E v (t, x) of the system (2.11) are decoupled from the other equations, it is sufficient to explore the system
Before studying the global existence, go back to system (2.11) to make more information. We impose the compatibility conditions based on the biological meaning of τ 1 (t) and τ 2 (t): and Next, for any ψ ∈ D, system (2.11) has the unique solution U(t, ·) = (M(t, ·), S v (t, ·), E v (t, ·), I v (t, ·), I h (t, ·)) satisfying U 0 = ψ. Corollary 4 in reference [36] implies that S v (t, x) ≥ 0, I v (t, x) ≥ 0 and I h (t, x) ≥ 0 on maximal existence interval. According to the uniqueness of the solutions and the compatibility conditions (3.4) and (3.5), one has M(t, ·) (see(2.6)) and Hence, M(t, x) ≥ 0 and E v (t, x) ≥ 0 on the maximal interval of existence.
∈ Ω, and Q := Q ω has a strong global attractor in W h .
The proof of Lemma 3.2 is given in Appendix A.

Threshold dynamics
Here, the basic reproduction number are introduced by the theory in reference [41], and then investigate the dynamics of system (2.12).

Threshold dynamics
The current main work is to explore the extinction and persistence of malaria in a spatially heterogeneous environment. The basic reproduction number R 0 is one of the most necessary concepts in epidemiology. It is usually defined as the mean number of secondary infections during the entire period when a typical infected person is introduced into a completely susceptible population. From a biological point of view, the basic reproduction number is a characteristic of the risk of infection. More meaningfully, it describes the threshold behavior of many infectious disease models.
For every fixed t ≥ 0, letQ t be the solution map of system (4.2) on the space E 1 , i.e.,Q t (ψ) =w t (ψ), for t ≥ 0 and ∀ψ ∈ E 1 . Thus,Q :=Q ω is the Poincaré map related to system (4.2) and r(Q) is the spectral radius ofQ. The next lemma indicates that the periodic semiflowQ t is eventually strongly positive.
The proof of Lemma 4.3 is given in Appendix C. Fix an integer n 0 satisfying n 0 ω > 3τ. Obviously,Q n 0 =Q n 0 ω is strongly monotone, according to proof of Lemma 4.3. In addition, through similar arguments in reference [43, Lemma 2.6], one can get thatQ n 0 is compact. Apply the Krein-Rutman theorem to the linear operatorQ n 0 , and because r(Q n 0 ) = (r(Q)) n 0 , one hasλ = r(Q) > 0, among themλ is the simple eigenvalue ofQ and there is a eigenvector ϑ ∈ Int(E + 1 ) that is strongly positive. Based on reference [22, Lemma 3.8], one has the following lemma.
Lemma 4.4. The spectral radii of the two Poincaré mapsQ : X 2 → X 2 andQ : E 1 → E 1 are the same i.e., r(Q) = r(Q). Furthermore, R 0 − 1 has the same sign as r(Q) − 1. The fact that system has a special solution is given by the following lemma, which can be used to study long-term dynamics. The following argument is inspired by references [29,Lemma 5] and [44, Proposition 1.1].
Lemma 4.6. SetQ t to be the solution map of system (2.12) on the space E 2 , i.e.,Q t (φ) = u t (φ), for t ≥ 0 and ∀φ ∈ E 2 . Thus,Q t is an ω-periodic semiflow on E 2 , in this sense, one has (i)Q 0 = I; Furthermore,Q :=Q ω admits a strong global attractor in E 2 .
Next, we prove that the solution of system (2.12) is strictly positive.
Then ω(φ) is an internally chain transitive set for {Q n } on E 2 .
With the above claim, it is easy to obtain that H is an isolated invariant set forQ in E 2 , and W s (H) ∩ P 0 = ∅, where W s (H) is the stable set of H forQ. For every integer n and nω >τ 2 ,Q n :=Q nω is compact, thenQ is compact, andQ is asymptotically smooth on E 2 . Moreover, similar to Lemma 4.6, we get thatQ has a global attractor on E 2 . According to reference [45,Theorem 3.7],Q has a global attractor A 0 in P 0 . Based on the acyclicity theorem on uniform persistence for map [32, Theorem 1.3.1 and Remark 1.3.1], we obtain thatQ is uniformly persistent with respect to (P 0 , ∂P 0 ), in this sense, there exists an γ 1 > 0, such that lim inf n→∞ d(Q n (φ), ∂P 0 ) ≥ γ 1 , ∀φ ∈ P 0 .

Numerical
In this part, we use numerical simulations to clarify the impact of spatial heterogeneity, periodic delays, seasonality, and vector-bias on the spread of malaria, to show how to derive some understanding of epidemiology from our analysis results.

Numerical confirmation of theoretical results
Here, some simulations to support the analysis results built in Section 4 are provided. More precisely, the potential dynamics outcomes of system (2.11) are given. For convenience, we will focus our attention on the domain Ω, which is one-dimensional, denoted by [0, π] to simulate the long-time behavior of system (2.11). The basic reproduction number R 0 often (but not always) controls the extinction and persistence of the disease. Specifically, system (2.11) applies to the transmission of malaria in Maputo Province, Mozambique, a sub-Saharan African country where malaria is at risk of spreading in the country. The WHO report shows that the number of confirmed malaria cases has been increasing in recent years, and the condition deteriorated severely in 2014.
According to reference [24], suppose the total population density N(x) =N = 53 (km 2 ) −1 . The unit of time is month. Set the periodic ω = 1 year = 12 months. We set ϱ = a 1 (2.05 − cos(2x)) Month −1 , where a 1 = 0.0187 means that people living in urban areas (near the center area) can get better medical services (number of doctors and hospitals, advanced medical equipment, supply of medicines) than people in rural areas. This results in a higher recovery rate near the center. See Appendix E for the values and definitions of other parameters.
To compute R 0 , we use the numerical method proposed in reference [41, Lemma 2.5 and Remark 3.2]. Using the above parameters set, the basic reproduction number can be numerically calculate as R 0 = 1.2233 > 1, which implies that in mosquito and human populations, the disease persist. Figure 3 shows numerical plots of humans and mosquitoes compartments with initial datâ . This conforms to the conclusion of Theorem 4.9(ii). The time interval can be truncated [20,30] to prove that the existence of the positive periodic solution.
Reduce the bite rate to 0.5β(t), and increase the death rate of mosquito to 1.5µ v (t), which can be achieved by using insecticide-treated nets, then R 0 = 0.2983 < 1. Thus, immature mosquitoes and susceptible adult mosquitoes show periodic fluctuations. The exposed adult mosquitoes, infectious adult mosquitoes and infectious humans converge to 0 (Figure 4). This is coincident with the conclusion of Theorem 4.9(i).   . The long-term behavior of system (2.11) solution, when R 0 < 1.

Effects of parameters on R 0
We are interested in how system parameters affect the disease risk R 0 . In this section, we choose the parameters that are consistent with those in Section 5.1.

The impact of periodic maturity delay on disease transmission
In this subsection, we explore the impact of mature period affected by seasonal climate change on malaria transmission. The results are given in Figure 5, the green one represents the time-averaged value of maturation delay [τ 1 ] in system (2.11), where [τ 1 ] = 1 ω ω 0 τ 1 (t)dt = 0.5605 Month, and the red one indicates system (2.11) is under a time-periodic maturation delay. An intuitive result is that in these four cases, the green line is always below the red line, which fully indicates that the time-average maturity delay will underestimate the risk of malaria transmission, so the impact of climate change cannot be ignored. In the following subsections, we analyze in detail the sensitivity of R 0 to these four parameters.

The relationship between R 0 and l/p
Looking at Figure 5, we let l/p vary in [0, 1], and other parameters are same as those in Figure 3 and explore the impact of the vector-bias effect, l/p are used to represent the relative attractiveness of susceptible to infected humans. From Figure 5(a), it can be easily seen that R 0 is an increasing function of l/p. Therefore, ignoring the vector-bias effect will ultimately underestimate the risk of disease transmission. And this conclusion is valid regardless of whether the periodic or the time-average maturation delay is considered. (c) The relationship between R 0 and a 1 (d) The relationship between R 0 and and δ 1 Figure 5. The relationship between R 0 and parameters. There are two curves. The green one represents the time-average maturity delay [τ 1 ], and the red one represents τ 1 (t) affected by the periodic temperature.

The impact of interventions on R 0
In this subsection, we examine the impact of interventions on R 0 . Note that in Figure 3(e), the population density of infected individuals in urban (central) areas is less than in rural (border) areas, which prompts us to explore the relationship between human recovery rate ϱ and disease risk R 0 in a heterogeneous space. Since medical resources (advanced medical equipment, the number of hospitals and doctors, medicine supply) are abundant in urban areas than in rural areas, humans can get better medical services in urban areas with a higher recovery rate.
In this subsection, we analyze the impact of medical resources on the spread of malaria, including two aspects: (i) the impact of increasing medical resources; (ii) for fixed medical resources, the impact of different allocation measures on disease transmission. In Figure 5(b), we let a 1 in ϱ vary in [0.02, 0.5], and the other parameters are same as those in Figure 3. It proposes a strategy, i.e., in the current allocation of the medical resources, the recovery rate can be improved through efforts to develop new drugs. Figure 5(b) shows that increasing medical resources can reduce R 0 to less than 1, which means that the disease can be controlled. In this paper, we use the recovery rate π 0 ϱ(x)dx of the entire region to represent medical resources.
In the case of fixed medical resources, we explore the impact of different allocation on the spread of disease. Choose a 1 = 0.058 and a new parameter δ 1 is added to ϱ, i.e., ϱ(x) = 0.058 × (1 − (1 − δ 1 ) cos(2x)), where δ 1 ∈ [0, 1], which measures the difference in the allocation of urban and rural medical resources. Other parameters are the same as those in Figure 3. When δ 1 = 0, the gap between rural and urban medical resources is the largest, and the recovery rate in urban area is the highest (near the center of space area, i.e., x = π 2 ). As δ 1 changes from 0 to 1, medical resources are allocated to nearby rural area (near x = 0 and x = π), and the gap in medical conditions between rural and urban areas gradually narrows, and is eventually evenly distributed in space. Through calculation, it can be found that if δ 1 = 0, in order to make R 0 less than 1, a 1 ≥ 0.0637 is needed, but if δ 1 = 1, the medical resources only need 0.0576, which can make R 0 < 1. From Figure 5(c), we come to an interesting conclusion that in the case of limited medical resources, the spatial heterogeneity of the distribution of medical resources increases the risk of disease transmission. This result indicates that the distribution of medical resources dose play an important role in the design of malaria control programs.

The influence of D h on the R 0
We explore the impact of human mobility on the transmission of malaria in a heterogeneous space. Let D h vary in [0.02, 0.5] and other parameters are consistent with those in Figure 3. The result is shown in Figure 5(d), which describes the relationship between R 0 and D h . It can be directly observed from Figure 5(d) that R 0 decreases with respect to D h . When D h is very small, R 0 will drop sharply and as D h continues to increase, the rate of decline of R 0 will slow down. This seems to indicate that malaria cannot be controlled if no protective measures are taken and only by improving a high diffusion to avoid mosquito bites without taking any protective measures [49] systematically analyzes this problem.

The impact of diffusion on the final size of malaria
In Section 5.2.4, we explore the impact of diffusion on R 0 , and study the impact of diffusion on the final size of disease transmission. Set ϱ(x) = 0.04(2.05 − cos(2x)) and other parameters are consistent with those in Figure 3. When D h = 0.1 and D v = 0.0125, calculate R 0 = 0.8448, the disease is extinct. However, without considering the movement of mosquitoes and humans, that is, when the diffusion coefficients are D h = 0 and D v = 0, it is found that R 0 = 1.1635 is greater than 1, and the disease is persistent. Then, we will get the result shown in Figure 6, which shows some interesting phenomena. Under this set of parameters, if the diffusion of mosquitoes and humans is not considered, then the disease is extinct in urban areas, and is persistent in rural areas. The reason for this result is that we study spatial heterogeneity. Due to urban environment is not suitable for mosquitoes to survive, the number of mosquitoes is relatively small, and the medical resources are relatively large, so the urban disease control is better. Compared with the right column of Figure 6, the movement of mosquitoes and humans will reduce the risk of disease transmission to a certain extent, which is consistent with the conclusions in Figure 5(d) and reference [24]. There are two explanations: (i) The existence of diffusion indicates that the infected humans in rural areas can be treated in cities with better medical conditions. (ii) As mentioned in reference [24], to a certain extent, it is difficult for mosquitoes to take a blood meal among moving humans. Figure 6. The evolution of infection compartments of humans and mosquitoes with/without diffusion.

Discussion
As we all know, in the study of malaria transmission, the combined effects of temperature-dependent EIP and maturation delay, vector-bias, spatial structure and seasonal variation are worth studying.
Here, we establish and analyze a reaction-diffusion model of malaria, including periodic delays and vector-bias in a heterogeneous environment.
Using the theory proposed by Liang et al. [41] and Zhao [42], the basic reproduction number R 0 can be introduced. We first prove that in a defined phase space, linear system (4.2) generates a periodic semiflow which is eventually strongly monotonic. Then, we describe R 0 as a threshold for the extinction and persistence of malaria. It is proved that if R 0 < 1, then malaria will be eliminated. When R 0 > 1, there admits at least one positive ω-periodic solution.
In the numerical part, data for exploring the transmission of malaria in Maputo Provence, Mozambique have been published. The numerical simulation is divided into three parts. First, we give some simulation conclusions to support the theoretical analysis. Second, we discuss the sensitivity of R 0 to the parameters in the model. (I) Our numerical results show that if we only consider the time-average maturation delay, we can underestimate the risk of malaria transmission assessed by R 0 . (II) It is easy to know that R 0 is decreasing with respect to l/p. Thence, if the impact of vector-bias effect on malaria transmission is ignored, the risk of transmission will be underestimated. (III) We analyze the impact of intervention measures on R 0 from two aspects: (i) the impact of increasing medical resources. This case is mainly due to the uneven distribution of medical resources in rural and urban areas, which increases the impact of medical resources on R 0 ; (ii) for fixed medical resources, the impact of different allocation measures on disease transmission. Through the analysis of this case, we find that in the case of limited medical resources, narrowing the gap in the allocation of medical resources between rural and urban areas can better reduce the risk of malaria transmission. Third, we explore the impact of diffusion on R 0 in a spatially heterogeneous environment, assuming that the distribution of medical resources and the initial data of mosquitoes and humans depend on space. Biologically speaking, we understand this spatial heterogeneity as the difference between rural and urban areas. In addition, we also study the impact of diffusion on the final size of malaria. We find that under the set parameters, if the diffusion of mosquitoes and humans is not considered, the disease will become extinct in urban areas, and persist in rural areas. But under this set of parameters, the disease will become extinct if the diffusion is not considered.
Analyze the following time-periodic reaction-diffusion equation By reference [38, Lemma 2.1], system (A.2) has a unique globally attractive positive ω-periodic solution u * 1 (t, x) in Y + . Since system (A.2) can dominate the first equation of (2.12), there is a B 1 > 0 such that for any ϕ ∈ W h , a positive integer l 1 = l 1 (ϕ) > 0 exists and meets u 1 (t, x, ϕ) ≤ B 1 for all t ≥ l 1 ω and x ∈Ω.
Next, similar to reference [46, Theorem 2.1], the ultimate boundedness of solutions can be proved.