784
Views
0
CrossRef citations to date
0
Altmetric
Research Article

Dynamics of a multi-strain malaria model with diffusion in a periodic environment

, &
Pages 766-815 | Received 01 Jun 2022, Accepted 31 Oct 2022, Published online: 22 Nov 2022

Abstract

This paper mainly explores the complex impacts of spatial heterogeneity, vector-bias effect, multiple strains, temperature-dependent extrinsic incubation period (EIP) and seasonality on malaria transmission. We propose a multi-strain malaria transmission model with diffusion and periodic delays and define the reproduction numbers Ri and R^i (i = 1, 2). Quantitative analysis indicates that the disease-free ω-periodic solution is globally attractive when Ri<1, while if Ri>1>Rj (ij,i,j=1,2), then strain i persists and strain j dies out. More interestingly, when R1 and R2 are greater than 1, the competitive exclusion of the two strains also occurs. Additionally, in a heterogeneous environment, the coexistence conditions of the two strains are R^1>1 and R^2>1. Numerical simulations verify the analytical results and reveal that ignoring vector-bias effect or seasonality when studying malaria transmission will underestimate the risk of disease transmission.

1. Introduction

Malaria is a vector-borne disease that is spread from person to person owing to the bite of a female Anopheles mosquito and is epidemic in more than 100 countries [Citation32]. According to available information, there are 219 million cases and approximately 3.3 billion people in the world are at risk of contracting this disease [Citation38]. The spread of malaria directly threatens public health and has a huge negative shock on the local economy. Hence, it is critical to investigate the transmission of malaria in the population. Ross [Citation28] first introduced the mathematical model of malaria transmission and then Macdonald [Citation23] extended it, which provides useful insights into the transmission dynamics [Citation41]. Whereas the Ross–Macdonald model is exceedingly simplistic and ignores many key elements of real-world biology and epidemiology [Citation29]. For example, the distribution of rural and urban areas in society. Therefore, in epidemiological models of vector-borne, spatial heterogeneity must be considered [Citation11,Citation13]. People often use reaction–diffusion equation to understand the impact of the movements of human and mosquito populations in heterogeneous environments on disease transmission [Citation2,Citation21]. Furthermore, the following several important biological factors related to malaria transmission cannot be ignored.

A typical feature of malaria transmission is ‘vector-bias’ effect, which expounds the distinction in the probability of mosquitoes selecting humans. In 2005, Lacroix proved that mosquitoes are more likely to be attracted to infected persons, a phenomenon known as the vector-bias effect [Citation18]. To understand this preference, Ref. [Citation5] combined vector bias into the model, selecting humans based on the different probability that mosquitoes randomly reach humans and whether he is infected or not. Subsequently, people begin to show solicitude for this issue [Citation5,Citation34]. Assessing the effect of vector bias in diverse regions is of great significance for understanding the spread of malaria.

Trade-off mechanisms for coexistence of strains have been a long-standing question in the evolution of pathogens. It is well know that malaria is caused by Plasmodium parasites, and there exists five species of Plasmodium that can infect humans: (1) P. vivax, (2) P. falciparum, (3) P. malariae, (4) P. ovale, and (5) P. knowlesi. P. falciparum causes the highest mortality rate and is malignant, while the others are benign [Citation12]. In fact, there are many clinical cases of malaria caused by the other four species of Plasmodium. In the mathematical model with multiple strains, the main outcome is competitive exclusion [Citation3], but multiple mechanisms have been found to cause strain coexistence [Citation30,Citation46]. The WHO reports that both P. falciparum and P. vivax are endemic in countries, for example, Brazil, Cambodia, Guyana, Pakistan, Afghanistan, and Ethiopia, indicating that humans in these areas are at risk of infection with these two species [Citation37]. In fact, there are many clinical cases of malaria caused by the other four species of Plasmodium. Whereas we find that most of the existing studies only considered transmission of a single strain of malaria between humans and mosquitoes. Simultaneously, different malaria parasites have variable infection and recovery rates [Citation48]. It is necessary to study multi-strain malaria transmission models.

Temperature directly affects the survivability, lifespan, size, blood supply and reproductive capacity of mosquitoes [Citation7], which is the main vector of malaria. More interestingly, there is considerable evidence that the EIP is exceedingly sensitive to seasonal changing temperature [Citation9,Citation22]. EIP means that mosquitoes may not transmit disease to humans for some time after picking up infected blood. The longevity of a mosquito is generally less than 100 days, and the EIP can arrive 30 days [Citation36]. Understanding its role in malaria spread is critical on account of seasonal changing temperature [Citation26], but few papers on malaria transmission think out the seasonality. Especially the combined effect of the aforesaid elements on the spatial spread of malaria seems to receive little attention.

This paper modifies the Ross model [Citation28] to include vector-bias effect, two strains, temperature-dependent EIPs and diffusion in a heterogeneous space. This work is motivated by the following biological questions: (1) How do seasonal temperature changes and vector bias affect the spread of malaria? (2) Do mosquito and human movements have different effects on the spread of malaria in different regions? (3) Considering seasonal factors, can these two strains coexist in a heterogeneous space? If yes, what conditions need to be met? As these issues are strongly associated with the spread and control of malaria, further research is needed. The challenging point of this work is to explore the asymptotic behaviour of a complex model that considers multiple factors that influence malaria transmission. In this paper, we extend the theoretical analysis method in Refs. [Citation22,Citation40]. Considering the interaction and coexistence conditions of the two strains in a spatially heterogeneous environment is what distinguishes us from other papers. Since two strains are considered, the system will have multiple boundary periodic solutions, which is different from previous models with only one strain, which leads to the complexity of proofs, especially the proof of uniform persistence.

The remaining parts of the paper are arranged. In Section 2, we develop a reaction–diffusion model in a heterogeneous space that differs from the previously mentioned models. The meaning of the parameters in the model is also explained in this section. Section 3 is dedicated to the well-posedness. In Section 4, we introduce Ri, R0 and R^i (i=1,2), and explore the solution maps of correlated linear reaction–diffusion subsystems with periodic delay. Section 5 establishes threshold dynamics based on Ri and R^i. Section 6 verifies the previous theoretical results through simulations. In addition, the impact of the diffusion of humans and mosquitoes on the prevention and control of malaria transmission is explored and, how to adjust the allocation of medical resources is analysed. The last is a brief conclusion.

2. The model

The objective of this section is to establish a model of malaria transmission in a spatially heterogeneous environment, incorporating temperature-dependent delays, vector-bias effect and two strains. One hypothesis is that susceptible individuals and mosquitoes can only be infected by one strain and that individuals become susceptible again after recovering, but the mosquitoes cannot recover due to their relatively short lifespan. Let Nh(t,x) be the total population consisting of three classes: susceptible population (Sh(t,x)), individuals infected with strain 1 (I1(t,x)), and individuals infected with strain 2 (I2(t,x)), where t denotes time and x stands for position. Then for t>0, Nh(t,x)=Sh(t,x)+I1(t,x)+I2(t,x), obeys (1) Nh(t,x)t=DhΔNh(t,x)+bh(x)μh(x)Nh(t,x),(1) and xΩ, where Ω is a bounded domain with smooth boundary Ω; Dh means the diffusion coefficient of human, which is Hölder continuous and positive on Ω¯; Δ is the usual Laplacian operator. Assume that the form of diffusion is random, it means that individual walkers walk randomly on the solid line using a fixed step [Citation4]. μh(x) and bh(x) are the natural death rate and the recruitment rate of human at position x, respectively. The Neumann boundary condition means that no population flux crosses the boundary Ω, (2) Nh(t,x)ν=0,t>0,xΩ,(2)

where ν is the normal derivative along the outward unit normal vector ν on Ω. We know that system (Equation1) with (Equation2) has a positive steady state N(x)C(Ω¯,R+){0} under appropriate assumptions which are globally attractive [Citation45, Theorem 3.15 and Theorem 3.16], where R+ is the set of positive real numbers and C(Ω¯,R+) marks the Banach space of continuous functions from Ω¯ to R+. For the sake of simplicity, the total number of human at time t and location x stabilizes at N(x), namely, t0, xΩ, Nh(t,x)N(x), which is motivated by Refs. [Citation2,Citation21,Citation35]. Then Sh(t,x)=N(x)I1(t,x)I2(t,x).

Similar to the discussion in Ref. [Citation27], temperature-dependent EIPs are introduced. Let Sv(t,x), Evi(t,x), and Ivi(t,x) (i = 1, 2) be the number of susceptible, exposed and infected mosquitoes, respectively. Here, the subscript i indicates strain i (i=1,2). For t>0, we set M(t,x)=Sv(t,x)+Ev1(t,x)+Ev2(t,x)+Iv1(t,x)+Iv2(t,x) to be the overall number of mosquitoes, and to obey (3) {M(t,x)t=DvΔM(t,x)+Λ(t,x)μv(t,x)M(t,x),xΩ,M(t,x)ν=0,xΩ,(3) where Λ(t,x)0 and μv(t,x) denote the recruitment and the natural mortality rates of mosquitoes, respectively. Particularly, based on Ref. [Citation40], we concentrate on the case where this disease has no discernible effect on the movement of humans and mosquitoes. Mathematically, it is assumed that all individuals have the identical diffusion coefficient Dh>0, and all mosquitoes have the same diffusion coefficient Dv>0. Mosquito biting rate β(t,x) refers to the number of bites of a mosquito per unit time at t and x. Furthermore, Λ(t,x), μv(t,x) and β(t,x) are Hölder continuous and nonnegative nontrivial on R×Ω¯ and is ω-periodic in t. Following the approach in Refs. [Citation2,Citation5,Citation36], at time t and location x, the numbers of newly infectious humans and newly infected mosquitoes per unit, are ciβ(t,x)l[N(x)I1(t,x)I2(t,x)]Ivi(t,x)p[I1(t,x)+I2(t,x)]+l[N(x)I1(t,x)I2(t,x)]andαiβ(t,x)pIi(t,x)Sv(t,x)p[I1(t,x)+I2(t,x)]+l[N(x)I1(t,x)I2(t,x)],i = 1, 2, respectively, where p and l are defined as the probability that mosquitoes randomly arrive at humans and choose infect and susceptible humans, respectively [Citation5]. Here ci is the transmission probability from mosquitoes infected with strain i to susceptible humans per bite, and αi is the transmission probability from humans infected with strain i to susceptible mosquitoes per bite (i = 1, 2). The incubation period of humans is relatively short (about 20 days), but the infection period can last for several months, and the lifespan of humans is decades, so the incubation period of humans can be ignored here.

On the basis of the above discussions, one has (4) {I1(t,x)t=DhΔI1(t,x)+c1β(t,x)l[N(x)I1(t,x)I2(t,x)]Iv1(t,x)p[I1(t,x)+I2(t,x)]+l[N(x)I1(t,x)I2(t,x)](μh(x)+ϱ1(x))I1(t,x),xΩ,I2(t,x)t=DhΔI2(t,x)+c2β(t,x)l[N(x)I1(t,x)I2(t,x)]Iv2(t,x)p[I1(t,x)+I2(t,x)]+l[N(x)I1(t,x)I2(t,x)](μh(x)+ϱ2(x))I2(t,x),xΩ,Sv(t,x)t=DvΔSv(t,x)+Λ(t,x)α1β(t,x)pI1(t,x)Sv(t,x)p[I1(t,x)+I2(t,x)]+l[N(x)I1(t,x)I2(t,x)]α2β(t,x)pI2(t,x)Sv(t,x)p[I1(t,x)+I2(t,x)]+l[N(x)I1(t,x)I2(t,x)]μv(t,x)Sv(t,x),xΩ,Ev1(t,x)t=DvΔEv1(t,x)+α1β(t,x)pI1(t,x)Sv(t,x)p[I1(t,x)+I2(t,x)]+l[N(x)I1(t,x)I2(t,x)]Mv1(t,x)μv(t,x)Ev1(t,x),xΩ,Ev2(t,x)t=DvΔEv2(t,x)+α2β(t,x)pI2(t,x)Sv(t,x)p[I1(t,x)+I2(t,x)]+l[N(x)I1(t,x)I2(t,x)]Mv2(t,x)μv(t,x)Ev2(t,x),xΩ,Iv1(t,x)t=DvΔIv1(t,x)+Mv1(t,x)μv(t,x)Iv1(t,x),xΩ,Iv2(t,x)t=DvΔIv2(t,x)+Mv2(t,x)μv(t,x)Iv2(t,x),xΩ,I1(t,x)ν=I2(t,x)ν=Sv(t,x)ν=Iv1(t,x)ν=Iv2(t,x)ν=0,xΩ,(4) for t>0, where Mvi(t,x) measures the number of newly occurred infected mosquitoes with strain i (i=1,2). Further, ϱi(x) represents the recovery rate of strain i (i = 1, 2), which is Hölder continuous and positive on Ω¯.

The expression of Mvi is derived. Let τi(t) denote the length of temperature-dependent EIP infected by strain i (i = 1, 2), because it is assumed that the temperature T varies with time t. Set q as the development level of infection, it is easy to see that q describes the completeness of the parasite during the developmental stage of the mosquito (in other words, the completeness of the incubation period), and γ(t) is a temperature-dependent rate increase in q, i.e. dqdt=γ(t). Let ρ(t,q,x) be the density of mosquitoes infected by strain 1, with the infection development level q at time t and position x.

According to the discussions in Appendix 1 and Refs. [Citation17,Citation40], we can derive the expression of Mvi (i = 1, 2) is Mvi(t,x)=(1τi(t))ΩΓ(t,tτi(t),x,y)αiβ(tτi(t),y)pIi(tτi(t),y)Sv(tτi(t),y)p[I1(tτi(t),y)+I2(tτi(t),y)]+l[N(y)I1(tτi(t),y)I2(tτi(t),y)]dy.Γ(t,t0,x,y) represents Green function related to u(t,x)t=DvΔu(t,x)μv(t,x)u(t,x) obeys the Neumann boundary condition. Γ(t,t0,x,y) is also biologically meaningful and represents the probability that a mosquito at position y at time t0 will arrive at position x after time tt0. Emphasise that Γ(t+ω,t0+ω,x,y)=Γ(t,t0,x,y) for t>t00 and x,yΩ by virtue of μv(t,)=μv(t+ω,). Here, we assume that the natural death rate of mosquitoes infected with two strains is the same.

Substituting Mvi(t,x) (i=1,2) into system (Equation4), one has

for t>0, where all constant parameters are positive and τi(t) is positive, continuous and ω-periodic function in t. Based on Ref. [Citation40], we get that 1τi(t)=γ(t)γ(tτi(t))>0, (i=1,2). See Table  for the biological interpretations.

Since Evi(t,x) (i=1,2) is decoupled from the other equations in system (5), the following system should be adequate, (6) {u1(t,x)t=DhΔu1(t,x)+c1β(t,x)l[N(x)u1(t,x)u2(t,x)]u4(t,x)p[u1(t,x)+u2(t,x)]+l[N(x)u1(t,x)u2(t,x)](μh(x)+ϱ1(x))u1(t,x),xΩ,u2(t,x)t=DhΔu2(t,x)+c2β(t,x)l[N(x)u1(t,x)u2(t,x)]u5(t,x)p[u1(t,x)+u2(t,x)]+l[N(x)u1(t,x)u2(t,x)](μh(x)+ϱ2(x))u2(t,x),xΩ,u3(t,x)t=DvΔu3(t,x)+Λ(t,x)α1β(t,x)pu1(t,x)u3(t,x)p[u1(t,x)+u2(t,x)]+l[N(x)u1(t,x)u2(t,x)]α2β(t,x)pu2(t,x)u3(t,x)p[u1(t,x)+u2(t,x)]+l[N(x)u1(t,x)u2(t,x)]μv(t,x)u3(t,x),xΩ,u4(t,x)t=DvΔu4(t,x)μv(t,x)u4(t,x)+(1τ1(t))ΩΓ(t,tτ1(t),x,y)α1β(tτ1(t),y)pu1(tτ1(t),y)u3(tτ1(t),y)p[u1(tτ1(t),y)+u2(tτ1(t),y)]+l[N(y)u1(tτ1(t),y)u2(tτ1(t),y)]dy,xΩ,u5(t,x)t=DvΔu5μv(t,x)u5(t,x)+(1τ2(t))ΩΓ(t,tτ2(t),x,y)α2β(tτ2(t),y)pu2(tτ2(t),y)u3(tτ2(t),y)p[u1(tτ2(t),y)+u2(tτ2(t),y)]+l[N(y)u1(tτ2(t),y)u2(tτ2(t),y)]dy,xΩ,u1(t,x)ν=u2(t,x)ν=u3(t,x)ν=u4(t,x)ν=u5(t,x)ν=0,xΩ,(6) for t>0, where (u1(t,x),u2(t,x),u3(t,x),u4(t,x),u5(t,x))(I1(t,x),I2(t,x),Sv(t,x),Iv1(t,x),Iv2(t,x)).

3. Well-posedness

Let X1:=C(Ω¯,R5) be the Banach space of continuous functions from Ω¯ to R5 with the supremum norm X1, and X1+:=C(Ω¯,R+5). Let τ^1:=maxt[0,ω]τ1(t), τ^2:=maxt[0,ω]τ2(t) and τ^:=max{τ^1,τ^2}. Define X1:=C([τ^,0],X1) to be a Banach space, with the norm ϕ:=maxθ[τ^,0]ϕ(θ)X1, ϕX1, and X1+:=C([τ^,0],X1+). Then (X1,X1+) and (X1,X1+) are ordered spaces. Considering a function z:[τ^,0)X1 for σ>0, we define ztX1 by zt(θ)=(z1(t+θ),z2(t+θ),z3(t+θ),z4(t+θ),z5(t+θ)),θ[τ^,0],for any t[0,σ).

Set Y:=C(Ω¯,R), Y+:=C(Ω¯,R+), and Yh:={φY+:0φ(x)N(x),xΩ¯}. Presume that Ti(t,s):YY (i=1,2) is the evolution operator related to ui(t,x)t=DhΔui(t,x)(μh(x)+ϱi(x))ui(t,x):=Aiui(t,x),and T3(t,s):YY is the evolution operator related to u3(t,x)t=DvΔu3(t,x)μv(t,x)u3(t,x):=A3(t)u3(t,x),obeys the Neumann boundary condition, respectively. Noting that Tj(t,s)=Tj(ts) and μv(t,) is ω-periodic in t, then Tj(t+ω,s+ω)=Tj(t,s) (j=1,2,3) for any (t,s)R2 with ts [Citation16, Lemma 6.1]. Moreover, from Section 7.1 and Corollary 7.2.3 in Ref. [Citation31], we know that Tj(t,s) is compact and strongly positive for (t,s)R2 with t>s. Set T(t,s)=diag{T1(t,s),T2(t,s),T3(t,s),T3(t,s),T3(t,s)}:X1X1, t0, to be a strongly continuous semigroup, and A(t)=diag{A1,A2,A3(t),A3(t),A3(t)}. Define Ai D(Ai):={φ~C2(Ω¯):φ~ν=0onΩ},Aiφ~=DhΔφ~(μh(x)+ϱi(x))φ~,φ~D(Ai),i=1,2,and A3(t) is defined by D(A3(t)):={φ~C2(Ω¯):φ~ν=0onΩ},A3(t)φ~=DvΔφ~μv(t,x)φ~,φ~D(A3(t)).Then T(t,s):X1X1 is a semigroup generated by the operator A(t) defined on D(A(t))=D(A1)×D(A2)×D(A3(t))×D(A3(t))×D(A3(t)). Set Wh:={ϕX1+:0ϕ1(x)+ϕ2(x)N(x),xΩ¯}, and Wh:=C([τ^,0],Wh). We define F=(F1,F2,F3,F4,F5):[0,+)×WhX1 by {F1(t,ϕ)=c1β(t,)l[N()ϕ1(0,)ϕ2(0,)]ϕ4(0,)p[ϕ1(0,)+ϕ2(0,)]+l[N()ϕ1(0,)ϕ2(0,)],F2(t,ϕ)=c2β(t,)l[N()ϕ1(0,)ϕ2(0,)]ϕ5(0,)p[ϕ1(0,)+ϕ2(0,)]+l[N()ϕ1(0,)ϕ2(0,)],F3(t,ϕ)=Λ(t,)α1β(t,)pϕ1(0,)ϕ3(0,)p[ϕ1(0,)+ϕ2(0,)]+l[N()ϕ1(0,)ϕ2(0,)]α2β(t,)pϕ2(0,)ϕ3(0,)p[ϕ1(0,)+ϕ2(0,)]+l[N()ϕ1(0,)ϕ1(0,)],F4(t,ϕ)=(1τ1(t))ΩΓ(t,tτ1(t),,y)α1β(tτ1(t),y)pϕ1(τ1(t),y)ϕ3(τ1(t),y)p[ϕ1(τ1(t),y)+ϕ2(τ1(t),y)]+l[N(y)ϕ1(τ1(t),y)ϕ2(τ1(t),y)]dy,F5(t,ϕ)=(1τ2(t))ΩΓ(t,tτ2(t),,y)α2β(tτ2(t),y)pϕ2(τ2(t),y)ϕ3(τ2(t),y)p[ϕ1(τ2(t),y)+ϕ2(τ2(t),y)]+l[N(y)ϕ1(τ2(t),y)ϕ2(τ2(t),y)]dy,for ϕ=(ϕ1,ϕ2,ϕ3,ϕ4,ϕ5)TWh, t0 and xΩ¯. As a consequence, system (Equation6) can be rewritten as the following abstract functional differential equation {u(t,x)t=A(t)u(t,x)+F(t,ut),u(θ,x)=ϕ(θ,x),where u(t,)=(u1(t,),u2(t,),u3(t,),u4(t,),u5(t,)), t>0, xΩ and θ[τ^,0]. Owing to the ω-periodicity of β(t,), τi(t), Λ(t,), τi(t) and μv(t,) (i = 1, 2), we know that A(t) and F(t,ut) are ω-periodic in t. Then, from Corollary 4 in Ref. [Citation25] and Corollary 8.1.3 in Ref. [Citation39], for each ϕWh, a mild solution can be obtained as a continuous solution of the following integral equation u(t,ϕ)=T(t,0)ϕ(0)+0tT(ts)F(s,us)ds,t>0,u0=ϕWh.Using Corollary 7.3.2 of Ref. [Citation31], we have the following proposition on the unique global solution that exists in the system (Equation6).

Lemma 3.1

For ϕWh, system (Equation6) admits the unique solution, remarked as z(t,,ϕ) on its maximal existence interval [0,t¯ϕ) with z0=ϕ, where t¯ϕ. Additionally, for t[0,t¯ϕ), z(t,,ϕ)Yh×Yh×Y+×Y+×Y+ and for all t>τ^, z(t,,ϕ) is a classical solution to system (Equation6).

We first return to system (5) for more observations. Considering the biological significance of τi(t) (i=1,2), the following compatibility condition are imposed (7) Evi(0,)=τi(0)0T3(0,s)αiβ(s,)pIi(s,)Sv(s,)p[I1(s,)+I2(s,)]+l[N()I1(s,)I2(s,)]ds.(7) Define D:={τ2(0)0ψ~C([τ^,0],C(Ω¯,R+7)):ψ~1(θ,)+ψ~2(θ,)N(),θ[τ^,0],ψ~4(0,)=τ1(0)0T3(0,s)α1β(s,)pψ~1(s,)ψ~3(s,)p[ψ~1(s,)+ψ~2(s,)]+l[N()ψ~1(s,)ψ~2(s,)]ds,ψ~5(0,)=τ2(0)0T3(0,s)α2β(s,)pψ~2(s,)ψ~3(s,)p[ψ~1(s,)+ψ~2(s,)]+l[N()ψ~1(s,)ψ~2(s,)]ds}.Therefore, for any ψ~D, system (5) admits a unique solution U(t,,ψ~)=(I1(t,),I2(t,),Sv(t,),Ev1(t,),Ev2(t,), Iv1(t,),Iv2(t,)) satisfying U0=ψ~. From [Citation25, Corollary 4], we get Ii(t,)0, Sv(t,)0 and Ivi(t,)0 (i=1,2). According to the uniqueness of solution and the compatibility conditions (Equation7), it follows that (8) Evi(t,)=tτi(t)tT3(t,s)αiβ(s,)pIi(s,)Sv(s,)p[I1(s,)+I2(s,)]+l[N()I1(s,)I2(s,)]ds,i=1,2.(8) Hence, Evi(t,)0.

Next, we proof the ultimate boundedness of the solution of system (Equation6).

Lemma 3.2

There is a G>0, such that any solution (u1(t,x),u2(t,x),u3(t,x),u4(t,x),u5(t,x)) of (Equation6) meets (9) 0ui(t,x)G,i=1,2,3,4,5,t>0,xΩ¯.(9)

Proof.

Based on (Equation1), one has Nh(t,x)tDhΔNh(t,x)+b~hμ¯hNh(t,x),where b~h=maxxΩ¯bh(x) and μ¯h=minxΩ¯μh(x). Similarly, for (Equation3), we also have M(t,x)tDvΔM(t,x)+Λ~μ¯vM(t,x), where Λ~=maxt[0,ω],xΩ¯Λ(t,x) and μ¯v=mint[0,ω],xΩ¯μv(t,x). Choose G=max{b~hμ¯h,Λ~μ¯v}, therefore, (Equation9) holds.

Based on the above discussion, the solution of (Equation6) in Wh exists globally on [0,) and is ultimately bounded. Through similar arguments with Lemma 2.6 in [Citation15], Lemma 2.1 in [Citation43], Theorem 2.1 and Theorem 7.3.1 of [Citation31], together with Theorem 2.9 in [Citation24], we have the existence of continuous semi-flow.

Lemma 3.3

For each ϕWh, system (Equation6) admits a unique solution u(t,,ϕ) on [0,+) with u0=ϕ. Moreover, system (Equation6) generates an ω-periodic semi-flow Φ~t:WhWh, defined by Φ~t:=ut(ϕ),t0, and Φ~:=Φ~ω has a strong global attractor in Wh.

4. Reproduction numbers

An extremely important concept in epidemiology is basic reproduction number, which is generally defined as the average number of secondary infections produced by a type infected human when join a completely susceptible population during the entire infection period [Citation8]. Furthermore, we also introduce invasion reproduction numbers, the average number of secondary infections in a population where an infected individual is susceptible to this strain, but another strain is already an endemic disease, which together with the basic reproduction numbers define the threshold behaviours for the epidemic model.

4.1. Basic reproduction numbers

We first explore two subsystems involving merely one strain, for t>0, (10) {u¯1(t,x)t=DhΔu¯1(t,x)+c1β(t,x)l[N(x)u¯1(t,x)]u¯3(t,x)pu¯1(t,x)+l[N(x)u¯1(t,x)](μh(x)+ϱ1(x))u¯1(t,x),xΩ,u¯2(t,x)t=DvΔu¯2(t,x)+Λ(t,x)α1β(t,x)pu¯1(t,x)u¯2(t,x)pu¯1(t,x)+l[N(x)u¯1(t,x)]μv(t,x)u¯2(t,x),xΩ,u¯3(t,x)t=DvΔu¯3(t,x)μv(t,x)u¯3(t,x)+(1τ1(t))ΩΓ(t,tτ1(t),x,y)α1β(tτ1(t),y)pu¯1(tτ1(t),y)u¯2(tτ1(t),y)pu¯1(tτ1(t),y)+l[N(y)u¯1(tτ1(t),y)]dy,xΩ,u¯1(t,x)ν=u¯2(t,x)ν=u¯3(t,x)ν=0,xΩ,(10) which is the subsystem of strain 1, where (u¯1(t,x),u¯2(t,x),u¯3(t,x))(u1(t,x),u3(t,x),u4(t,x)), and for t>0, (11) {u^1(t,x)t=DhΔu^1(t,x)+c2β(t,x)l[N(x)u^1(t,x)]u^3(t,x)pu^1(t,x)+l[N(x)u^1(t,x)](μh(x)+ϱ2(x))u^1(t,x),xΩ,u^2(t,x)t=DvΔu^2(t,x)+Λ(t,x)α2β(t,x)pu^1(t,x)u^2(t,x)pu^2(t,x)+l[N(x)u^1(t,x)]μv(t,x)u^2(t,x),xΩ,u^3(t,x)t=DvΔu^3(t,x)μv(t,x)u^3(t,x)+(1τ2(t))ΩΓ(t,tτ2(t),x,y)α2β(tτ2(t),y)pu^1(tτ2(t),y)u^2(tτ2(t),y)pu^1(tτ2(t),y)+l[N(y)u^1(tτ2(t),y)]dy,xΩ,u^1(t,x)ν=u^2(t,x)ν=u^3(t,x)ν=0,xΩ,(11) which is the subsystem of strain 2, where (u^1(t,x),u^2(t,x),u^3(t,x))(u2(t,x),u3(t,x),u5(t,x)). Notice that the positive ω-periodic solution or disease-free ω-periodic solution for each subsystem can be used as a boundary periodic solution of system (Equation6). According to Ref. [Citation10], E0=(0,0,u3(t,),0,0), E01=(0,u3(t,),0) and E02=(0,u3(t,),0) denote the disease-free ω-periodic solution of system (Equation6), (Equation10) and (Equation11), respectively. E1=(u¯1(t,),u¯3(t,),u¯4(t,)), is the positive ω-periodic solution of (Equation10) (and the boundary periodic solution of system (Equation6)). When E1 is the nontrivial boundary periodic solution of system (Equation6), we call it the periodic solution of strain 1, represented by (u¯1(t,),0,u¯3(t,),u¯4(t,),0). E2=(u~2(t,),u~3(t,),u~5(t,)) represents the positive ω-periodic solution of (Equation11) (and the boundary periodic solution of system (Equation6)). When E2 is the nontrivial boundary periodic solution of system (Equation6), we call it the periodic solution of strain 2, denoted by (0,u~2(t,),u~3(t,),0,u~5(t,)).

Set Ri to be basic reproduction number of strain i without other types of strains (i=1,2). The parameters are space-dependent and time-periodic, leading to the complexity of model analysis. We will make use of the results in Refs. [Citation19,Citation44] to define basic reproduction numbers R1 and R2 of subsystems (Equation10) and (Equation11). Let strain 1 as an exemplar, and the outcomes also apply to strain 2. The next generation generator of subsystem (Equation10) is defined, whose spectral radius is R1.

Set u¯1(t,x)=u¯3(t,x)=0 in system (Equation10), for t>0, we get (12) {u¯2(t,x)t=DvΔu¯2(t,x)+Λ(t,x)μv(t,x)u¯2(t,x),xΩ,u¯2(t,x)ν=0,xΩ.(12) According to [Citation43, Lemma 2.1], we can get the following lemma.

Lemma 4.1

System (Equation12) has a unique globally attractive positive ω-periodic solution u3(t,) in Y+.

Linearizing system (Equation10) at E01, and for t>0, considering the equations (13) {v1(t,x)t=DhΔv1(t,x)+c1β(t,x)v2(t,x)(μh(x)+ϱ1(x))v1(t,x),xΩ,v2(t,x)t=DvΔv2(t,x)μv(t,x)v2(t,x)+(1τ1(t))ΩΓ(t,tτ1(t),x,y)α1β(tτ1(t),y)pu3(tτ1(t),y)v1(tτ1(t),y)lN(y)dy,xΩ,v1(t,x)ν=v2(t,x)ν=0,xΩ.(13) Let v¯(v1(t,x),v2(t,x))=(u¯1(t,x),u¯3(t,x)). Similarly, linearizing (Equation11) around E02, and considering the equations of infective compartments (14) {v3(t,x)t=DhΔv3(t,x)+c2β(t,x)v4(t,x)(μh(x)+ϱ2(x))v3(t,x),t>0,xΩ,v4(t,x)t=DvΔv4(t,x)μv(t,x)v4(t,x)+(1τ2(t))ΩΓ(t,tτ2(t),x,y)α2β(tτ2(t),y)pu3(tτ2(t),y)v3(tτ2(t),y)lN(y)dy,t>0,xΩ,v3(t,x)ν=v4(t,x)ν=0,t>0,xΩ.(14) Then v^(v3(t,x),v4(t,x))=(u^1(t,x),u^3(t,x)).

Let X2:=C(Ω¯,R2), X2+:=C(Ω¯,R+2), and Cω(R,X2) be the Banach space consisting of all ω-periodic and continuous functions from R to X2, where ψCω(R,X2):=maxθ[0,ω]ψ(θ)X2 for any ψCω(R,X2). Denote X2:=C([τ^,0],X2) and X2+:=C([τ^,0],X2+). Define F1(t):X2X2 by F1(t)(ϕ1ϕ4)=(c1β(t,)ϕ4(0,)(1τ1(t))ΩΓ(t,tτ1(t),,y)α1β(tτ1(t),y)pu3(tτ1(t),y)ϕ1(τ1(t),y)lN(y)dy),for tR, (ϕ1,ϕ4)X2, and V1(t)v¯=D~Δv¯W1(t)v¯, where D~=diag(Dh,Dv) and [W1(t)](x)=((μh(x)+ϱ1(x))00μv(t,x)),xΩ¯.System (Equation14) can be marked as dv¯dt=F1(t)v¯tV1(t)v¯,t0.Set Ψ1(t,s)=diag(T1(t,s),T3(t,s)), ts, to be the evolution operator, related to dv¯dt=V1(t)v¯,where T1(t,s) and T3(t,s) are defined in Section 3, and subject to the Neumann boundary condition. Then [Citation33, Theorem 3.12] implies that V1(t) is resolvent positive.

Define the exponential growth bound of Ψ1(t,s) as ω¯(Ψ1)=inf{ω~:L1suchthatΨ1(t+s,s)Leω~t,sR,t0}.By [Citation33, Proposition A.2], we have ω¯(Ψ1)=lnr(Ψ1(ω,0))ω=lnr(Ψ1(ω+s,s))ω,s[0,ω].Refer to Krein–Rutman theorem and [Citation14, Lemma 14.2], 0<r(Ψ1(ω,0))=max{r(T1(ω,0)),r(T3(ω,0))}<1,where r(Ψ1(ω,0)) is the spectral radius of Ψ1(ω,0). By [Citation33, Proposition 5.6] with s = 0, we obtain ω¯(Ψ1)<0. Emphasise that Ψ1(t,s) is a positive operator, in that sense, Ψ1(t,s)X2+X2+ for all ts. It is easy to see that F1(t) maps X2+ into X2+, for any t0, W1(t) is cooperative and ω¯(Ψ1)<0.

Using [Citation19,Citation44], we can define R1 for subsystem (Equation10). Suppose that humans and mosquitoes are near (0,u3(t,x),0). Let v¯Cω(R,X2) and v¯(t) be the initial distribution of infectious humans and mosquitoes with strain 1 at tR. Note that for s0, F1(ts)v¯ts is the density distribution of newly infected humans and mosquitoes at ts, which is caused by the infected humans and mosquitoes who were introduced over the interval [tsτ^,ts]. Thus Ψ1(t,ts)F1(ts)v¯ts denotes the distribution of those infected humans and mosquitoes who newly became infectious at time ts and still remain in the infectious compartments at time t. Thus the integral 0+Ψ1(t,ts)F1(ts)v¯tsds=0+Ψ1(t,ts)F1(ts)v¯(ts+)ds represents the distribution of cumulative new infected humans and mosquitoes at time t caused by all those infectious introduced at all prior time to t. Define the next generation operator L1 associated with subsystem (Equation10) as [L1v¯](t):=0+Ψ1(t,ts)F1(ts)v¯(ts+)ds,tR,v¯Cω(R,X2).Then L1 is a positive and continuous operator, which maps v¯(t) to the distribution of the total infected members generated among infection period.

Motivated by Refs. [Citation8,Citation33], the spectral radius of L1 is defined as the basic reproduction number of (Equation10), (15) R1:=r(L1).(15) R2 is the same. Define the basic reproduction number of system (Equation6) as R0=max{R1,R2}.

For any ψ in X2, let P1t be the solution map of (Equation13) on X2, i.e. P1t(ψ)=v¯t(ψ), t0, where v¯t(ψ)(θ)(x)=v¯(t+θ,x,ψ)=(v1(t+θ,x,ψ),v2(t,x,ψ)), θ[τ^,0], and v¯(t,x,ψ) is the unique solution of system (Equation13) with v¯(θ,x)=ψ(θ,x) for all θ[τ^,0], xΩ¯. Hence, P1:=P1ω is the Poincaré map related to subsystem (Equation13). Denote P2t to be the solution map of the linear periodic Equation (Equation14) on X2, it means P2t(ϕ~)=v^t(ϕ~), t0, and v^t(ϕ~)(θ)(x)=v^(t+θ,x,ϕ~)=(v3(t+θ,x,ϕ~),v4(t,x,ϕ~)), θ[τ^2,0] with v^(θ,x)=ϕ~(θ,x) for all θ[τ^,0], and ϕ~ in X2. P2:=P2ω is the Poincaré map associated with subsystem (Equation14). r(Pi) denotes the spectral radius of Pi, (i=1,2).

In light of Ref. [Citation19, Theorem 3.7], we have the following results.

Lemma 4.2

Ri1 and r(Pi)1 (i = 1, 2) have the same sign.

Define X3:=C(Ω¯,R3), X3+:=C(Ω¯,R+3), then (X3,X3+) is an order Banach space. Set X3:=C([τ^,0],Yh×Y+×Y+). Let Qˇt(ϕ¯):=u¯t(ϕ¯), t0, is an ω-periodic semi-flow of system (Equation10) in X3, and Qˇ:=Qˇω.

Lemma 4.3

For ϕ¯X3, system (Equation10) has a unique solution u¯(t,,ϕ¯) with u¯0=ϕ¯. Additionally, system (Equation10) begets an ω-periodic semi-flow Qˇt(ϕ¯):=u¯t(ϕ¯) in X3, t0, and Qˇ:=Qˇω has a strong attractor in X3.

In order to investigate the dynamics behaviour of system (Equation6), first explore the system (Equation13) to generates a periodic semi-flow on the phase space E1 that is eventually strongly monotone, E1:=C([τ1(0),0],Y)×Y.Set E1+:=C([τ1(0),0],Y+)×Y+, E2:=C([τ2(0),0],Y)×Y, and E2+:=C([τ2(0),0],Y+)×Y+. Thus, (E1,E1+) and (E2,E2+) are ordered Banach spaces. For φE1, let ϖ(t,x,φ)=(ϖ1(t,x,φ),ϖ2(t,x,φ)) be the unique solution of (Equation13) with ϖ0(φ)(θ,x)=φ(θ,x) for all θ[τ1(0),0], xΩ¯, where (16) ϖt(φ)(θ,x)=ϖ(t+θ,x,φ)=(ϖ1(t+θ,x,φ),ϖ2(t,x,φ)),t0,(θ,x)[τ1(0),0]×Ω¯.(16)

Lemma 4.4

For each φE1+, system (Equation13) has a unique nonnegative solution ϖ(t,,φ) on [0,+) with ϖ0=φ.

Proof.

Set τ¯1=mint[0,ω]τ1(t). For t[0,τ¯1], since tτ1(t) is strictly increasing in t, one has τ1(0)=0τ1(0)tτ1(t)τ¯1τ1(τ¯1)τ¯1τ¯1=0.Thus, u¯1(tτ1(t),)=φ1(tτ1(t),). Therefore, for any t[0,τ¯1], system (Equation13) becomes {u¯1(t,x)t=DhΔu¯1(t,x)+c1β(t,x)u¯3(t,x)(μh(x)+ϱ1(x))u¯1(t,x),xΩ,u¯3(t,x)t=DvΔu¯3(t,x)μv(t,x)u¯3(t,x)+(1τ1(t))ΩΓ(t,tτ1(t),x,y)α1β(tτ1(t),y)pu3(tτ1(t),y)φ1(tτ1(t),y)lN(y)dy,xΩ,u¯1(t,x)ν=u¯3(t,x)ν=0,xΩ.Fix φE1+, the solution (u¯1(t,),u¯3(t,)) of the above system uniquely exists for t[0,τ¯1]. In other words, u¯1(θ,)=ϖ1(θ,) for θ[τ1(0),τ¯1] and u¯2(θ,)=ϖ2(θ,), for θ[0,τ¯1].

By repeating the above step to t[nτ¯1,(n+1)τ¯1], n=1,2,3, we see that the solution ϖ(t,,φ) with initial date φE1+ exists uniquely for all t0.

Remark 4.1

The uniqueness of solutions stated in Lemmas 4.3 and 4.4, consequently, for ψX2+ and φE1+ with ψ1(θ,)=φ1(θ,), θ[τ1(0),0] and ψ2()=φ2(), then v(t,,ψ)=ϖ(t,,φ), t0, where v(t,,ψ) and ϖ(t,,φ) are solutions of system (Equation13) meeting v0=ψ and ϖ0=φ, respectively.

Set P¯1t to be the solution map of (Equation13) on space E1, then P¯1:=P¯1ω is corresponding Poincaré map, and r(P¯1) is the spectral radius of P¯1. The next lemma indicates that P¯1t is eventually strongly monotone motivated by Ref. [Citation40] and [Citation20, Lemma 3.5 and Lemma 3.6].

Lemma 4.5

For each φE1+ with φ0, the solution ϖ(t,) of (Equation13) with ϖ0=φ satisfies ϖ(t,)>0 for all t>2τ^1, therefore, t>3τ^1, P¯1tφ0.

The proof of Lemma 4.5 is given in Appendix 2.

Fixed an integer n0 satisfying n0ω>3τ^1. Based on Lemma 4.5, one can easily see that P¯1n0 is strongly monotone. In addition, through similar arguments in [Citation15, Lemma 2.6], one can prove that P¯1n0 is compact. Apply the Krein–Rutman theorem to P¯1n0, and r(P¯1n0)=(r(P¯1))n0, one has λ=r(P¯1)>0, where λ is the simple eigenvalue of P¯1, and has a strongly positive eigenfunction ϑint(E+). Denote P¯2t to be the solution map of system (Equation14) on space E2, and let P¯2:=P¯2ω be the Poincaré map related to (Equation14). r(P¯2) is the spectral radius of P¯2. According to [Citation22, Lemma 3.8], one has the following lemma.

Lemma 4.6

Two Poincaré maps Pi:X2X2 and P¯i:=EiEi have the same spectral radius, i.e. r(Pi)=r(P¯i), (i=1,2).

To explore long-term dynamics, we give the fact that the system has a special solution. The next argument is motivated by the treatment in [Citation2, Lemma 5] and [Citation42, Proposition 1.1].

Lemma 4.7

There exists a v¯(t,x), which is positive ω-periodic function, then eμ1tv¯(t,x) is a solution of (Equation13), where μ1=lnr(P¯1)ω.

4.2. Invasion reproduction numbers

Use R^i to denote the invasion reproduction number, which means the ability of strain i to intrude another strain (i=1,2). When E2 is regarded as the non-trivial boundary periodic solution of (Equation6), u~2(t,x)>0 and u~5(t,x)>0 are non-zero. The periodic solution of strain 2 meets (17) {0=DhΔu~2(t,x)+c2β(t,x)l[N(x)u~2(t,x)]u~5(t,x)pu~2(t,x)+l[N(x)u~2(t,x)](μh(x)+ϱ2(x))u~2(t,x),t>0,xΩ,0=DvΔu~3(t,x)u~3(t,x)+Λ(t,x)α2β(t,x)pu~2(t,x)u~3(t,x)pu2(t,x)+l[N(x)u2(t,x)]μv(t,x)u~3(t,x),t>0,xΩ,0=DvΔu~5(t,x)μv(t,x)u~5(t,x)+(1τ2(t))ΩΓ(t,tτ2(t),x,y)α2β(tτ2(t),y)pu~2(tτ2(t),y)u~3(tτ2(t),y)pu~2(tτ2(t),y)+l[N(y)u~2(tτ2(t),y)]dy,t>0,xΩ,u~2(t,x)ν=u~3(t,x)ν=u~5(t,x)ν=0,t>0,xΩ.(17) Based on Section 4.1, there always exists E2=(0,u~2(t,x),u~3(t,x),0,u~5(t,x)) of (Equation6). Let (18) u1(t,x)=w1(t,x),u2(t,x)=u~2(t,x)+w2(t,x),u3(t,x)=u~3(t,x)+w3(t,x),u~4(t,x)=w4(t,x),u5(t,x)=u~5(t,x)+w5(t,x).(18) Linearizing (Equation6) at (0,u~2(t,x),u~3(t,x),0,u~5(t,x)), we gain the system which contains only w1 and w4 decoupled from other equations (19) {w1(t,x)t=DhΔw1(t,x)+c1β(t,x)l[N(x)u~2(t,x)]w4(t,x)pu~2(t,x)+l[N(x)u~2(t,x)](μh(x)+ϱ1(x))w1(t,x),xΩ,w4(t,x)t=DvΔw4(t,x)μv(t,x)w4(t,x)+(1τ1(t))ΩΓ(t,tτ1(t),x,y)α1β(tτ1(t),y)pu~3(tτ1(t),y)w1(tτ1(t),y)pu~2(tτ1(t),y)+l[N(y)u~2(tτ1(t),y)]dy,xΩ,w1(t,x)ν=w4(t,x)ν=0,xΩ,(19) for t>0. Define F^1(t):X2X2 by F^1(t)(ϕ1ϕ4)=(c1β(t,)l[N()u~2(t,)]ϕ4(0,)pu~2(t,)+l[N()u~2(t,)](1τ1(t))ΩΓ(t,tτ1(t),,y)α1β(tτ1(t),y)pu~3(tτ1(t),y)ϕ1(τ1(t),y)pu~2(tτ1(t),y)+l[N(y)u~2(tτ1(t),y)]dy),for tR, ϕ1,ϕ4E. Then we have [L^1v^](t):=0Ψ1(t,ts)F^1(ts)v^(ts+)ds,tR,v^Cω(R,X2),where the definition of Ψ1(t,ts) is consistent with that in Section 4.1 and v^(t) is the spatial distribution of infective humans and mosquitoes with strain 1 introduced at time tR. Obviously, L^1 on Cω(R,X2) is positive and bounded. From the Section 4.1, we know that (20) R^1:=r(L^1).(20) Let Q1t be the solution map of (Equation19) on X2 for t0, that is, Q1t(ϕ)=wt(x,ϕ), where wt(ϕ)(θ,x)=(w1(t+θ,x,ϕ),w4(t,x,ϕ)), θ[τ^,0] and w(t,ϕ) is the unique solution of (Equation19) with w(θ,x)=ϕ(θ,x) for all θ[τ^,0] and xΩ¯. Subsequently, Q1:=Q1ω is Poncaré map related to (Equation19), i.e. Q1(ϕ)=wω(ϕ). r(Q1) denotes the spectral radius of Q1.

Analogously, linearizing system (Equation6) at (u¯1(t,x),0,u¯3(t,x),u¯4(t,x),0), one has the following system containing only w2 and w5 decoupled from other equations (21) {w2(t,x)t=DhΔw2(t,x)+c2β(t,x)l[N(x)u¯1(t,x)]w5(t,x)pu¯1(t,x)+l[N(x)u¯1(t,x)](μh(x)+ϱ2(x))w2(t,x),t>0,xΩ,w5(t,x)t=DvΔw5(t,x)μv(t,x)w5(t,x)+(1τ2(t))ΩΓ(t,tτ2(t),x,y)α2pβ(tτ2(t),y)u¯3(tτ2(t),y)w2(tτ2(t),y)pu¯1(tτ2(t),y)+l[N(y)u¯1(tτ2(t),y)]dy,t>0,xΩ,w2(t,x)ν=w5(t,x)ν=0,t>0,xΩ.(21) For t0, denote Q2t to be the solution map of (Equation21) on X2, that is, Q2t(ϕ)=wt(ϕ,x), where wt(ϕ)(θ,x)=(w2(t+θ,x,ϕ),w5(t,x,ϕ)), θ[τ^,0], xΩ¯ and w(t,x,ϕ) is the unique solution of system (Equation21) with w(θ,x)=ϕ(θ,x) for all θ[τ^,0]. Let Q2:=Q2ω be the Poncaré map associated with system (Equation22), i.e. Q2(ϕ)=wω(ϕ,x), and r(Q2) denote the spectral radius of Q2. Use Q¯1t to be the solution map of (Equation19) on space E1, then Q¯1:=Q¯1ω be Poincaré map related to (Equation19), and r(Q¯1) be the spectral radius of Q¯1. Set Q¯2t to be the solution map of (Equation22), then Q¯2:=Q¯2ω is Poincaré map related to (Equation21) on space E2, and r(Q¯2) is the spectral radius of Q¯2. We know that there exists a positive linear operator L^2 on Cω(R,X2), defined by [L^2vˇ](t):=0Ψ2(t,ts)F^2(ts)vˇ(ts+)ds,tR,vˇCω(R,X2),where Ψ2(t,s)=diag(T2(t,s),T3(t,s)), vˇ(t) is the spatial distribution of infected humans and mosquitoes with strain 2 introduced at time tR and F^2(t)(ϕ2ϕ5)=(c2β(t,)l[N()u¯1(t,)]ϕ5(0,)pu¯1(t,)+l[N()u¯1(t,)](1τ2(t))ΩΓ(t,tτ2(t),,y)α2β(tτ2(t),y)pu¯3(tτ2(t),y)ϕ2(τ2(t),y)pu¯1(tτ2(t),y)+l[N(y)u¯1(tτ2(t),y)]dy).Then R^2:=r(L^2).Similar to [Citation43, Lemma 3.4] and [Citation19, Theorem 3.7], one has the following result.

Lemma 4.8

R^i and r(Qi)1 (i=1,2) have the same sign.

According to [Citation22, Lemma 3.8], [Citation2, Lemma 5] and [Citation42, Proposition 1.1], the following results can be obtained.

Lemma 4.9

Two Poincaré maps Qi:X2X2 and Q¯i:=EiEi have the same spectral radius, i.e. r(Qi)=r(Q¯i) (i=1,2).

We wish to explore the long-time behaviour of (Equation6), based on the evidence that the system has a special solution.

Lemma 4.10

There is a v~(t,x), which is positive ω-periodic function, in a manner that eμ2tv~(t,x) is a solution of (Equation22), where μ2=lnr(Q¯2)ω.

According to [Citation46, Proposition 3.11], we have the following lemma.

Lemma 4.11

If R^i>1, then Ri>1 for i = 1, 2.

5. Threshold dynamics

The main focus of this position is to explore the threshold dynamics of (5) based on Ri and R^i (i=1,2). The stability analysis and persistent results contribute to a discernment of the interactions between the disease dynamics as well as possible coexistence of the two strains.

Theorem 5.1

If R1<1 and R2<1, then (0,0,u3(t,x),0,0) of (Equation6) is globally attractive.

The conclusion of this theorem can be obtained directly from the conclusions of the following two theorems.

Let X4=C([τ1(0),0],Yh)×C([τ1(0),0],Y+)×Y+, and X5=C([τ2(0),0],Yh)×C([τ2(0),0],Y+)×Y+. In the light of Lemma 4.4, for any ϕ¯X3 and ψ^X4 with ϕ¯1(θ,)=ψ^1(θ,), θ[τ1(0),0], ϕ¯2(θ,)=ψ^2(θ,), θ[τ1(0),0], and ϕ¯3()=ψ^3(0,), u¯(t,,ϕ¯)=z¯(t,,ψ^) holds, t0, where u¯(t,,ϕ¯) and z¯(t,,ψ^) are solutions of system (Equation10) meeting u¯0=ϕ¯ and z¯0=ψ^, respectively. Consequently, the solution of system (Equation10) on X4 exists globally on [0,+) and ultimately bounded.

Theorem 5.2

  1. As R1<1, E01 is globally attractive for (Equation10);

  2. As R1>1, system (Equation10) has at least one positive ω-periodic solution E1, and there is a constant γ1>0, making φ¯X4 with φ¯1(0,)0 and φ¯3()0, one has lim inftminxΩ¯u¯i(t,x,ϕ˘)γ1,i=1,2,3.

Proof of Theorem 5.2 is given in Appendix 3. Similarly, the following conclusion holds.

Theorem 5.3

  1. When R2<1, E02 is globally attractive for system (Equation11);

  2. When R2>1, system (Equation11) has at least one positive ω-periodic solution E2, and there is a constant γ5>0 such that, for any φˇX5 with φˇ1(0,)0 and φˇ3()0, one has lim inftminxΩ¯u^i(t,x,φˇ)γ5,i=1,2,3.

Furthermore, by [Citation22, Lemma 3.5], [Citation15] and [Citation24, Theorem 2.9], the following result is given.

Lemma 5.1

Denote Dt to be the solution map of system (Equation10) on X4, that is, for each ψ^X4, Dt(ψ^)=u¯t(ψ^), t0. Let Dt be an ω-periodic semi-flow on X4, in this sense, one has: (i) D0=I; (ii) t0, Dt+ω=DtDω holds; and (iii) Dt(ψ^) is continuous in (t,ψ^)[0,)×X4. Additionally, in X4, D:=Dω has a strong global attractor.

Lemma 5.2

Use (u¯1(t,,ψ^),u¯2(t,,ψ^),u¯3(t,,ψ^)) to be the solution of (Equation10) with initial data ψ^X4, if it exists t10, such that u¯i(t1,x,ψ^)0, (i=1,3), in that way, the solution of (Equation10) meets u¯i(t,x,ψ^)>0,t>t1,xΩ¯.Besides, for nontrivial data ψ^X4, one has u¯2(t,x,ψ^)>0, t>0, xΩ¯, and lim inftu¯2(t,x,ψ^)ε1,uniformlyforxΩ¯,where ε1 is a ψ^-independent positive constant.

The proof of Lemma 5.2 is given in Appendix 4.

5.1. Competitive exclusion and coexistence

Denote τ¯=min{τ1(0),τ2(0)}, X6=C([τ¯,0],Yh)×C([τ¯,0],Yh)×C([τ¯,0],Y+)×Y+×Y+.Next, we prove the competitive exclusion and coexistence of (Equation6).

For each ϕX1+ and ψˇX6 with ϕ1(η1,)=ψˇ1(η1,), ϕ2(η1,)=ψˇ2(η1,), ϕ3(η1,)=ψˇ3(η1,), ϕ4(0,)=ψˇ4(0,), ϕ5(0,)=ψˇ5(0,), for η1[τ¯,0]. We have u(t,,ϕ)=uˇ(t,,ψˇ), for t0, where u(t,,ϕ) and uˇ(t,,ψˇ) are solutions of system (Equation6) satisfying u0=ϕ and uˇ0=ψˇ, respectively. Consequently, the solution of system (Equation6) on X6 exists globally on [0,+) and ultimately bounded. Furthermore, in view of those in [Citation22, Lemma 3.5] and [Citation15] together with [Citation24, Theorem 2.9], we obtain the following result.

Lemma 5.3

Denote Φt to be the solution map of system (Equation6) on X6, that is, for ψˇX6, Φt(ψˇ)=ut(ψˇ), t0. Then Φt is an ω-periodic semi-flow on X6 in this sense, one has: (i) Φ0=I; (ii) Φt+ω=ΦtΦω, t0; and (iii) Φt(ψˇ) is continuous in (t,ψˇ)[0,)×X6. Additionally, Φ:=Φω has a strong global attractor in X6.

Based on the comparison principle and Lemma 5.3, we see that the solution of system (Equation6) is strictly positive.

Lemma 5.4

Use (u1(t,,ψˇ),u2(t,,ψˇ),u3(t,,ψˇ),u4(t,,ψˇ),u5(t,,ψˇ)) to be the solution of (Equation6) with the initial data ψˇX6, if there is some t30 such that uj(t3,,ψˇ)0, in that way, the solution of (Equation6) meets uj(t,x,ψˇ)>0,foranyt>t3,j=1,2,4,5,xΩ¯.Besides, for any nontrivial initial data ψˇX6, one has u3(t,x,ψˇ)>0, for t>0, xΩ¯, and lim inftu3(t,x,ψˇ)ε2,uniformlyforxΩ¯,where ε2>0 is a ψˇ-independent constant.

Theorem 5.4

Under the condition of R2>1>R1, suppose (u1(t,x,ψˇ),u2(t,x,ψˇ),u3(t,x,ψˇ),u4(t,x,ψˇ),u5(t,x,ψˇ)) is the solution of (Equation6) with the initial data ψˇ=(ψˇ1,ψˇ2,ψˇ3,ψˇ4,ψˇ5)X6. Then there is a P1>0, such that, for each ψˇX6 with ψˇ2(0,)0 and ψˇ5()0, one has limtu1(t,x,ψˇ)=0, limtu4(t,x,ψˇ)=0 and lim inftminxΩ¯u2(t,x,ψˇ)P1,lim inftminxΩ¯u5(t,x,ψˇ)P1.

Proof.

As that in Theorem 5.3(a), based on R1<1, we can prove limtu1(t,x,ψˇ)=0 and limtu4(t,x,ψˇ)=0. According to the proof of Theorem 5.3(b), one has lim inftminxΩ¯u2(t,x,ψˇ)P1 and lim inftminxΩ¯u5(t,x,ψˇ)P1 by virtue of R2>1.

Theorem 5.5

Suppose that R1>1>R2, and (u1(t,x,ψˇ),u2(t,x,ψˇ),u3(t,x,ψˇ),u4(t,x,ψˇ),u5(t,x,ψˇ)) is the solution of (Equation6) with the initial data ψˇ=(ψˇ1,ψˇ2,ψˇ3,ψˇ4,ψˇ5)X6. Then there is a P2>0 such that for each ψˇX6 with ψˇ1(0,)0 and ψˇ4()0, one has limtu2(t,x,ψˇ)=0, limtu5(t,x,ψˇ)=0 and lim inftminxΩ¯u1(t,x,ψˇ)P2,lim inftminxΩ¯u4(t,x,ψˇ)P2.

The persistence of (Equation6) is shown. Given ψˇX6, set u(t,x,ψˇ) to be the unique solution of system (Equation6) with u0=ψˇ. The following result shows the results of (Equation6) based on R^i(i=1,2).

Theorem 5.6

Presume that R1>1 and R2>1, so that the boundary periodic solutions of system (Equation6) exist.

  1. If R^1>1, then (0,u~2(t,x),u~3(t,x),0,u~5(t,x)) is unstable. There exists a P3, then for any ψˇX6 with ψˇ1(0,)0 and ψˇ4()0, one has limtu2(t,x,ψˇ)=0, limtu5(t,x,ψˇ)=0 and lim inftminxΩ¯u1(t,x,ψˇ)P3,lim inftminxΩ¯u4(t,x,ψˇ)P3.

  2. If R^2>1, then the ω-periodic solution of strain 1, (u¯1(t,x),0,u¯3(t,x),u¯4(t,x),0), is unstable. There is a P4 such that for any ψˇX6 with ψˇ2(0,)0 and ψˇ5()0, we have limtu1(t,x,ψˇ)=0, limtu4(t,x,ψˇ)=0 and lim inftminxΩ¯u2(t,x,ψˇ)P4,lim inftminxΩ¯u5(t,x,ψˇ)P4.

Theorem 5.7

As R1>1, R2>1, R^1>1 and R^2>1, system (Equation6) admits at least one positive solution (u¯1(t,x),u¯2(t,x), u¯3(t,x),u¯4(t,x),u¯5(t,x)), which is ω-periodic, and there is a constant η~>0 such that, for any ψˇ=(ψˇ1,ψˇ2,ψˇ3,ψˇ4,ψˇ5)X6, meets ψˇ1(0,)0, ψˇ2(0,)0, ψˇ4()0, ψˇ5()0, we have lim inftminxΩ¯ui(t,x,ψ)η~,i=1,2,3,4,5.

Proof.

By virtue of Lemma 4.11 and R^i>1 (i=1,2), one has Ri>1 (i=1,2). Let Z0:={ψˇX6:ψˇ1(0,)0andψˇ2(0,)0andψˇ4()0andψˇ5()0},and Z0:=X6Z0={ψˇX6:ψˇ1(0,)0orψˇ2(0,)0orψˇ4()0orψˇ5()0}.Notice that for any ψˇZ0, Lemma 5.4 indicates that ui(t,x,ψˇ)>0 (i=1,2,4,5) t>0, xΩ¯. Thus, we get Φn(Z0)Z0, nN. Lemma 5.3 reveals that Φ has a strong global attractor in X6.

Define Z:={ψˇZ0:Φn(ψˇ)Z0,nN},and ω(ψˇ) to be the omega limit set and the corresponding orbit is γ+(ψˇ):={Φn(ψˇ):nN}. Denote M1={(0^,0^,u3(t,),0,0)}, M2={(u¯1(t,),0^,u¯3(t,),u¯4(t,),0)} and M3={(0^,u~2(t,),u~3(t,),0,u~5(t,))}, t>0, where 0^ represents the constant function identically zero, 0^(η1,)=0, η1[τ¯,0]. Next, we show the following claims.

Claim 1. ϑ¯Zω(ϑ¯)=M1M2M3, ϑ¯Z, where ω(ϑ¯) is the omega limit set and the corresponding orbit is γ+={Φn(ϑ¯):nN} of system (Equation6) for ϑ¯Z.

According to the definition of Z, we obtain Φn(ϑ¯)Z0 for ϑ¯Z. Then, either u1(nω,,ϑ¯)0 or u2(nω,,ϑ¯)0 or u4(nω,,ϑ¯)0 or u5(nω,,ϑ¯)0, for nN. What's more, by contradiction and Lemma 5.4, it is evident that for t0, either u1(t,,ϑ¯)0 or u2(t,,ϑ¯)0 or u4(t,,ϑ¯)0 or u5(t,,ϑ¯)0. In the case where u1(t,,ϑ¯)0 for t0, we get that the forth equation in (Equation6) satisfies u4(t,x,ϑ¯)tDvΔu4(t,x,ψ)μ¯vu4(t,x,ϑ¯),where μ¯v is defined in Section 3. Based on the comparison principle, one has limtu4(t,x,ϑ¯)=0 uniformly for xΩ¯. In this instance, u2(t,x,ϑ¯) and u5(t,x,ϑ¯) have the following possibilities:

  1. If u2(t,,ϑ¯)0, then limtu5(t,x,ϑ¯)=0 uniformly for xΩ¯. Consequently, u3 equation abies by an autonomous system, which is asymptotic to the system (Equation12). Again, according to the theory of internal chain transitive sets [Citation45], one can obtain that limt(u3(t,x,ϑ¯)u3(t,x,ϑ¯))=0 uniformly for xΩ¯.

  2. If there is some t40, such that u2(t4,,ϑ¯)0, then u2(t,,ϑ¯)>0, tt4, based on Lemma 5.4. Thus, we get u5(t,,ϑ¯)0 or u5(t5,,ϑ¯)0 for some t50. Let's discuss the following situations.

    1. If u5(t,,ϑ¯)0, then we obtain limtu2(t,x,ϑ¯)=0 uniformly for xΩ¯, which contracts u2(t4,,ϑ¯)0, tt4.

    2. If there is some t50, such that u5(t5,,ϑ¯)0, then based on Lemma 5.4, one has u5(t,,ϑ¯)>0, for tt5. Then, equations u2, u3 and u5 satisfy a nonautonomous system, which is asymptotic to system (Equation17). Additionally, according to the theory of internal chain transitive sets [Citation45], one has limt((u2(t,x,ϑ¯),u3(t,x,ϑ¯),u5(t,x,ϑ¯))(u~2(t,x),u~3(t,x),u~5(t,x)))=0. uniformly for xΩ¯.

Under the above discussion, we get ω(ϑ¯)=M1M3.

If there is some t60, such that u1(t6,,ϑ¯)0, then according to Lemma 5.4, we obtain u1(t,,ϑ¯)>0, tt6. Consequently, one has u2(t,,ϑ¯)0 or u4(t,,ϑ¯)0 or u5(t,,ϑ¯)0.

If u4(t,,ϑ¯)0, from the fifth equation of system (Equation6), we can get limtu1(t,x,ϑ¯)=0, which is contradict with the fact that u1(t,,ϑ¯)>0, tt6. Thus, there is t70 such that u4(t7,,ϑ¯)0, then we get u4(t,,ϑ¯)>0 from Lemma 5.4 for tt7. Therefore, u2(t,,ϑ¯)0 or u5(t,,ϑ¯)0 must be established. Here we only discuss the case that u2(t,,ϑ¯)0 holds, and the same is true for u5(t,,ϑ¯)0. If u2(t,,ϑ¯)0, for t0, according to comparison principle, one has limtu5(t,x,ϑ¯)=0 uniformly for xΩ¯. So, u3 satisfies a nonautonomous system, which is asymptotic to the second equation of periodic system (Equation22).

Consequently, we get ω(ϑ¯)=M2. Thus, Claim 1 holds.

Claim 2. M1 is a uniformly weak repeller for Z0, in this sense, on has lim suptΦn(ψˇ)M1σ3,ψˇZ0,for σ3 small enough.

Claim 3. Mi is uniformly weak repeller for Z0 (i=2,3), in this sense, there is a sufficiently small σ4>0 meeting lim suptΦn(ψˇ)Miσ4,ψˇZ0,We only give the proof of M2, the same proof step is also applicable to M3. For ϵ~>0, consider the following equation with parameter ϵ~ (22) {v1ϵ~(t,x)t=DhΔv1ϵ~(t,x)+c2β(t,x)l[N(x)u¯1(t,x)2ϵ~]v2ϵ~(t,x)p[u¯1(t,x)+2ϵ~]+l[N(x)u¯1(t,x)+ϵ~](μh(x)+ϱ2(x))v1ϵ~(t,x),xΩ,v2ϵ~(t,x)t=DvΔv2ϵ~(t,x)μv(t,x)v2ϵ~(t,x)+(1τ2(t))ΩΓ(t,tτ2(t),x,y)α2β(tτ2(t),y)p[u¯3(tτ2(t),y)ϵ~]v1ϵ~(tτ2(t),y)p[u¯1(tτ2(t),y)+2ϵ~]+l[N(y)u¯1(tτ2(t),y)+ϵ~]dy,xΩ,v1ϵ~(t,x)ν=v2ϵ~(t,x)ν=0,xΩ.(22) Set vϵ~(t,x,φ)=(v1ϵ~(t,x,φ),v2ϵ~(t,x,φ)) to be the unique solution of (Equation22), for φE2, with v0ϵ~(φ)(θ2,x)=(φ1(θ2,x),φ2(0,x)) for all θ2[τ2(0),0] xΩ¯, where vtϵ~(φ)(θ2,x)=vϵ~(t+θ2,x,φ)=(v1ϵ~(t+θ2,x,φ),v2ϵ~(t,x,φ)),t0,(θ2,x)[τ2(0),0]×Ω¯.Let Q¯2ϵ~:E2E2 be the Poincaré map of (Equation22), i.e. Q¯2ϵ~(φ)=vωϵ~(φ), φE2, and r(Q¯2ϵ~) be the spectral radius of Q¯2ϵ~. Since limϵ~0r(Q¯2ϵ~)=r(Q¯2)>1, fixing a sufficiently small constant ϵ~>0, such that ϵ~<min{mint[0,ω],xΩ¯u¯3(t,x),mint[0,ω],xΩ¯N(x)u¯1(t,x)2}andr(Q¯2ϵ~)>1.Based on the continuous dependence of solutions on the initial data, for fixed ϵ~>0, there is ϵ~>0, then for ψˇ with ψˇM2<ϵ~, and t[0,ω], we arrive that Φt(ψˇ)Φt(M2)<ϵ~. Now, by contradiction, we state that M2 is uniformly weak repeller for Z0.

Assuming contradictions, for some ψˇ0Z0, one has lim supnΦn(ψˇ0)M2<ϵ~. There is n21, then Φn(ψˇ0)M2<ϵ~ for nn2. For each tn2ω, let t=nω+t with n=[t/ω] and t[0,ω), one has (23) Φt(ψˇ0)Φt(M2)=Φt(Φn(ψˇ0))Φt(M2)<ϵ~.(23) Following (Equation23) and Lemma 5.4, (24) u¯1(t,x)ϵ~<u1(t,x)<u¯1(t,x)+ϵ~,0<u2(t,x)<ϵ~,u¯3(t,x)ϵ~<u3(t,x)<u¯3(t,x)+ϵ~,u¯4(t,x)ϵ~<u4(t,x)<u¯4(t,x)+ϵ~,0u5(t,x)<ϵ~,(24) for each tn2ωτ^ and xΩ¯. As a consequence, as tn2ω, u2(t,x,ψˇ0) and u5(t,x,ψˇ0) meet (25) {u2(t,x)tDhΔu2(t,x)+c2β(t,x)l[N(x)u¯1(t,x)2ϵ~]u5(t,x)p[u¯1(t,x)+2ϵ~]+l[N(x)u¯1(t,x)+ϵ~](μh(x)+ϱ2(x))u2(t,x),xΩ,u5(t,x)tDvΔu5(t,x)μv(t,x)u5(t,x)+(1τ2(t))ΩΓ(t,tτ2(t),x,y)α2β(tτ2(t),y)p[u¯3(tτ2(t),y)ϵ~]u2(tτ2(t),y)p[u¯1(tτ2(t),y)+2ϵ~]+l[N(y)u¯1(tτ2(t),y)+ϵ~]dy,xΩ,u2(t,x)ν=u5(t,x)ν=0,xΩ.(25) Based on u(t,x,ψˇ0)0 for xΩ¯ and t0, there must be a constant σ5>0 then (26) (u2(t,x,ψˇ0),u5(t,x,ψˇ0))σ5eμ2ϵ~tv~ϵ~(t,x),t[n2ωτ^,n2ω],xΩ¯,(26) where v~ϵ~(t,x) is a positive ω-periodic function then eμ2ϵ~tv~ϵ~(t,x) is a solution of system (Equation23), where μ2ϵ~=lnr(Q¯2ϵ~)ω. It follows the comparison theorem, that (u2(t,x,ψˇ0),u5(t,x,ψˇ0))σ5eμ2ϵ~tv~ϵ~(t,x),tn2ω,xΩ¯.Since μ2ϵ~>0, then u2(t,,ψˇ0) and u5(t,,ψˇ0) as t can be obtained directly, which leads to a contradiction. As a consequence, Claim 3 holds.

With the above discussion, it is easy to obtain that ϝ:=M1M2M3 is isolated and invariant for Φ in X6, and Ws(ϝ)Z0=, where Ws(ϝ) is the stable set of Γ for Φ. It then follows from [Citation45, Theorem 1.3.1 and Remark 1.3.1] that Φ is uniformly persistent with respect to (Z0,Z0), in this sense, there is an η^>0, (27) lim inftd(Φ(ψˇ),Z0)η^,ψˇZ0.(27) Since Φn:=Φnω is compact, for n with nω>τ^, then, Φ is asymptotically smooth. Lemma 5.3 also indicates that Φ has a global attractor on X6. According to [Citation24, Theorem 3.7], Φ admits a global attractor A^0 in Z0.

Now we prove practical persistence required. Following A^0=Φ(A^0), we obtain that ψˇ1(0,)>0, ψˇ2(0,)>0, ψˇ4()>0 and ψˇ5()>0 for any ψˇA^0. Set B^:=t[0,ω]Φt(A^0). Obviously, we have B^Z0 and limtd(Φt(ψˇ),B^)=0, ψˇZ0. A continuous function q2:X6R+ is defined as q2(ψ):=min{minxΩ¯ψˇ1(0,x),minxΩ¯ψˇ2(0,x),minxΩ¯ψˇ4(x),minxΩ¯ψˇ5(x)},ψˇ=(ψˇ1,ψˇ2,ψˇ3,ψˇ4,ψˇ5)X6.Owing to B^ is compact subset of Z0, then infψˇB^q2(ψˇ)=minψˇB^q2(ψˇ)>0. Thence, there is η>0 such that lim inftq2(Φt(ψˇ))=lim inftmin(minxΩ¯u1(t,x,ψˇ),minxΩ¯u2(t,x,ψˇ),minxΩ¯u4(t,x,ψˇ),minxΩ¯u5(t,x,ψˇ))η,ψˇZ0.Furthermore, according to Lemma 5.4, there must be an η~(0,η), then lim inftminxΩ¯ui(t,x,ψˇ)η~,i=1,2,3,4,5.The existence of a positive periodic steady state remains to be proved. Based on [Citation2, Lemma 8] and [Citation45, Theorem 3.5.1], one has that for t>0, solution map Φ~t:WhWh of (Equation6), defined in Lemma 3.3, is a κ-contraction. Define Q0:={ψˇWh:ψˇ1(0,)0andψˇ2(0,)0andψˇ4()0andψˇ5()0},and Q0:=Wh/Q0={ψˇWh:ψˇ1(0,)=0orψˇ2(0,)=0orψˇ4()=0orψˇ5()=0}.Then Φ~ is ρ~2-uniformly persistent with ρ~2(ψˇ)=d(ψˇ,Q0), ψˇQ0, is easily attainable. According to [Citation24, Theorem 4.5], system (Equation6) has an ω-periodic solution (z1#(t,),z2#(t,),z3#(t,),z4#(t,),z5#(t,)) with (z1t#,z2t#,z3t#,z4t#,z5t#)Q0. Set u¯i(η1,)=zi#(η1,), (i=1,2,3), and u¯j(0,)=zj#(0,), (j=4,5) with θ[τ¯,0]. Again, due to the uniqueness of solutions, we get that (u1(t,,ψˇ),u2(t,,ψˇ),u3(t,,ψˇ),u4(t,,ψˇ),u5#(t,,ψˇ)) is a positive periodic solution of system (Equation6).

Next, we prove the asymptomatic behaviour of Evi(t,x) (i=1,2) in system (5). When R1<1 and R2<1, we obtain limt(I1(t,x),I2(t,x),Sv(t,x),Iv1(t,x),Iv2(t,x)(0,0,u3(t,x),0,0))=0,uniformly for xΩ¯. By (Equation8), limtEv1(t,x)=limtEv2(t,x)=0,uniformlyforxΩ¯.Presume that R1>1 and R2>1.

(a) If R^1>1, for each ψˇX6 with ψˇ1(0,)0 and ψˇ4()0, we have limtI2(t,x,ψˇ)=0, limtIv2(t,x,ψˇ)=0 and lim inftminxΩ¯I1(t,x,ψˇ)P3,lim inftminxΩ¯Iv1(t,x,ψˇ)P3.Based on the integral forms of (Equation8), there is P5>0, lim inftminxΩ¯Ev1(t,x)P5,lim inftminxΩ¯Ev2(t,x)=0.(b) If R^2>1, for each ψˇX6 with ψˇ2(0,)0 and ψˇ5()0, we have limtI1(t,x,ψˇ)=0, limtIv1(t,x,ψˇ)=0 and lim inftminxΩ¯I2(t,x,ψˇ)P4,lim inftminxΩ¯Iv2(t,x,ψˇ)P4.Based on the integral forms of (Equation8), there is P6>0, such that lim inftminxΩ¯Ev1(t,x)=0,lim inftminxΩ¯Ev2(t,x)P6.(c) As R^1>1 and R^2>1, for each ψˇX6 with ψˇi(0,)0, i = 1, 2, and ψˇj()0 j = 4, 5, we get lim inftIi(t,x,ψˇ)η~,lim inftSv(t,x,ψˇ)η~,lim inftIvi(t,x,ψˇ)η~,i=1,2,uniformly for xΩ¯. On the basis of the integral forms of (Equation8), there exists a constant η2>0, such that lim inftminxΩ¯Ev1(t,x)η2,lim inftminxΩ¯Ev2(t,x)η2,with Ev1(0,x) and Ev2(0,x) meeting compatibility conditions (Equation7). Besides, if (I1(t,x),I2(t,x),Sv(t,x),Iv1(t,x),Iv2(t,x)) is ω-periodic in t, then Ev1(t,x) and Ev2(t,x) are also ω-periodic in t.

6. Numerical simulations

In this position, numerical simulations exhibit how some epidemiological insights can be gained from our analytical results. Without loss of generality, we make use of one-dimensional domain Ω=[0,π] to simulate the long-time behaviour which is inspired by Lou and Zhao [Citation21], Bai et al. [Citation2] and Wu et al. [Citation40], and apply system (Equation6) to the spread of malaria in Maputo Province, Mozambique. Set the periodic ω=12 months.

6.1. Numerical verification of theoretical analysis

First, we illustrate the theoretical results obtained in the previous sections. The reproduction numbers usually (but not always) control the outcome of strain competition. The basic reproduction numbers are calculated based on the theory provided by [Citation19] and the expected results are summarized in Table .

Table 1. The potential dynamical outcomes of system (Equation6).

In the next subsections, the theoretical results are verified one by one using the parameter values defined in Appendix 5.

Figure 1. The evolution of infection compartments of humans and mosquitoes when R1<1 and R2<1. (a) The evolution of u1. (b) Then evolution of u4. (c) Then evolution of u2. (d) Then evolution of u5.

Figure 1. The evolution of infection compartments of humans and mosquitoes when R1<1 and R2<1. (a) The evolution of u1. (b) Then evolution of u4. (c) Then evolution of u2. (d) Then evolution of u5.

Figure 2. The evolution of infection compartments of humans and mosquitoes when R1<1 and R2>1. (a) The evolution of u1. (b) The evolution of u4. (c) The evolution of u2. (d) The evolution of u5.

Figure 2. The evolution of infection compartments of humans and mosquitoes when R1<1 and R2>1. (a) The evolution of u1. (b) The evolution of u4. (c) The evolution of u2. (d) The evolution of u5.

Figure 3. The evolution of infection compartments of humans and mosquitoes when R1>1 and R2>1. (a) The evolution of u1. (b) Then evolution of u4. (c) Then evolution of u2. (d) Then evolution of u5.

Figure 3. The evolution of infection compartments of humans and mosquitoes when R1>1 and R2>1. (a) The evolution of u1. (b) Then evolution of u4. (c) Then evolution of u2. (d) Then evolution of u5.

Figure 4. The evolution of infection compartments of humans and mosquitoes when R1>1 and R2>1. (a) The evolution of u1. (b) The evolution of u4. (c) The evolution of u2. (d) The evolution of u5.

Figure 4. The evolution of infection compartments of humans and mosquitoes when R1>1 and R2>1. (a) The evolution of u1. (b) The evolution of u4. (c) The evolution of u2. (d) The evolution of u5.

Figure 5. The evolution of infection compartments of humans and mosquitoes when R1>1 and R2>1. (a) The evolution of u1. (b) The evolution of u4. (c) The evolution of u2. (d) The evolution of u5.

Figure 5. The evolution of infection compartments of humans and mosquitoes when R1>1 and R2>1. (a) The evolution of u1. (b) The evolution of u4. (c) The evolution of u2. (d) The evolution of u5.

6.1.1. Case 1 in Table 1

We choose c1=0.0151, c2=0.01454, α1=0.0146, α2=0.012. Let ϱ1=a1(1.05cos(2x))Month1 and ϱ2=a2(1.05cos(2x))Month1, where a1=0.085, a2=0.081, reflect the fact that people living in urban areas (around the centre) can enjoy added medical treatment (the number of physicians and hospitals, supply of medicines, state-of-the-art medical equipment) than people in rural area, resulting in a higher recovery rate around the centre, and maintaining the other parameters listed in Table . To compute R1 and R2, the numerical scheme recently proposed in [Citation19, Lemma 2.5 and Remark 3.2] is used. For this set of parameters, we can numerically calculate R1=0.3644<1, R2=0.3733<1. Figure  shows numerical plots of infected compartments following the initial data u(θ,x)=(11346540000cos2x134655386cos2x6807900113465cos2x890000100000cos2x177206158cos2x),θ[τ^,0],x[0,π].Then Figure  implies that two strains die out, which is coincident with Theorem 5.1. Infectious humans and mosquitoes go to 0, indicating that the disease will be eliminated under this set of parameters.

6.1.2. Case 2 and Case 3 in Table 1

Here, we simulate the result of Case 2 in Table , and then we can handle Case 3 similarly by adjusting the parameters. We choose ϱ1=a1(1.05cos(2x))Month1, ϱ2=a2(1.05cos(2x))Month1, where a1=0.06, a2=0.05, c1=0.02, c2=0.025, α1=0.03, α2=0.04, R1=0.7045<1, R2=1.1136>1. Figure  shows that if R1<1 and R2>1, strain 1 becomes extinct and strain 2 persists, which is consistent with Theorem 5.4. Note that we truncate the time interval by [30, 60]. Adjust the values of parameters befittingly to satisfy R1>1 and R2<1, and the conclusion of Case 3 in Table  is obtained. Details are not mentioned here.

6.1.3. Case 4 in Table 1

For Case 4 in Table , there are three results for R1>1 and R2>1. For the first conclusion, we choose ϱ1=a1(1.05cos(2x))Month1, ϱ2=a2(1.05cos(2x))Month1, where a1=0.07, a2=0.05, c1=0.09, c2=0.1, α1=0.04, α2=0.05, and maintain other parameters recorded in Table . Based on this set of parameters, numerically compute R1=1.6079>1 and R2=2.2901>1. Figure  indicates that strain 1 becomes extinct and strain 2 is persistent. Note that we truncate the time interval by [30, 60]. Because under this set of parameters, we assume that individuals and mosquitoes infected with strain 2 are more infectious and have a smaller recovery rate. Therefore, once this phenomenon occurs, there will be greater challenges to the control and eradication of malaria.

By adjusting the parameters, we can also get the result that strain 1 is persistent, and strain 2 is extinct under the condition of R1>1 and R2>1. We choose ϱ1=a1(1.05cos(2x))Month1, ϱ2=a2(1.05cos(2x))Month1, where a1=0.07, a2=0.06, c1=0.15, c2=0.1, α1=0.055, α2=0.05, and maintain other parameters recorded in Table . Then, compute R1=2.4342>1 and R2=2.2901>1. Figure  exhibits that strain 1 is persistent and strain 2 is extinct. Note that we truncate the time interval by [30, 60].

Next, the coexistence of two strains under the conditions of R1>1 and R2>1 is simulated. We choose ϱ1=a1(1.05cos(2x))Month1, ϱ2=a2(1.05cos(2x))Month1, where a1=0.025, a2=0.0505, c1=0.025, c2=0.04836, α1=0.035, α2=0.04, and keep other parameters listed in Table . For this set of parameters, we can compute R1=1.2759>1 and R2=1.5417>1. Figure  shows that two strains can coexist, which is consistent with Theorem 5.7. Note that we truncate the time interval [50,80] to demonstrate the existence of periodic solution, in which two strains coexist.

6.2. The impact of diffusion

Except for Dh=0 and Dv=0, the other parameters are in accordance with in Case 1. When Dh=0.1 and Dv=0.0125 in Case 1, then R1=0.3644 and R2=0.3733, and the disease dies out. Whereas, if Dh=0 and Dv=0, then R1=1.3858, R2=1.4100 are exceeding 1, and the disease persists. As a consequence, the result is shown in Figure .

Figure 6. The evolution of infection compartments of humans and mosquitoes when Dh=0 and Dv=0. (a) The evolution of u1. (b) The evolution of u4. (c) The evolution of u2. (d) The evolution of u5.

Figure 6. The evolution of infection compartments of humans and mosquitoes when Dh=0 and Dv=0. (a) The evolution of u1. (b) The evolution of u4. (c) The evolution of u2. (d) The evolution of u5.

Figure  shows some intriguing phenomena. In this case, malaria becomes extinct in urban areas, but persists in rural areas. It is caused by spatial heterogeneity. Because the urban environment is not apposite to mosquitoes to survive, there are relatively few mosquitoes, and there are more medical resources, the disease control in urban is relatively better. Compare with Figure , the movement will reduce disease spread risk, to a limited degree, which is consistent with the conclusion in [Citation40]. The implications include: (i) The presence of spread suggests that infected humans in rural areas can enjoy better treatment in cities. (ii) As stated in Ref. [Citation40], it is challenging for mosquitoes to collect blood meal in moving humans.

In order to understand the effect of human and mosquito movement on malaria transmission, under a set of the parameter values in Case 1, we do the following work:

(i) We assume that the diffusion coefficient of humans is Dh=0, and the diffusion coefficient of mosquitoes is Dv=0.0125 (see Figure ). Trough calculation, R1=1.1369 and R2=1.1719. Compared with the above Dv=0 and Dh=0 (Figure ), the diffusion of mosquitoes reduces the risk of disease transmission to some extent. Biologically, this may be due to the fact that there are more mosquitoes in rural areas, and infected mosquitoes move to urban areas. With the same bite rate, city has better medical conditions and a greater recovery rate, so the risk of malaria transmission will be decreased.

Figure 7. The evolution of infection compartments of humans and mosquitoes with Dh=0 and Dv=0.0125. (a) The evolution of u1. (b) The evolution of u4. (c) The evolution of u2. (d) The evolution of u5.

Figure 7. The evolution of infection compartments of humans and mosquitoes with Dh=0 and Dv=0.0125. (a) The evolution of u1. (b) The evolution of u4. (c) The evolution of u2. (d) The evolution of u5.

(ii) We assume that the diffusion coefficient of humans is Dh=0.1, but the diffusion coefficient of mosquitoes is Dv=0 (see Figure ). Then R1=0.3645 and R2=0.3734. Compared with the case of Dh=0 and Dv=0 (Figure ), humans diffusion reduces the risk of disease transmission. Biologically, it can be explained that the diffusion of humans gives them the opportunity to find a better medical environment, which will increase the recovery rate and reduce the risk of disease transmission. This may be due to the better medical conditions in the city, and the infected people will go to the city for treatment, which will increase the recovery rate and reduce the risk of disease transmission. Compared with Dh=0 and Dv=0.0125, humans diffusion has a greater impact on disease transmission. In summary, the risk of malaria transmission will be overestimated if the diffusion is ignored in the model analysis.

Figure 8. The evolution of infection compartments of humans and mosquitoes with Dh=0.1 and Dv=0. (a) The evolution of u1. (b) The evolution of u4. (c) The evolution of u2. (d) The evolution of u5.

Figure 8. The evolution of infection compartments of humans and mosquitoes with Dh=0.1 and Dv=0. (a) The evolution of u1. (b) The evolution of u4. (c) The evolution of u2. (d) The evolution of u5.

6.3. The impact of medical resources on disease transmission

In this position, the influence of medical resources on malaria transmission is analysed from two aspects: (i) The impact of the amount of medical resources; (ii) For fixed medical resources, the influence of different allocation measures on malaria transmission. In this subsection, we take strain 1 as an example and draw conclusions about strain 2.

Since there are more medical resources in urban areas than in rural areas, people can get better medical services, and therefore have higher recovery rates. We already know that a1 varies between [0.042,0.57], here, take a1=0.085,0.15,0.25,0.35,0.45 to get the spatial distribution of the recovery rate in Figure . The corresponding medical resources in the whole region are 0πϱ1(x)dx=0.08925,0.1575,0.2625,0.3675,0.4725. In this paper, we use the recovery rate of the entire region 0πϱ1(x)dx to represent medical resources. Choose c1=0.09, α1=0.04, and other parameters are consistent with those in Table . Under different medical resources, through calculation, the corresponding R1 is 1.4724, 1.1484, 0.9347, 0.8251, 0.7556, respectively. By observation, in a heterogeneous environment, increasing medical resources will not have a great impact on improving the recovery rate in rural areas. However, the analysis in Subsection 6.2 shows that, due to human diffusion, some patients in rural areas will go to the cities for medical treatment and enjoy better medical care, so the overall transmission risk is reduced.

Figure 9. Increasing medical resources in the environment of spatial heterogeneity and its influence on the R1.

Figure 9. Increasing medical resources in the environment of spatial heterogeneity and its influence on the R1.

First of all, we hope to see whether maintaining a balanced distribution of urban and rural medical resources will help control the disease for a fixed recovery rate. When fixed medical resources, exploring the effect of different allocation on malaria spread. Fix a1=0.085 and introduce a new parameter σ6 into ϱ1, that is, ϱ1(x)=0.085×(1.05σ6cos(2x)), where σ6[0,1], to measure the difference in medical resources allocation between urban and rural areas. As σ6=1, there is the largest gap between urban and rural medical resources, and the urban areas have the highest recovery rates. With the change of σ6 from 1 to 0, medical resources are delivered to nearby rural areas, and the gap in medical conditions between rural and urban areas gradually is narrowed, and finally evenly distributed in space. Assume that the treatment capacity of the entire medical system is certain, 0πϱ1(x)dx=0.28, and dose not vary with σ6. Select c1=0.045, α1=0.04, and other parameters are kept in Table . By calculation, when σ6 takes 0,14,12,34,1, R1 is 0.9879, 0.9911, 1.0006, 1.0171, 1.0411, respectively. Looking at Figure , the spatial heterogeneity of medical resources distribution increases the malaria spread risk, which demonstrates that the distribution of medical resources dose act a critical purpose in the pattern of malaria control. By calculation, if σ6=1, in order to make R1 less than 1, medical resources need at least 0.2205, but if σ6=0, the medical resources only need 0.1733, which can make R1<1. In this way, we come to an interesting conclusion that under the condition 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.

Figure 10. For fixed medical resources, the spatial distribution of ϱ1 and the corresponding value of R1.

Figure 10. For fixed medical resources, the spatial distribution of ϱ1 and the corresponding value of R1.

6.4. The impact of seasonal changes in temperature and vector-bias on disease transmission

Here, we are interested in the impact of temperature-dependent EIP and vector-bias on the spread of disease. Now take single-strain subsystem (Equation10) as an example.

6.4.1. The role of vector-bias

To understand the influence of vector-bias on the size of malaria transmission, we make two comparison figures. The parameter value is consistent with the parameter value of Case 1, in this case, p = 0.8 and l = 0.6, i.e. l/p=3/4, Figure (a,c), show the evolutions of infected mosquitoes and individuals at each position in [0,π] as time evolves. Adjust the values of p = 0.8 and l = 0.8, so that l/p=1, that is, there is no vector-bias effect, shown in Figure (b,d). Through observation and comparison, if vector-bias effect is not taken into account, the duration of disease extinction will be shortened, i.e. if there is no vector-bias effect, the risk of malaria transmission will be underestimated.

Figure 11. The evolution of infection compartments of humans and mosquitoes with different vector-bias parameter. Among them, green indicates that there is a vector-bias effect (l/p=3/4), and red indicates that there is no vector-bias effect (l/p=1). (a) With vector-bias effect. (b) Without vector-bias effect. (c) With vector-bias effect. (d) Without vector-bias effect.

Figure 11. The evolution of infection compartments of humans and mosquitoes with different vector-bias parameter. Among them, green indicates that there is a vector-bias effect (l/p=3/4), and red indicates that there is no vector-bias effect (l/p=1). (a) With vector-bias effect. (b) Without vector-bias effect. (c) With vector-bias effect. (d) Without vector-bias effect.

6.4.2. The impact of seasonal changes in temperature on disease transmission

To understand the influence of temperature-dependent EIP on malaria transmission, two comparative figures of evolution of infectious humans and mosquitoes are produced, [τ1]:=1ω0ωτ1(t)dt=17.25/30.4Month, [β]:=1ω0ωβ(t)dt=6.983Month1, [μv]:=1ω0ωμv(t)dt=3.086Month1, [Λ]:=1ω0ωΛ(t)dt=k^×[β]. The parameters are constant, that is, they are not affected by temperature, and the other parameter values are consistent with Case 1.

Figure  reveals that if we do not consider seasonal changes in temperature when analysing the spread of malaria, the duration of disease extinction will be shortened, that is, the risk of disease transmission will be underestimated. This will lead to deviations in the prevention and control departments when formulating measures to control and eliminate malaria.

Figure 12. The evolution infection compartments of humans and mosquitoes with seasonal temperature changes and without seasonal temperature changes. Among them, red indicates that the parameters depend on temperature, and green indicates that the parameters do not depend on temperature.

Figure 12. The evolution infection compartments of humans and mosquitoes with seasonal temperature changes and without seasonal temperature changes. Among them, red indicates that the parameters depend on temperature, and green indicates that the parameters do not depend on temperature.

7. Discussion

This paper proposes and analyses a reaction–diffusion malaria model with vector-bias effect and multi-strains in a periodic environment. With the purpose of exploring the dynamics, we first use the existing theory to derive Ri, R0, and then introduce R^i ((i=1,2)). Using the comparison arguments and persistence theory of periodic semi-flow, we prove that malaria disease will be eliminated if R1<1 and R2<1, and the strain i is extinct and the strain j is persistent if Rj>1>Ri (i,j=1,2,ij). When R1>1 and R2>1, there are three cases: (i) strain 1 is die out and strain 2 persists; (ii) strain 1 is persistent and strain 2 is extinct; (iii) the two strains coexist under the condition of R^1>1, R^2>1 and exhibit spatial and seasonal fluctuations. This phenomenon is consistent with the multi-strain coexistence problem considered in Ref. [Citation30,Citation46], but they do not consider the periodic delay. Wu and Zhao [Citation40] propose a reaction–diffusion model of vector-borne disease with periodic delays to investigate the multiple effects of the temperature sensitivity of incubation periods, spatial heterogeneity, and the seasonality on disease transmission. However, they did not consider the case of multiple strains.

The numerical analysis is divided into four parts. First, the simulation results support the theoretical analysis. Second, the impact of diffusion on malaria transmission in a heterogeneous space is explored. Biologically, we figure out this spatial heterogeneity to be the difference between rural and urban regions. Simulation results show that under the condition of uneven distribution of medical resources, ignoring the diffusion of humans and mosquitoes, the risk of malaria transmission will be overestimated, which is consistent with the result in Ref. [Citation40,Citation47]. The implications include: (i) Infected humans in rural areas can seek better treatment in cities. (ii) it is challenging for mosquitoes to pick up a blood meal with moving humans. (iii) The impact of medical resources allocation on the spread of malaria is explored from two aspects: (a) increasing medical resources; (b) fixed medical resources, different allocation measures. Third, through research, we see that in the case of uneven distribution of medical resources, increasing medical resources has little impact on the spread of disease in rural areas, but it will reduce the risk of disease transmission in the entire region. Besides, when medical resources is limited, reducing the gap in the allocation of medical resources between rural and urban areas can further decrease the malaria spread risk, which is consistent with the result in Ref. [Citation40]. Therefore, when formulating measures to eliminate malaria, attention should be paid to the input of rural medical resources. Numerical results indicate that the risk of malaria may be underestimated if ignoring vector-bias effect or season. [Citation1] established a single strain ordinary differential equation model to study the effect of vector-bias on malaria transmission. The result taken by Abboubakar et al. [Citation1] in Figure  shows that increasing vector-bias parameter decreases the equilibrium prevalence of infectious humans which is consistent with our conclusion.

In our study, there are some shortcomings that need to be improved. In the process of theoretical analysis, due to the complexity of periodic delays, we have not discussed all situations clearly. We have not solved the global stability of the positive periodic solution of this model. In malaria transmission model, few researchers have studied the influence of infection age and spatial diffusion on disease transmission. In fact, after malaria infection, the intensity of infectivity varies at different stages of infection. The time after infection is called the age of infection, affects the number of secondary infections. This important factor needs to be considered in the process of malaria transmission. These problems are very important for exploring the law of malaria transmission and formulating control measures. In future work, we hope to better explore them.

Acknowledgments

The authors are grateful for the support and useful comments provided by Professor Xiaoqiang Zhao (Memorial University of Newfoundland).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This research is supposed by National Natural Science Foundation of China [grant number 11971013] and the Postgraduate Research and Practice Innovation Program of Jiangsu Province [grant number KYCX21_0176].

References

  • H. Abboubakar, B. Buonomo, and N. Chitnis, Modelling the effects of malaria infection on mosquito biting behaviour and attractiveness of humans, Ric. Mat. 65(1) (2016), pp. 329–346.
  • Z. Bai, R. Peng, and X.-Q. Zhao, A reaction–diffusion malaria model with seasonality and incubation period, J. Math. Biol. 77 (2018), pp. 201–228.
  • H.J. Bremermann and H.R. Thieme, A competitive exclusion principle for pathogen virulence, J. Math. Biol. 27 (1989), pp. 179–190.
  • R.S. Cantrell and C. Cosner, Spatial Ecology via Reaction–Diffusion Equations, John Wiley and Sons Ltd, Chichester, England, 2003.
  • F. Chamchod and N.F. Britton, Analysis of a vector-bias model on malaria transmission, Bull. Math. Biol. 73(3) (2011), pp. 639–657.
  • N. Chitnis, J.M. Hyman, and J.M. Cushing, Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model, Bull. Math. Biol. 70(5) (2008), pp. 1272–1296.
  • A.T. Ciota, A.C. Matacchiero, A.M. Kilpatrick, and L.D. Kramer, The effect of temperature on life history traits of Culex mosquitoes, J. Med. Entomol. 51(1) (2014), pp. 55–62.
  • O. Diekmann, J.A.P. Heesterbeek, and J.A.J. Metz, On the definition and the computation of the basic reproduction ratio R0 in models for infectious-diseases in heterogeneous populations, J. Math. Biol. 28(4) (1990), pp. 365–382.
  • D.A. Ewing, C.A. Cobbold, B.V. Purse, M.A. Nunn, and S.M. White, Modelling the effect of temperature on the seasonal population dynamics of temperate mosquitoes, J. Theor. Biol. 400 (2016), pp. 65–79.
  • Z. Feng, Z. Qiu, Z. Sang, C. Lorenzo, and J. Glasser, Modeling the synergy between HSV-2 and HIV and potential impact of HSV-2 therapy, Math. Biosci. 245(2) (2013), pp. 171–187.
  • J. Ge, K.I. Kim, Z. Lin, and H. Zhu, A SIS reaction–diffusion-advection model in a low-risk and high-risk domain, J. Differ. Equ. 259 (2013), pp. 5486–5509.
  • M.A. Gregson, Communicable disease control in emergencies: a field manual, Ann. Internal Med. 146(7) (2007), pp. 544–544.
  • T.J. Hagenaars, C.A. Donnelly, and N.M. Ferguson, Spatial heterogeneity and the persistence of infectious diseases, J. Theor. Biol. 229(3) (2004), pp. 349–359.
  • P. Hess, Periodic-Parabolic Boundary Value Problems and Positivity, Longman Scientific and Technical, Harlow, UK, 1991.
  • Y. Jin and X.-Q. Zhao, Spatial dynamics of a nonlocal periodic reaction–diffusion model with stage structure, SIAM J. Math. Anal. 40(6) (2009), pp. 2496–2516.
  • P. Koch-Mesina and D. Daners, Abstract Evolution Equations, Periodic Problems and Applications, Longman Scientific and Technical, Harlow, UK, 1992.
  • M. Kot, Elements of Mathematical Ecology, Cambridge University Press, Cambridge, UK, 2001.
  • R. Lacroix, W. R Mukabana, L. C. Gouagna, J. C Koella, and T. Smith, Malaria infection increases attractiveness of humans to mosquitoes, PLoS Biol. 3(9) (2005), Article ID e298.
  • X. Liang, L. Zhang, and X.-Q. Zhao, Basic reproduction ratios for periodic abstract functional differential equations (with application to a spatial model for Lyme disease), J. Dyn. Differ. Equ. 31 (2019), pp. 1247–1278.
  • X. Liu, Y. Wang, and X.-Q. Zhao, Dynamics of a climate-based periodic Chikungunya model with incubation period, Appl. Math. Model. 80 (2020), pp. 151–168.
  • Y. Lou and X.-Q. Zhao, A reaction–diffusion malaria model with incubation period in the vector population, J. Math. Biol. 62(4) (2011), pp. 543–568.
  • Y. Lou and X.-Q. Zhao, A theoretical approach to understanding population dynamics with seasonal developmental durations, J. Nonlinear Sci. 27 (2017), pp. 573–603.
  • G. Macdonald, The analysis of equilibrium in malaria, Trop. Dis. Bull. 49(9) (1952), pp. 813–829.
  • P. Magal and X.-Q. Zhao, Global attractors and steady states for uniformly persistent dynamical systems, SIAM J. Math. Anal. 37(1) (2005), pp. 251–275.
  • R.H. Martin and H.L. Smith, Abstract functional differential equations and reaction–diffusion systems, Trans. Am. Math. Soc. 321(1) (1990), pp. 1–44.
  • E.T. Ngarakana-Gwasira, C.P. Bhunu, and E. Mashonjowa, Assessing the impact of temperature on malaria transmission dynamics, Afr. Mat. 25(4) (2014), pp. 1095–1112.
  • R.M. Nisbet and W.S.C. Gurney, The systematic formulation of population models for insects with dynamically varying instar duration, Theor. Popul. Biol. 23(1) (1983), pp. 114–135.
  • R. Ross, The Prevention of Malaria, John Murray, London, 1911.
  • S. Ruan, D. Xiao, and J.C. Beier, On the delayed Ross–Macdonald model for malaria transmission, Bull. Math. Biol. 70(4) (2008), pp. 1098–1114.
  • Y. Shi and H. Zhao, Analysis of a two-strain malaria transmission model with spatial heterogeneity and vector-bias, J. Math. Biol. 82(4) (2021), pp. 1-42.
  • H.L. Smith, Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems, American Mathematical Society, Providence, 1995.
  • W.J. Tabachnick, Challenges in predicting climate and environmental effects on vector-borne disease episystems in a changing world, J. Exp. Biol. 213 (2010), pp. 946–954.
  • H.R. Thieme, Spectral bound and reproduction number for infinite-dimensional population structure and time heterogeneity, SIAM J. Appl. Math. 70 (2009), pp. 188–211.
  • C. Vargas-De-León, Global analysis of a delayed vector-bias model for malaria transmission with incubation period in mosquitoes, Math. Biosci. Eng. 9(1) (2011), pp. 165–174.
  • X. Wang and X.-Q. Zhao, A malaria transmission model with temperature-dependent incubation period, Bull. Math. Biol. 79 (2017), pp. 1155–1182.
  • X. Wang and X.-Q. Zhao, A periodic vector-bias malaria model with incubation period, SIAM J. Appl. Math. 77(1) (2017), pp. 181–201.
  • WHO, Available at: Malaria, http://www.who.int/ith/ith_country_list.pdf
  • WHO, (2018). Availabl at: World Malaria Report, https://www.who.int/malaria/media/world-malaria-report-2018/zh/
  • J. Wu, Theory and Applications of Partial Functional Differential Equations, Springer, New York, 1996.
  • R. Wu and X.-Q. Zhao, A reaction–diffusion model of vector-borne disease with periodic delays, J. Nonlinear Sci. 29 (2019), pp. 29–64.
  • Y. Xiao and X. Zou, Transmission dynamics for vector-borne diseases in a patchy environment, J. Math. Biol. 69(1) (2014), pp. 113–146.
  • D. Xu and X.-Q. Zhao, Dynamics in a periodic competitive model with stage structure, J. Math. Anal. Appl. 311 (2005), pp. 417–438.
  • L. Zhang, Z. Wang, and X.-Q. Zhao, Threshold dynamics of a time periodic reaction–diffusion epidemic model with latent periodic, J. Differ. Equ. 258(9) (2015), pp. 3011–3036.
  • X.-Q. Zhao, Basic reproduction ratios for periodic compartmental models with time delay, J. Dyn. Differ. Equ. 29 (2017), pp. 67–82.
  • X.-Q. Zhao, Dynamical Systems in Population Biology, Springer-Verlag, London, 2017.
  • L. Zhao, Z.-C. Wang, and S. Ruan, Dynamics of a time-periodic two-strain SIS epidemic model with diffusion and latent period, Nonlinear Anal.: Real World Appl. 51 (2020), Article ID 102966.
  • H. Zhao, Y. Shi, and X. Zhang, Dynamic analysis of a malaria reaction–diffusion model with periodic delays and vector bias, Math. Biosci. Eng. 19(3) (2022), pp. 2538–2574.
  • T. Zheng, L.-F. Nie, Z. Teng, and Y. Luo, Competitive exclusion in a multi-strain malaria transmission model with incubation period, Chaos Solitons Fractals 131 (2020), Article ID 109545.

Appendices

 

Appendix 1. Derivation of Mvi (i = 1, 2)

Once again, ρ(t,q,x)q is the number of female mosquitoes of development level infection (q,q+q) at time t and location x. If we consider the rate of change of the number of female mosquitoes in a given development level of infection interval (q,q+q), we may write t[Ωρ(t,q,x)qdV]=[+rateofentryatqrateofdepartureat(q+q)deathsmoveout],or tΩρ(t,q,x)qdV=Ωf(t,q,x)qdV+ΩJ1(t,q,x)dVΩJ1(t,q+q,x)dVΩJ2(t,q,x)qdS,where f(t,q,x) represents the rate of the population change per unit level of infection development, J1(t,q,x) denotes the positive (left to right) ‘flux’ of individuals at the level of infection development q, time t and position x, and J2(t,q,x) is the diffusive flux. dS is the vector element of surface. Dividing by q gives us tΩρ(t,q,x)dV=Ωf(t,q,x)dVΩJ1(t,q+q,x)J1(t,q,x)qdVΩJ2(t,q,x)dS.Use Gauss's divergence theorem to transform the third integral into a volume integral to simplify this equation, ΩJ2(t,q,x)dS=ΩJ2(t,q,x)dV.Taking the limit as q approaches zero produces a conservation law for the density of individuals (A1) tΩρ(t,q,x)dV=Ωf(t,q,x)dVΩJ1(t,q,x)qdVΩJ2(t,q,x)dV.(A1) This causes Equation (EquationA1) reducing to Ω[ρ(t,q,x)tf(t,q,x)+J1(t,q,x)q+J2(t,q,x)]dV=0.Since it is for arbitrary Ω, this integral is zero, then (A2) ρ(t,q,x)tf(t,q,x)+J1(t,q,x)q+J2(t,q,x)=0.(A2) For more clarity, the form of f(t,q,x) needs to be discussed. According to Ref. [Citation17], f(t,q,x) depends on the independent variables, the level of infection development, time and space. However, this is usually achieved through the dependent variable, population density. Here, let f(t,q,x)=μv(t,x)ρ(t,q,x). According to Fick's law and the hypothesis of the mosquito diffusion coefficient Dv, write the diffusive flux as J2(t,q,x)=Dvρ(t,q,x).General balance Equation (EquationA2) now reduces to (A3) ρ(t,q,x)t=μv(t,x)ρ(t,q,x)J1(t,q,x)q+Dv2ρ(t,q,x),(A3) where 2ρ(t,q,x) is the Laplacian, i.e. 2ρ(t,q,x)=Δρ(t,q,x). So (EquationA3) can write as (A4) ρ(t,q,x)t=DvΔρ(t,q,x)J1(t,q,x)qμv(t,x)ρ(t,q,x),(A4) To progress any more, we must say something about the flux of mosquitoes J1(t,x,q). This is not a flux in space but a development level of infection. All mosquito infection levels develop, and then J1(t,q,x)=γ(t)ρ(t,q,x).Write Equation (EquationA4) as (A5) ρ(t,q,x)t=DvΔρ(t,q,x)γ(t)ρ(t,q,x)qμvρ(t,q,x).(A5) Then by the same arguments as in Ref. [Citation40] and Ref. [Citation17], we can derive the expression of Mvi (i = 1, 2) as Mvi(t,x)=(1τi(t))ΩΓ(t,tτi(t),x,y)αiβ(tτi(t),y)pIi(tτi(t),y)Sv(tτi(t),y)p[I1(tτi(t),y)+I2(tτi(t),y)]+l[N(y)I1(tτi(t),y)I2(tτi(t),y)]dy.Γ(t,t0,x,y) is the Green function associated with u(t,x)t=DvΔu(t,x)μv(t,x)u(t,x), for tt00 and x,yΩ, subject to the Neumann boundary condition.

Appendix 2.

 

Proof of Lemma 4.5

Proof.

Similar to the proof of Lemma 4.4, on interval [nτ¯1,(n+1)τ¯1], nN, one can easily obtain that ϖi(t,)0, (i=1,2) for all t0. Choose a large number K>max{μ~h+ϱ~1,μ~v}, where μ~h=maxxΩ¯μh(x), μ~v=maxt[0,ω],xΩ¯μv(t,x) and ϱ~1=maxxΩ¯ϱ1(x), such that for each tR, g1(t,,ϖ1):=(μh()+ϱ1())ϖ1+Kϖ1 is increasing in ϖ1, and g2(t,,ϖ2):=μv(t,)ϖ2+Kϖ2 is increasing in ϖ2. Then, both ϖ1(t,x) and ϖ2(t,x) satisfy the followingsystem {ϖ1(t,x)t=DhΔϖ1(t,x)+c1β(t,x)ϖ2(t,x)+g1(t,x,ϖ1)Kϖ1(t,x),xΩ,ϖ2(t,x)t=DvΔϖ2(t,x)+g2(t,x,ϖ3)Kϖ2(t,x)+(1τ1(t))ΩΓ(t,tτ1(t),x,y)α1β(tτ1(t),y)pu3(tτ1(t),y)ϖ1(tτ1(t),y)lN(y)dy,xΩ,ϖ1(t,x)ν=ϖ2(t,x)ν=0,xΩ.Thus, for a given ξE1+, one has (A6) {ϖ1(t,,ξ)=T~1(t,0)ξ1(0)+0tT~1(t,s)g1(s,,ϖ1(s,))ds+0tT~1(t,s)c1β(s,)ϖ2(s,)ds,ϖ2(t,,ξ)=T~2(t,0)ξ2(0)+0tT~2(t,s)g2(s,,ϖ2(s,))ds+0tT~2(t,s)(1τ1(s))ΩΓ(s,sτ1(s),,y)α1β(sτ1(s),y)pu3(sτ1(s),y)ϖ1(sτ1(s),y)lN(y)dyds,(A6) where T~1(t,s),T~2(t,s):=YY are the evolution operators associated with ϖ1t=DhΔϖ1Kϖ1 and ϖ2t=DvΔϖ2Kϖ2 subject to the Neumann boundary condition, respectively. Since Lˇ(t):=tτ1(t) is increasing in tR, it easily follows that [τ1(0),0]Lˇ([0,τ^1]). We assume that φ1>0. Then there exists an (θ,x0)[τ1(0),0]×Ω such that ϖ1(θ,x0)>0. According to the second equation of (EquationA6), we have ϖ2(t,,φ)>0 for all t>τ^1. Note that if s>2τ^1, then sτ1(s)>2τ^1τ^1=τ^1. It is easy to know that ϖ1(t,,φ)>0 for all t>2τ^1 from the first equation of (EquationA6). This shows that ϖi(t,)>0 (i=1,2) for all t>2τ^1, therefore, the solution map P¯t is strongly positive whenever t>3τ^1.

Appendix 3.

 

Proof of Theorem 5.2

Proof.

(a) In this case of R1<1, it follows from Lemmas 4.2 and 4.6 that r(P¯1)<1, and therefore μ1=lnr(P¯1)ω<0. Discuss the following system with parameter ϵ>0 (A7) {wˇ1(t,x)t=DhΔwˇ1(t,x)+c1β(t,x)wˇ2(t,x)(μh(x)+ϱ1(x))wˇ1(t,x),t>0,xΩ,wˇ2(t,x)t=DvΔwˇ2(t,x)μv(t,x)wˇ2(t,x)+(1τ1(t))ΩΓ(t,tτ1(t),x,y)α1β(tτ1(t),y)p[u3(tτ1(t),y)+ϵ]wˇ1(tτ1(t),y)lN(y)dy,t>,xΩ,wˇ1(t,x)ν=wˇ2(t,x)ν=0,t>0,xΩ.(A7) Denote wˇ(t,x,φ)=(wˇ1(t,x,φ),wˇ2(t,x,φ)) to be the unique solution of system (EquationA7) for any φE1, with wˇ0(φ)(θ,x)=φ(θ,x) for all θ[τ1(0),0], xΩ¯, where wˇt(φ)(θ,x)=wˇ(t+θ,x,φ)=(wˇ1(t+θ,x,φ),wˇ2(t,x,φ)),t0,(θ,x)[τ1(0),0]×Ω¯.Let P¯1ϵ:=E1E1 be the Poincaré map of system (EquationA7), i.e. P¯1ϵ(φ)=wˇω(φ), φE1, and make r(P¯1ϵ) to be the spectral radius of P¯1ϵ. From the continuity of the principle eigenvalue about parameters, we know that limϵ0r(P¯1ϵ)=r(P¯1). Allow ϵ>0 to be small enough, one has r(P¯1ϵ)<1. Based on Lemma 4.7, there exists a positive ω-periodic function wˇ(t,x) such that wˇ(t,x)=eμ1ϵtwˇ(t,x) is a solution of system (EquationA7), where μ1ϵ=lnr(P¯1ϵ)ω<0. For ϵ>0 fixed above, by the global attractivity of u3(t,x) of system (Equation12) and the comparison principle, there must be an integer n1>0 large enough such that n1ω>τ^ and u¯2(t,x)u3(t,x)+ϵ,tn1ωτ^,xΩ¯.Thus, as t>n1ω, one has (A8) {u¯1(t,x)tDhΔu¯1(t,x)+c1β(t,x)u¯3(t,x)(μh(x)+ϱ1(x))u¯1(t,x),xΩ,u¯3(t,x)tDvΔu¯3(t,x)μv(t,x)u¯3(t,x)+(1τ1(t))ΩΓ(t,tτ1(t),x,y)α1β(tτ1(t),y)p[u3(tτ1(t),y)+ϵ]u¯1(tτ1(t),y)lN(y)dy,xΩ,u¯1(t,x)ν=u¯3(t,x)ν=0,xΩ.(A8) There exists some σ1>0, for any designated φ¯X4, such that (u¯1(t,x,φ¯),u¯3(t,x,φ¯))σ1(wˇ1(t,x),wˇ2(t,x)),t[n1ωτ^,n1ω],xΩ¯.Using (EquationA7), (EquationA8) and the comparison theorem for abstract functional differential equation [Citation25, Proposition 1], we get (u¯1(t,x,φ¯),u¯3(t,x,φ¯))σ1eμ1ϵtwˇ(t,x),tn1ω,xΩ¯.Accordingly, limt(u¯1(t,x,φ¯),u¯3(t,x,φ¯))=(0,0) uniformly for xΩ¯.

Using the theory of internal chain transitive sets [Citation45, Section 1.2], for xΩ¯, limt(u¯2(t,x,φ¯)u3(t,x))=0 uniformly, where u3(t,x) is a globally attractive solution of system (Equation12). Based on the above argument, it is easy to see that u¯2(t,x,φ) meets a nonautonomous system, which is asymptotic to the periodic system (Equation12).

Clearly, Dt and D are solution map and Poincaré map associated with system (Equation10) on X4, respectively. Let J=ω(φ¯) be the omega limit set of φ¯=(φ¯1,φ¯2,φ¯3)X4 for D. Since limtu¯1(t,x,φ)=0 and limtu¯3(t,x,φ¯)=0 uniformly for xΩ¯, one has J={0^}×J¯×{0}. Here, 0^ represents a function which is identical to 0, i.e. 0^(θ,)=0, θ[τ1(0),0]. Following Lemma 5.2, one has 0^J¯.

For any ζC([τ1(0),0],Y+), let w(t,x,ζ(0,)) be the solution of (Equation12) with initial data w(0,x)=ζ(0,x). A solution semi-flow of (Equation12) on C([τ1(0),0],Y+) described as wt(θ,x,ζ)={w(t+θ,x,ζ(0)),ift+θ>0,θ[τ1(0),0],t>0,ζ(t+θ,x),ift+θ0,θ[τ1(0),0],t>0.Set D¯(ζ)=wω(ζ). Following [Citation45, Lemma 1.2.1], J is an internal chain transitive set for D, and hence J¯ is an internal chain transitive set for D¯. Define u30C([τ1(0),0],Y+) by u30(θ,)=u3(θ,) for θ[τ1(0),0]. Due to J¯{0^} and u30 being globally attractive in C([τ1(0),0],Y+){0^}, we know J¯Ws(u30), where Ws(u30) is the stable set of u30. By [Citation45, Theorem 1.2.1], one has J¯={u30}, which leads to J={0^,u30,0}, thence, limt(u¯1(t,,φ¯),u¯2(t,,φ¯),u¯3(t,,φ¯))(0,u3(t,),0)X3=0.(b) As R1>1, Lemmas 4.2 and 4.6 indicate r(P¯1)>1, consequently, μ1=lnr(P¯1)ω>0. Let P0:={φ¯X4:φ¯1(0,)0andφ¯3()0},and P0:=X4P0={φ¯X4:φ¯1(0,)0orφ¯3()0}.For any φ¯P0, Lemma 5.2 reveals that u¯1(t,x,φ¯)>0 and u¯3(t,x,φ¯)>0, t>0, xΩ¯. Thus knowing Dn(P0)P0, nN. Lemma 5.1 evidences that D:X4X4 has a strong global attractor in X4.

Set M:={φ¯P0:Dn(φ¯)P0,nN},and ω(φ¯) to be the omega limit set of the orbit γ+:={Dn(φ¯):nN}. Denote Q=(0^,u3(t,x),0). Next, we prove that Q cannot form a cycle for D in P0.

Claim 1. For any φ¯M, the omega limit set ω(φ¯)=Q.

For any designated φ¯M, Dn(φ¯)P0, nN. Then, either u¯1(nω,,φ¯)0 or u¯3(nω,,φ¯)0, for each nN. Furthermore, by contradiction and Lemma 5.2, it is obvious that for each t0, either u¯1(t,,φ¯)0 or u¯3(t,,φ¯)0. If u¯1(t,,φ¯)0 for all t0, then Lemma 4.1 guarantees limtu¯2(t,x,φ¯)=u3(t,x) uniformly for xΩ¯. Thereby, the u¯3 equation in (Equation10) meets (A9) u¯3(t,x,φ¯)tDvΔu¯3(t,x,φ¯)μ¯vu¯3(t,x,φ¯),(A9) where μ¯v=mint[0,ω],xΩ¯μv(t,x). By the comparison principle, we have limtu¯3(t,x)=0 uniformly for xΩ¯. Suppose that u¯1(t2,,φ¯)0 for some t20, in the light of Lemma 5.2, we obtain u¯1(t,,φ¯)>0, tt2. Accordingly, u¯3(t,x)0, tt2. From the u¯1 equation in (Equation10), one has limtu¯1(t,x,φ¯)=0 uniformly for xΩ¯. Consequently, the u¯2 equation in (Equation10) abides by a nonautonomous system, which is asymptomatic to the periodic system (Equation12). The theory of asymptomatically periodic system [Citation45, Section 3.2] implies that limt(u¯2(t,x)u3(t,x))=0 uniformly for xΩ¯. Therefore, ω(φ¯)=Q for any φ¯M, and Q cannot form a cycle for D in P0.

Consider the following time-periodic parabolic system with parameter δ>0 (A10) {v1δ(t,x)t=DhΔv1δ(t,x)+c1β(t,x)l[N(x)δ]v2δ(t,x)pδ+lN(x)(μh(x)+ϱ1(x))v1δ(t,x),t>0,xΩ,v2δ(t,x)t=DvΔv2δ(t,x)μv(t,x)v2δ(t,x)+(1τ1(t))ΩΓ(t,tτ1(t),x,y)α1β(tτ1(t),y)p[u3(tτ1(t),y)δ]v1δ(tτ1(t),y)pδ+lN(y)dy,t>0,xΩ,v1δ(t,x)ν=v2δ(t,x)ν=0,t>0,xΩ.(A10) For any φE1, set vδ(t,x,φ)=(v1δ(t,x,φ),v2δ(t,x,φ)) to be the unique solution of system (EquationA10) with v0δ(φ)(θ,x)=(ϕ1(θ,x),ϕ2(0,x)) for all θ[τ1(0),0], xΩ¯, where vtδ(φ)(θ,x)=vδ(t+θ,x,φ)=(v1δ(t+θ,x,φ),v2δ(t,x,φ)),t0,(θ,x)[τ1(0),0]×Ω¯.Let P¯1δ:E1E1 be the Poincaré map of system (EquationA10), i.e. P¯1δ(φ)=vωδ(φ), φE1. Make r(P¯1δ) to be the spectral radius of P¯1δ. Since limδ0r(P¯1δ)=r(P¯1)>1, choose a sufficiently small δ>0 such that δ<min{mint[0,ω],xΩ¯u3(t,x),minxΩ¯N(x)},andr(P¯δ)>1.For the above fixed δ>0, according to the continuous dependence of the solution on the initial data, there must be a constant δ>0 such that for all φ¯ with φ¯Q<δ, which leads to Dt(φ¯)Dt(Q)<δ for all t[0,ω]. We now prove the following claim.

Claim 2. For any φ¯P0, lim supnDn(φ¯)Qδ.

Assuming contradictions, for some φ¯0P0, one has limsupnDn(φ¯0)Q<δ. So there is n21 such that Dn(φ¯0)Q<δ for all nn2. For each tn2ω, setting t=nω+t with n=[t/ω] and t[0,ω), one has (A11) Dt(φ¯0)Dt(Q)=Dt(Dn(φ¯0))Dt(Q)<δ.(A11) Following (EquationA11) and Lemma 5.2, we can receive u3(t,x)δ<u¯2(t,x,φ¯0)and0<u¯i(t,x,φ¯0)<δ,i=1,3,for each tn2ωτ^, and xΩ¯. Thus, as tn2ω, u¯1(t,x,φ¯0) and u¯3(t,x,φ¯0) satisfy (A12) {u¯1(t,x)tDhΔu¯1(t,x)+c1β(t,x)l[N(x)δ]u¯3(t,x)pδ+lN(x)(μh(x)+ϱ1(x))u¯1(t,x),xΩ,u¯3(t,x)tDvΔu¯3(t,x)μv(t,x)u¯3(t,x)+(1τ1(t))ΩΓ(t,tτ1(t),x,y)α1β(tτ1(t),y)p[u3(tτ1(t),y)δ]u¯1(tτ1(t),y)pδ+lN(y)dy,xΩ,u¯1(t,x)ν=u¯3(t,x)ν=0,xΩ.(A12) Based on u¯(t,x,φ¯0)0 for all t>0 and xΩ¯, there must be a constant σ2>0 such that (u¯1(t,x,φ¯0),u¯3(t,x,φ¯0))σ2eμ1δtvδ(t,x),t[n2ωτ^,n2ω],xΩ¯,where vδ(t,x) is a positive ω-periodic function such that eμ1δtvδ(t,x) is a solution of system (EquationA10), and μ1δ=lnr(P¯1δ)ω>1. It follows from (EquationA12) and the comparison theorem that (u¯1(t,x,φ¯0),u¯3(t,x,φ¯0))σ2eμ1δtvδ(t,x),tn2ω,xΩ¯.Obviously, u¯i(t,x,φ¯0)+ as t+, i = 1, 3, which leads to a contradiction.

With the above claim, it is easy to obtain that Q is an isolated invariant set for D in X4, and Ws(Q)P0=, where Ws(Q) is the stable set of Q for D. For any integer n with nω>τ^, Dn:=Dnω is compact, it follows that D is compact, then D is asymptotically smooth on X4. Moreover, by Lemma 5.1, we obtain that D has a global attractor on X4. According to [Citation24, Theorem 3.7], D has a global attractor A0 in P0. Based on the acyclicity theorem on uniform persistence for maps [Citation45, Theorem 1.3.1, Remark 1.3.1], then D is uniformly persistent with respect to (P0,P0) in the sense that there exists an γ2>0, such that (A13) lim infnd(Dn(φ¯),P0)γ2,φ¯P0.(A13) Now we derive the practical persistence required. Following A0=D(A0), then φ¯1(0,)>0 and φ¯3()>0 for any φ¯A0. Set B:=t[0,ω]Dt(A0). Obviously, BP0 and limtd(Dt(φ¯),B)=0, φ¯P0. Define a continuous function q1:X4R+ by q1(φ¯)=min{minxΩ¯φ¯1(0,x),minxΩ¯φ¯3(x)},φ¯=(φ¯1,φ¯2,φ¯3)X4.Because B is a compact subset of P0, we obtain that infφ¯Bq1(φ¯)=minφ¯Bq1(φ¯)>0. Consequently, there is an γ3>0 such that lim inftq1(Dt(φ¯))=lim inftmin(minxΩ¯u¯1(t,x,φ¯),minxΩ¯u¯3(t,x,φ¯))γ3,φ¯P0.Furthermore, according to Lemma 5.2, there must be an γ4(0,γ3) such that lim inftminxΩ¯u¯j(t,x,φ¯)γ4,φ¯P0(j=1,2,3).The existence of a positive periodic steady state remains to be proved. By [Citation2, Lemma 8] and [Citation45, Theorem 3.5.1], for each t>0, the solution map Qˇt of system (Equation10), is a κ-contraction with respect to an equivalent norm on X3, where κ is Kuratowski measure of noncompactness. Define W0={φ¯X3:φ¯1(0,)0andφ¯3()0},and W0=X3/W0={φ¯X3:φ¯1(0,)0orφ¯3()0}.For any φ¯W0, Qˇ is ρ~1-uniformly persistent with ρ~1(φ¯)=d(φ¯,W0) is easily attainable. It then follows from [Citation24, Theorem 4.5], as applied to Qˇ, that system (Equation10) has an ω-periodic solution (z1(t,),z2(t,),z3(t,)) with (z1t,z2t,z3t)W0. Let u¯1(θ,)=z1(θ,), u¯2(θ,)=z2(θ,) and u¯3(0,)=z3(0,), where θ[τ1(0),0]. Combing the uniqueness of solutions and Lemma 5.2, we have that (u¯1(t,),u¯2(t,),u¯3(t,)) is periodic solution of system (Equation10), and it is also strictly positive.

Appendix 4.

 

Proof of Lemma 5.2

Proof.

According to the comparison principle for cooperative systems, we have u¯1(t,x,ψ^)0, and u¯3(t,x,ψ^)0 for all t>0 and xΩ¯. Moreover, u¯1(t,x,ψ^) and u¯3(t,x,ψ^) satisfy {u¯1(t,x)tDhΔu¯1(t,x)(μ~h+ϱ~1)u¯1(t,x),t>0,xΩ,u¯3(t,x)tDvΔu¯3(t,x)μ~vu¯3(t,x),t>0,xΩ,u¯1(t,x)ν=u¯3(t,x)ν=0,t>0,xΩ,where μ~h, ϱ~1, μ~v are defined in Section 3. If there exists t30, such that u¯1(t3,,ψ^)0 and u¯3(t3,,ψ^)0, one has u¯1(t,,ψ^)>0 and u¯3(t,,ψ^)>0, t>t3, based on the strong maximum principle [Citation14, Proposition 13.1]. Set n¯(t,x,ψ^) to be the solution of (A14) {n¯(t,x)t=DvΔn¯(t,x)+Λ(t,x)(α1β(t,x)pl+μv(t,x))n¯(t,x),t>0,xΩ,n¯(t,x)ν=0,t>0,xΩ.(A14) Clearly, Λ(t,x) is Hölder continuous and nonnegative nontrivial on R×Ω¯. By the comparison principle, u¯2(t,x)n¯(t,x)>0,t>0,xΩ¯.In addition, from [Citation43, Lemma 2.1], system (EquationA14) admits a unique positive ω-periodic n¯(t,), which is globally attractive. Then lim inftu¯2(t,x,ψ^3)ε2:=inft[0,ω],xΩ¯n¯(t,x),uniformly for xΩ¯. The proof is complete.

Appendix 5.

 

Parameter definition and value

Wang et al. [Citation36] have studied the seasonal factors that affect malaria transmission, including seasonal forced bite rate (A15) β(t)=6.9831.993cos(πt/6)0.4247cos(πt/3)0.128cos(πt/2)0.04095cos(2πt/3)+0.0005486cos(5πt/6)1.459sin(πt/6)0.007642sin(πt/3)0.0709sin(πt/2)+0.05452sin(2πt/3)0.06235sin(5πt/6)Month1,(A15) periodic mortality rate of mosquitoes (A16) μv(t)=3.086+0.04788cos(πt/6)+0.01942cos(πt/3)+0.007133cos(πt/2)+0.0007665cos(2πt/3)0.001459cos(5πt/6)+0.02655sin(πt/6)+0.01819sin(πt/3)+0.01135sin(πt/2)+0.005687sin(2πt/3)+0.003198sin(5πt/6)Month1,(A16) periodic EIP (A17) τ1(t)=1/30.4[17.25+8.369cos(πt/6)+4.806sin(πt/6)+3.27cos(πt/3)+2.857sin(πt/3)+1.197cos(πt/2)+1.963sin(πt/2)+0.03578cos(2πt/3)+1.035sin(2πt/3)0.3505cos(5πt/6)+0.6354sin(5πt/6)0.3257cos(πt)+0sin(πt)]Month,(A17) (A18) τ2(t)=0.8/30.4[17.25+8.369cos(πt/6)+4.806sin(πt/6)+3.27cos(πt/3)+2.857sin(πt/3)+1.197cos(πt/2)+1.963sin(πt/2)+0.03578cos(2πt/3)+1.035sin(2πt/3)0.3505cos(5πt/6)+0.6354sin(5πt/6)0.3257cos(πt)+0sin(πt)]Month,(A18) and mosquitoes recruitment rate (A19) Λ(t)=k^×β(t)(km2Month)1,wherek^=5×1205709.(A19) In addition, due to the assumption of vector-bias effect, and according to Ref. [Citation35], in the following simulations, we take p = 0.8 and l = 0.6.

In order to better understand the periodicity and biological significance of the above parameters, we make the following figures. Looking at Figure , the distribution and seasonal dynamics of mosquito populations are affected by climate change. Analysing Figure (a), the peak of mosquito bites happens during warm or rainy seasons (May to September). In a year, the death rate may be lower in spring and summer when the humidity and temperature are suitable for species growth, and higher in winter due to low temperature and dry weather. In addition, higher temperature can shorten the duration of EIP.

Figure A1. A small interval of development level of infection.

Figure A1. A small interval of development level of infection.

Figure A2. The parameters change with time. (a) Change of contact rate β(t) with time. (b) Change of EIP τ1(t) with time. (c) Change of the death rate of mosquitoes μv(t) with time.

Figure A2. The parameters change with time. (a) Change of contact rate β(t) with time. (b) Change of EIP τ1(t) with time. (c) Change of the death rate of mosquitoes μv(t) with time.

Table F1. Definition and value of parameters.