Loading [MathJax]/jax/element/mml/optable/GeneralPunctuation.js
Research article

A novel analytical treatment for the Ambartsumian delay differential equation with a variable coefficient

  • Received: 17 November 2024 Revised: 06 December 2024 Accepted: 09 December 2024 Published: 24 December 2024
  • MSC : 34K06, 65L03

  • The Ambartsumian delay differential equation with a variable coefficient is considered in this paper. An effective transformation is produced to convert the extended Ambartsumian equation to the pantograph model. Two kinds of analytical solutions are determined. The first solution is expressed as an exponential function multiplied by an infinite power series. The second solution is obtained as an infinite series in terms of exponential functions. Several exact solutions are established for different forms of the extended Ambartsumian equation under specific relations. In addition, the convergence analysis is addressed theoretically. Moreover, numeric calculations are conducted to estimate the accuracy. The results reveal that the present analysis is efficient and accurate and can be further applied to similar delay models in a straightforward manner.

    Citation: Rana M. S. Alyoubi, Abdelhalim Ebaid, Essam R. El-Zahar, Mona D. Aljoufi. A novel analytical treatment for the Ambartsumian delay differential equation with a variable coefficient[J]. AIMS Mathematics, 2024, 9(12): 35743-35758. doi: 10.3934/math.20241696

    Related Papers:

    [1] Adel M. Al-Mahdi, Mohammad M. Al-Gharabli, Nasser-Eddine Tatar . On a nonlinear system of plate equations with variable exponent nonlinearity and logarithmic source terms: Existence and stability results. AIMS Mathematics, 2023, 8(9): 19971-19992. doi: 10.3934/math.20231018
    [2] Adel M. Al-Mahdi . The coupling system of Kirchhoff and Euler-Bernoulli plates with logarithmic source terms: Strong damping versus weak damping of variable-exponent type. AIMS Mathematics, 2023, 8(11): 27439-27459. doi: 10.3934/math.20231404
    [3] Mohammad M. Al-Gharabli, Adel M. Al-Mahdi, Mohammad Kafini . Global existence and new decay results of a viscoelastic wave equation with variable exponent and logarithmic nonlinearities. AIMS Mathematics, 2021, 6(9): 10105-10129. doi: 10.3934/math.2021587
    [4] Sarra Toualbia, Abderrahmane Zaraï, Salah Boulaaras . Decay estimate and non-extinction of solutions of p-Laplacian nonlocal heat equations. AIMS Mathematics, 2020, 5(3): 1663-1679. doi: 10.3934/math.2020112
    [5] Salim A. Messaoudi, Mohammad M. Al-Gharabli, Adel M. Al-Mahdi, Mohammed A. Al-Osta . A coupled system of Laplacian and bi-Laplacian equations with nonlinear dampings and source terms of variable-exponents nonlinearities: Existence, uniqueness, blow-up and a large-time asymptotic behavior. AIMS Mathematics, 2023, 8(4): 7933-7966. doi: 10.3934/math.2023400
    [6] Adel M. Al-Mahdi, Mohammad M. Al-Gharabli, Maher Nour, Mostafa Zahri . Stabilization of a viscoelastic wave equation with boundary damping and variable exponents: Theoretical and numerical study. AIMS Mathematics, 2022, 7(8): 15370-15401. doi: 10.3934/math.2022842
    [7] Abdelbaki Choucha, Salah Boulaaras, Asma Alharbi . Global existence and asymptotic behavior for a viscoelastic Kirchhoff equation with a logarithmic nonlinearity, distributed delay and Balakrishnan-Taylor damping terms. AIMS Mathematics, 2022, 7(3): 4517-4539. doi: 10.3934/math.2022252
    [8] Adel M. Al-Mahdi, Mohammad M. Al-Gharabli, Mohamed Alahyane . Theoretical and numerical stability results for a viscoelastic swelling porous-elastic system with past history. AIMS Mathematics, 2021, 6(11): 11921-11949. doi: 10.3934/math.2021692
    [9] Shuai Li, Tianqing An, Weichun Bu . Existence results for Schrödinger type double phase variable exponent problems with convection term in RN. AIMS Mathematics, 2024, 9(4): 8610-8629. doi: 10.3934/math.2024417
    [10] Yousif Altayeb . New scenario of decay rate for system of three nonlinear wave equations with visco-elasticities. AIMS Mathematics, 2021, 6(7): 7251-7265. doi: 10.3934/math.2021425
  • The Ambartsumian delay differential equation with a variable coefficient is considered in this paper. An effective transformation is produced to convert the extended Ambartsumian equation to the pantograph model. Two kinds of analytical solutions are determined. The first solution is expressed as an exponential function multiplied by an infinite power series. The second solution is obtained as an infinite series in terms of exponential functions. Several exact solutions are established for different forms of the extended Ambartsumian equation under specific relations. In addition, the convergence analysis is addressed theoretically. Moreover, numeric calculations are conducted to estimate the accuracy. The results reveal that the present analysis is efficient and accurate and can be further applied to similar delay models in a straightforward manner.



    Malaria is a lethal infectious disease caused by Plasmodium parasites, which infects humans through the bite of infected female geneus Anopheles mosquitoes. It is endemic in about 100 countries and an estimated 3.3 billion people are at risk of being infected in the word [1]. Therefore, it is necessary to study the spread of malaria in the population. The pioneering work was given by Ross [2] which showed that malaria is transmitted by mosquitoes and proposed the first mathematical model of this disease. Subsequently, Macdonald promoted this [3], which plays a particularly important role in understanding the disease. However, the Ross-Macdonald model is extremely simplified, ignoring many critical factors of real-world ecology and epidemiology [4]. One omission is spatial heterogeneity, for example, the distribution of urban and rural areas and the spatial movement of mosquitoes and people, which increases the spatial transmission of malaria. Hence, in the malaria epidemiological model, spatial heterogeneity must be considered [5,6,7]. In 1951, Skellam specially examined the reaction-diffusion model of species population density in a bounded habitat [8]. One advantage of reaction-diffusion models is that they can easily incorporate simple rules about edge effects. Use these features of reaction-diffusion model to study how environmental heterogeneity affects malaria [9,10,11]. However, in these studies, it seems that the following important biological aspects are rarely paid, simultaneously.

    First, "vector-bias" is a typical feature of malaria transmission, which describes the different probabilities of a mosquito bitting humans. Lacroix et al. [12] proved that mosquitoes are generally more inclined to bite infected individuals. To describe this preference, Chamchod and Britton incorporated the vector-bias effect into malaria model [13]. Subsequently, a wealth of papers take into account the effect of vector-bias when studying malaria transmission [14,15,16]. Due to the different densities of mosquitoes in various regions, the infectivity of malaria is also different. The analysis of vector-bias effects is of great significance for understanding the dynamics of malaria transmission.

    Secondly, in recent years, climate has become more fluctuating [17], and the rise in temperature is related to global warming. Climate warming is directly connected with the number of Anopheles, the main vector of malaria. Because temperature affects the immature development rate, immature stage survival rate, adult size, adult lifespan, blood feeding, and reproduction ability of mosquitoes [18]. It is widely accepted that both the life cycle of mosquitoes and the development of parasites are extremely sensitive to temperature. Empirical evidence indicates that the life cycle of mosquitoes experiences four different stages: eggs, larvae, pupae and adults. The first three stages are also called aquatic or immature stages [19]. In addition, a very important process in the development of parasites is the extrinsic incubation period (EIP), i.e., mosquitoes may not be able to transmit the disease to humans for a period of time after taking infected blood meal. The lifespan of female adult mosquitoes is generally 3 to 100 days, and EIP can reach 30 days [20]. These infected mosquitoes that live pass EIP will keep infectious until they die. There is considerable evidence that temperature rise will shorten the length of each immature stage [21] and EIP is extremely sensitive to seasonal changes in temperature [22,23]. In view of climate change, exploring the effect of temperature in the spread of malaria is particularly essential. Therefore, these seasonal forced EIP and mosquito life cycles need to be considered in the malaria model.

    Note that mathematical models of malaria transmission rarely consider key features such as vector-bias, climate/seasonality, temperature-dependent mosquito life cycle and EIP in a spatially heterogeneous environment, simultaneously. A reaction-diffusion model with periodic delays is proposed to study the influence of the above factors on the transmission of malaria. This work is inspired by the following biological questions: (I) What is the impact of vector-bias on malaria transmission? (II) What effect does the maturity and incubation period of seasonal temperature changes have on the spread of malaria? (III) Does the movement of mosquitoes and humans have different effects for the transmission of malaria in different areas? These issues that are closely related to the spread and control of malaria need to be explored in more depth. How to add the temperature-dependent maturity period to the model is the first issue we consider, and the detailed model derivation is also the difficulty of this paper. The existence of temperature-dependent periodic delays makes it difficult to understand the asymptotic behavior of the model. We consider the stage structure, because the mosquitoes in the aquatic stage cannot fly, which makes it impossible for us to directly use the theoretical results in references [22,24]. Our theoretical results indicate the model has a threshold parameter, the basic reproduction number R0, which determines wether the malaria can persist or become extinct in susceptible populations.

    The structure of the rest of this paper is shown below. Section 2 establishes human and mosquito populations models and malaria transmission model. This section also gives the biological significance of these parameters in models. Section 3 is devoted to well-posedness and some properties of the model. In Section 4, we first define the basic reproduction number R0, and then analyze the threshold dynamics of models in terms of R0. Section 5 gives the basic framework of the system threshold dynamics through numerical simulation, and discusses the influence of parameters on basic reproduction number, mosquitoes and humans in the process of malaria transmission. A brief discussion section concludes the paper.

    Motivated by models in references [19,25,26], we introduce a seasonal stage-structured model of malaria to explain the cross-infection between mosquitoes and individuals. In this subsection, we first consider the mosquito population. According to the biological cycle, mosquitoes can go through several stages and have different habits at various stages. We divide them into immature aquatic stage M(t,x) and mature stage V(t,x) at time t0 and location x in a spatial bounded domain ΩRn with smooth boundary Ω. Let η(t,a,x) denote the density of the mosquito population at time t, chronological age a0, and location x. Metz and Diekmann [27] gave

    {(t+a)η(t,a,x)=D(a)Δη(t,a,x)d(t,a,x)η(t,a,x),η(0,a,x)=ˆϕ(a,x),a0,xˉΩ,η(t,0,x)=Λv(t,x),t0,xˉΩ, (2.1)

    through the standard argument for structured populations and spatial diffusion. Here, Δ denotes the usual Laplacian operator. ˆϕ(,x) represents the non-negative initial age and location distribution of mosquitoes, the term Λv(t,x) is the recruitment rate of the immature stage, which is Hölder continuous and nonnegative function on R×ˉΩ, D(a) is the age-related diffusion coefficient, and d(t,a,x) indicates the death rate of mosquitoes with time t, age a and location x, which is Hölder continuous and positive on R×ˉΩ. Λv(t,) and d(t,,) are periodic in t and have the same period ω.

    Let us first assume some age thresholds so that the growth of mosquitoes can be divided into discrete age stages. When modeling species dynamics with a stage structure, we divide mosquitoes into immature stage (also known as aquatic stage) and mature stage, the time delay is the mature period. The temperature T(t) is an ω-periodic function of time t. Generally speaking, the maturity age is determined by seasonal changes in temperature. We will use a time-dependent positive function τ1(t) to describe the duration from newborn to adult. That is, an individual at time t becomes mature only if its age exceeds τ1(t), in other words, the adult mosquitoes that develop at time t enter the aquatic stage at tτ1(t). We assume that the maximum chronological age of aquatic mosquitoes at time t is A(t). Here and in the following, τ1(t) and A(t) are C1 (a section of continuous differentiable functions) periodic functions of time t with the same periodic ω. The relationship between time-dependent threshold age and time-dependent delay can be indicated as follows:

    A(t)=τ1(t).
    Figure 1.  At time t the mosquitoes reach the threshold age A(t) and develop into the mature stage.

    Here we introduce another quantity f to describe the maturity level of mosquitoes. Let f=fM=0, that is, the maturity level of mosquitoes that have just entered the aquatic stage is 0. At time t, the mosquitoes reach the maximum chronological age A(t) of the aquatic stage. At this time, the maturity level is also the maximum, fV. Then one has

    fVfM=ttτ1(t)ˆρ(r1)dr1

    where ˆρ(r1) is the proportion of maturity at time r1. Taking the derivative of the above formula with respect to t, we can get

    1τ1(t)=ˆρ(t)ˆρ(tτ1(t)),

    and hence 1τ1(t)>0. Recall that τ1(t) represents the time it takes for a mosquito to go from the aquatic stage to the adult stage at time t. For the case where δ>0 is small enough, if the mosquito is an adult mosquito at time t, then it is still an adult mosquito at time t+δ. Thus τ1(t+δ)<τ1(t)+δ, and hence, τ1(t)=limδ0τ1(t+δ)τ1(t)δ1. This condition also makes sense biologically. 1τ1(t)>0, is called maturation rate [26], which excludes the possibility of mosquitoes returning from the mature stage to the immature stage unless they are born. Mathematically speaking, tτ1(t) is strictly increased with respect to t, that is, since the developmental rate depends only on time, it is impossible for juveniles to reach maturity before they are born.

    M(t,x) and V(t,x) denote the density of mosquitoes in the aquatic and mature stages at time t and location x, respectively. The accumulative density between the two age thresholds can be recorded as the population size at each stage. In particular, for the aquatic stage M and the mature stage V (for the sake of simplicity, we use M and V instead of M(t,x) and V(t,x), respectively), we have the following mathematical expressions:

    M=A(t)0η(t,a,x)da,V=A(t)η(t,a,x)da. (2.2)

    Next, we put forward natural biological assumptions about birth rate, death rate and diffusion coefficient in Eq (2.1). Since immature mosquitoes cannot fly, we assume that the diffusion coefficient satisfies

    D(a)={0,0a<A(t),Dv,aA(t).

    At each age stage, all mosquitoes undergo the same birth rate and death rate which are age-independent, then the natural death rates in the aquatic and mature stages are μm(t,x) and μv(t,x), respectively, i.e.,

    d(t,a,x)={μm(t,x),0a<A(t),μv(t,x),aA(t),

    where μm(t,x) and μv(t,x) are non-negative non-trivial and Hölder continuous on R×ˉΩ and ω-periodic in t. Differentiating the equations in system (2.2) with respect to time t on both sides yields

    Mt=A(t)0η(t,a,x)tda+η(t,A(t),x)A(t)=μm(t,x)Mη(t,A(t),x)+η(t,0,x)+η(t,A(t),x)A(t),

    and

    Vt=A(t)η(t,a,x)tdaη(t,A(t),x)A(t)=DvΔVμv(t,x)Vη(t,,x)+(1A(t))η(t,A(t),x).

    In a biological sense, because no mosquito can live forever, η(t,,x)=0. η(t,0,x) represents the inflow rate of mosquitoes in the aquatic stage at time t, then η(t,0,x)=Λv(t,x) is a function of t.

    {Mt=Λv(t,x)μm(t,x)M(1A(t))η(t,A(t),x),Vt=DvΔVμv(t,x)V+(1A(t))η(t,A(t),x).

    Mathematically, assume that the delay τ1(t) is continuously differentiable in [0,), bounded and far away from zero and infinity. The expression of η(t,A(t),x) in the above system will be obtained by integrating along characteristics. For tA(t) and t<A(t), η(t,A(t),x) has different expressions. Without loss of generality, we study tA(t) (see, [28]), which is practicable due to the boundedness of A(t), because we are concerned with long-term behavior of population dynamics. Here, η(tτ1(t),0,)=Λv(tτ1(t),). Let t=t0+h, a=a0+h, P(h,x)=η(t0+h,a0+h,x) and ˆd(h,x)=d(t0+h,a0+h,x), then

    P(h,x)h=η(t0+h,a0+h,x)hdhdt+η(t0+h,a0+h,x)hdhda=(t+a)η(t,a,x)=DvΔP(h,x)ˆd(h,x)P(h,x). (2.3)

    Integrating Eq (2.3) from h1 to h2 yields

    P(h2,x)=P(h1,x)eh2h1ˆd(r2,)dr2.

    For tA(t), setting t0=tA(t), h=A(t) and a0=0, we have

    η(t,A(t),x)=η(tA(t),0,x)ettA(t)μm(r2,x)dr2=Λv(tτ1(t),x)ettτ1(t)μm(r2,x)dr2.

    Since we are concerned about the long-term behavior of population dynamics, the model for tA(t) is as follows

    {Mt=Λv(t,x)μm(t,x)M(1τ1(t))Λv(tτ1(t),x)ettτ1(t)μm(r2,x)dr2,xΩ,Vt=DvΔVμv(t,x)V+(1τ1(t))Λv(tτ1(t),x)ettτ1(t)μm(r2,x)dr2,xΩ,Vν=0,xΩ, (2.4)

    where ν is the outward unit normal vector on Ω and ν indicates the normal derivative along ν on Ω, ettτ1(t)μm(r2,x)dr2 represents probability surviving through natural death during aquatic stage. It is easy to see that the equation V(t,x) is decoupled, so we can get the following system

    {Vt=DvΔVμv(t,x)V+(1τ1(t))Λv(tτ1(t),x)ettτ1(t)μm(r2,x)dr2,xΩ,Vν=0,xΩ. (2.5)

    for any tA(t). Alternatively, according to biological significance, M and V can be written as integral form. Note that τ1(t) is the maturation time of aquatic mosquitoes at time t. Thence, aquatic mosquitoes at time t contain all newly born mosquitoes at a previous time ξ1 with ξ1(tτ1(t),t) and survive to time t. Hence, one has

    M=ttτ1(t)Λv(ξ1,x)etξ1μm(r2,x)dr2dξ1,xˉΩ. (2.6)

    Similarly, for any xˉΩ,

    V=ΩΓ(t,0,x,y)V(0,y)dy+t0ΩΓ(t,ξ1,x,y)Λv(ξ1τ1(ξ1),y)(1τ1(ξ1))eξ1ξ1τ1(ξ1)μm(r2,)dr2dydξ1, (2.7)

    where V(0,) is the initial data of V(t,).

    Figure 2.  Experience of mature mosquitoes at time t, age a.

    The first term in Eq (2.7) represents those mosquitoes that were adults at time 0 at position y and are still in V class at time t at position x by diffusion. The second term in Eq (2.7) denotes the total density of mosquitoes at time t. These mosquitoes were newly developed adults at time ξ1 and survived until time t. These new adults developed from the aquatic stage at previous time ξ1τ1(ξ1), which grown through τ1(ξ1) periodic of time and developed into adult stage at time ξ1 with "maturation rate" 1τ1(ξ1). Γ(t,t0,x,y) with tt00 is the Green function related to Vt=DvΔVμv(t,x)V that meets the Neumann boundary condition.

    Since the incubation period of the pathogen in the human is 1–12 days, which is short compared with the human life span, we do not consider the incubation period of the pathogen in the human body. In this subsection, inspired by [13,19,20,24,29], we give a model of malaria transmission, including the mobility of mosquitoes and individuals, vector-bias, temperature-dependent maturation delay and EIP. Because mosquitoes are the vectors of malaria transmission from person to person, the modeling takes into account the dynamics of human and mosquito populations. One supposes that individuals will become susceptible again after recovery. Thus, humans include two epidemiological types: susceptible (Sh) and infected (Ih) compartments. Assume that all populations walk unbiased and random. Since climate change has almost no impact on human birth and death, according to [29], the following reaction-diffusion equation is suitable:

    Nh(t,x)t=DhΔNh(t,x)+Λh(x)μh(x)Nh(t,x),t>0,xˉΩ. (2.8)

    Dh>0 is the human diffusion coefficient. Suppose the dispersal pattern is random diffusion, that is, individual walkers use a fixed step to walk randomly on the solid line [30,31]. Λh(x) and μh(x) are the recruitment rate and natural death rate of human at location x, respectively. In particular, we should note that the disease dose not have a significant impact on the mobility of individuals and mosquitoes. It is mathematically assumed that the diffusion coefficient Dh>0 of all individuals is the same, and the diffusion coefficient Dv>0 of all mosquitoes is recognized as the same. Assuming that no population flow crosses the boundary Ω, therefore, the Neumann boundary condition is imposed:

    Nh(t,x)ν=0,t>0,xΩ, (2.9)

    here, the definition of ν is same as Section 2.1. Following Theorems 3.1.5 and 3.1.6 in reference [32] and Section 2 in reference [29], we realize that the positive steady state N(x)C(ˉΩ,R+){0} of systems (2.8) and (2.9) is globally attractive, where C(ˉΩ,R+) is Banach space of continuous functions from ˉΩ to R+. According to Appendix 8.1 in reference [33], it is assumed that the total human density of at time t and location x is stable at N(x), i.e., Nh(t,x)N(x), for t0 and xΩ. Then we know that Sh=N(x)Ih. Here and later, for the sake of simplicity, we use Sh and Ih instead of Sh(t,x) and Ih(t,x), respectively.

    From a biological and epidemiological point of view, immature mosquitoes do not take part in infection, and are primarily in a waiting period, because they cannot fly and bite humans. Since adult male mosquitoes do not eat blood meal, they cannot transmit and acquire the virus, so only adult female mosquitoes are modeled. Mature female mosquitoes are divided into susceptible (Sv(t,x)), exposed (Ev(t,x)) and infectious (Ev(t,x)), t is time and x denotes location. The exposed mosquitoes are already infected, but not infectious. From Section 2.1, we know that V=Sv+Ev+Iv is the total density of female adult mosquitoes and satisfies reaction-diffusion system (2.5). Here and later, for the sake of simplicity, we use V, Sv, Ev and Iv instead of V(t,x), Sv(t,x), Ev(t,x) and Iv(t,x), respectively.

    Vector-bias effect is a situation in which mosquitoes display the weight of preference when choosing a host. In order to add a vector-bias term to the model, parameters p and l are introduced, which are the probability of mosquitoes arriving at humans randomly, and select infectious and susceptible humans, respectively [13]. Since mosquitoes are more attracted to infected humans, we assume that the vector-bias parameter χ=pl1. To comprehend the influences of spatial structure and climate change, we allow a space-dependent bite rate β(t,x) of female mosquitoes to be Hölder continuous non-negative non-trivial on R×ˉΩ with β(t,x)0, representing the number of mosquitoes biting per unit time at location x, and β(t,) is a ω-period function of time t. Let B represent the mosquitoes bite incident, and A be the event of mosquitoes reaching humans infected with malaria. Then P(A)=IhN(x), P(Ac)=1IhN(x), P(B|A)=p and P(B|Ac)=l. So the probability that an infectious human is chosen as P(A|B),

    then one has

    P(Ac|B)=l[N(x)Ih]pIh+l[N(x)Ih]andP(A|B)=pIhpIh+l[N(x)Ih],

    which are the probability of infected individuals and susceptible individuals are bitten. Then at time t and location x, the numbers of newly infected individuals and mosquitoes per unit time are

    cβ(t,x)l[N(x)Ih]IvpIh+l[N(x)Ih]andαβ(t,x)pIhSvpIh+l[N(x)Ih],

    respectively, where c is the probability that an infected mosquito will spread to susceptible individuals per bite, and α is probability that an infected individual will transmission to susceptible mosquitoes per bite.

    During the malaria transmission cycle, the EIP of parasite in mosquitoes is one of the most pivotal parameters. Mosquitoes can fly around during EIP. The malaria parasites undergo different growth stages before migrating to the salivary glands where can spread to humans. The degree of this development is related to environmental factors, such as temperature. EIP is exceedingly sensitive to temperature. Therefore, it is necessary to incorporate this seasonal forced incubation period when describing malaria transmission. We will use a time-dependent positive function τ2(t) to describe the duration from new infection to infectiousness. That is, a mosquito is infectious at time t only if its age of infection exceeds τ2(t), where τ2(t) is a C1 periodic function in [0,) with the periodic ω, and tτ2(t) is strictly increasing in t.

    Based on the above discussions, we will pay attention to the dynamics of the following malaria transmission reaction-diffusion model

    {Mt=Λv(t,x)(1τ1(t))Λv(tτ1(t),x)ettτ1(t)μm(r2,x)dr2μm(t,x)M,xΩ,Svt=DvΔSv+(1τ1(t))Λv(tτ1(t),x)ettτ1(t)μm(r2,x)dr2αβ(t,x)pIhSvpIh+l[N(x)Ih]μv(t,x)Sv,xΩ,Evt=DvΔEvEm(t,x)μv(t,x)Ev+αβ(t,x)pIhSvpIh+l[N(x)Ih],xΩ,Ivt=DvΔIv+Em(t,x)μv(t,x)Iv,xΩ,Iht=DhΔIh(μh(x)+ϱ(x))Ih+cβ(t,x)l[N(x)Ih]IvpIh+l[N(x)Ih],xΩ,Svν=Evν=Ivν=Ihν=0,xΩ, (2.10)

    for t>0. Similar to the derivation in Section 2.1, we get

    Em(t,x)=(1τ2(t))ΩΓ(t,tτ2(t),x,y)αβ(tτ2(t),y)pIh(tτ2(t),y)Sv(tτ2(t),y)pIh(tτ2(t),y)+l[N(y)Ih(tτ2(t),y)]dy.

    Substituting Em(t,x) into system (2.10), one has

    {Mt=Λv(t,x)(1τ1(t))Λv(tτ1(t),x)ettτ1(t)μm(r,x)drμm(t,x)M,xΩ,Svt=DvΔSv+(1τ1(t))Λv(tτ1(t),x)ettτ1(t)μm(r,x)drαβ(t,x)pIhSvpIh+l[N(x)Ih]μv(t,x)Sv,xΩ,Evt=DvΔEv+αβ(t,x)pIhSvpIh+l[N(x)Ih]μv(t,x)Ev(1τ2(t))ΩΓ(t,tτ2(t),x,y)αβ(tτ2(t),y)pIh(tτ2(t),y)Sv(tτ2(t),y)pIh(tτ2(t),y)+l[N(y)Ih(tτ2(t),y)]dy,xΩ,Ivt=DvΔIvμv(t,x)Iv+(1τ2(t))ΩΓ(t,tτ2(t),x,y)αβ(tτ2(t),y)pIh(tτ2(t),y)Sv(tτ2(t),y)pIh(tτ2(t),y)+l[N(y)Ih(tτ2(t),y)]dy,xΩ,Iht=DhΔIh(μh(x)+ϱ(x))Ih+cβ(t,x)l[N(x)Ih]IvpIh+l[N(x)Ih],xΩ,Svν=Evν=Ivν=Ihν=0,xΩ, (2.11)

    for t>0, where α, p, l and c are positive constants, and ϱ(x) denotes the recovery rate of humans in position x, which is Hölder continuous and positive on Ω. Neumann boundary condition indicates that all individuals and mosquitoes still remain in this domain Ω.

    Since M(t,x) and Ev(t,x) of the system (2.11) are decoupled from the other equations, it is sufficient to explore the system

    {u1t=DvΔu1μv(t,x)u1+(1τ1(t))Λv(tτ1(t),x)ettτ1(t)μm(r,x)drαβ(t,x)pu3u1pu3+l[N(x)u3],xΩ,u2t=DvΔu2μv(t,x)u2+(1τ2(t))ΩΓ(t,tτ2(t),x,y)αβ(tτ2(t),y)pu3(tτ2(t),y)u1(tτ2(t),y)pu3(tτ2(t),y)+l[N(y)u3(tτ2(t),y)]dy,xΩ,u3t=DhΔu3+cβ(t,x)l[N(x)u3]u2pu3+l[N(x)u3](μh(x)+ϱ(x))u3,xΩ,u1t=u2t=u3t=0,xΩ, (2.12)

    for t>0, where (u1,u2,u3)=(Sv,Iv,Ih).

    In this section, the well-posedness of the model is proved. Let X1:=C(ˉΩ,R3) be the Banach space of continuous functions from ˉΩ to R3 the supremum norm be , and \mathbb{X}_{1}^{+}: = C\left(\bar{\Omega}, \mathbb{R}_{+}^{3}\right) . Let \hat{\tau}_{1}: = \max\limits _{t \in[0, \omega]}\tau_{1}(t) , \hat{\tau}_{2}: = \max\limits _{t \in[0, \omega]}\tau_{2}(t) , and \hat{\tau}: = \max\left\{\hat{\tau}_{1}, \hat{\tau}_{2}\right\} . Define \mathcal{X}_{1}: = C([-\hat{\tau}, 0], \mathbb{X}_{1}) to be a Banach space, with the norm \|\phi\| = \max\limits_{\theta\in[-\hat{\tau}, 0]}\|\phi(\theta)\|_{\mathbb{X}_{1}} , \forall \phi\in \mathcal{X}_{1} , and \mathcal{X}_{1}^{+}: = C([-\hat{\tau}, 0], \mathbb{X}_{1}^{+}) . Then (\mathbb{X}_{1}, \mathbb{X}_{1}^{+}) and (\mathcal{X}_{1}, \mathcal{X}_{1}^{+}) are order Banach spaces. A function z:[-\hat{\tau}, \sigma)\rightarrow \mathbb{X}_{1} for \sigma > 0 is given and z_{t}\in \mathcal{X}_{1} is defined by

    \begin{equation} \nonumber \begin{aligned} z_{t}(\theta) = (z_{1}(t+\theta),z_{2}(t+\theta),z_{3}(t+\theta)), \; \; \; \forall \theta\in[-\hat{\tau},0], \end{aligned} \end{equation}

    for any t\in[0, \sigma) .

    Denote \mathbb{Y}: = C\left(\bar{\Omega}, \mathbb{R}\right) , \mathbb{Y}^{+}: = C\left(\bar{\Omega}, \mathbb{R}_{+}\right) , \mathbb{Y}_{h}: = \{\varphi \in \mathbb{Y}^{+}:0\leq \varphi(x)\leq N(x), \forall x\in \bar{\Omega}\} . Suppose that T_{1}(t, s), T_{2}(t, s): \mathbb{Y}\rightarrow \mathbb{Y} , are the evolution operators associated

    \begin{equation} \nonumber \begin{aligned} \frac{\partial u_{1}}{\partial t} = D_{v}\Delta u_{1}-\mu_{v}(t,x)u_{1}(t,x): = A_{1}(t)u_{1}, \end{aligned} \end{equation}

    and

    \begin{equation} \nonumber \begin{aligned} \frac{\partial u_{3}}{\partial t}& = D_{h}\Delta u_{3}-(\mu_{h}(x)+\varrho(x))u_{3} : = A_{2}u_{3}, \end{aligned} \end{equation}

    obey the Neumann boundary condition, respectively. Since T_{i}(t, s) = T_{i}(t-s) (i = 1, 2) and \mu_{v}(t, \cdot) is \omega -periodic in t , = [34,Lemma 6.1] indicates that T_{i}(t+\omega, s+\omega) = T_{i}(t, s) for (t, s)\in \mathbb{R}^{2} with t\geq s . Furthermore, T_{i}(t, s) is compact and strongly positive [35,Section 7.1 and Corollary 7.2.3] for (t, s)\in \mathbb{R}^{2} with t > s (i = 1, 2) . Let T(t, s) = {\rm{diag}}\{T_{1}(t, s), T_{1}(t, s), T_{2}(t, s)\}: \mathbb{X}_{1}\rightarrow \mathbb{X}_{1} , t\geq s , is a strongly continuous semigroup and A(t) = {\rm{diag}}\{A_{1}(t), A_{1}(t), A_{2}\} , defined on D(A(t)) = D(A_{1}(t))\times D(A_{1}(t))\times D(A_{2}) . A_{1}(t) is defined by

    \begin{equation} \nonumber \begin{aligned} D(A_{1}(t)):& = \left\{\varphi\in C^{2}\left(\bar{\Omega}\right):\frac{\partial \varphi}{\partial \nu} = 0\; {\rm{on}} \; \partial\Omega\right\},\\ A_{1}(t)\varphi& = D_{v}\Delta \varphi-\mu_{v}(t,x)\varphi,\; \; \; \; \forall \varphi\in D(A_{1}(t)), \end{aligned} \end{equation}

    and A_{2} is defined by

    \begin{equation} \nonumber \begin{aligned} D(A_{2}):& = \left\{\varphi\in C^{2}\left(\bar{\Omega}\right):\frac{\partial \varphi}{\partial \nu} = 0\; {\rm{on}} \; \partial\Omega\right\},\\ A_{2}\varphi& = D_{h}\Delta \varphi-(\mu_{h}(x)+\varrho(x))\varphi,\; \; \; \; \forall \varphi\in D(A_{2}). \end{aligned} \end{equation}

    Set

    \begin{equation} \nonumber \begin{aligned} W_{h} = C([-\hat{\tau},0],\mathbb{Y}^{+}\times \mathbb{Y}^{+}\times \mathbb{Y}_{h}). \end{aligned} \end{equation}

    Write F = (F_{1}, F_{2}, F_{3}):[0, \infty)\times W_{h}\rightarrow \mathbb{X}_{1} as

    \begin{equation} \left\{ \begin{aligned} F_{1}(t,\phi) = &(1-\tau_{1}^{\prime}(t))\Lambda_{v}(t-\tau_{1}(t),\cdot)e^{-\int_{t-\tau_{1}(t)}^{t}\mu_{m}(r_{2},\cdot) {\rm{d}}r_{2}} -\frac{\alpha\beta(t,\cdot)p\phi_{3}(0,\cdot)\phi_{1}(0,\cdot)}{p\phi_{3}(0,\cdot) +l[N(\cdot)-\phi_{3}(0,\cdot)]},\\ F_{2}(t,\phi) = &(1-\tau^{\prime}_{2}(t)) \int_{\Omega} \Gamma(t,t-\tau_{2}(t),\cdot,y)\frac{\alpha\beta(t-\tau_{2}(t),y)p\phi_{3}(-\tau_{2}(t),y)\phi_{1}(-\tau_{2}(t),y)}{p\phi_{3}(-\tau_{2}(t),y) +l[N(y)-\phi_{3}(-\tau_{2}(t),y)]} {\rm{d}}y,\\ F_{3}(t,\phi) = &\frac{c\beta(t,\cdot)l[N(\cdot)-\phi_{3}(0,\cdot)]\phi_{2}(0,\cdot)}{p\phi_{3}(0,\cdot)+l[N(\cdot)-\phi_{3}(0,\cdot)]}, \end{aligned} \right. \end{equation} (3.1)

    for t\geq 0 , x\in\bar{\Omega} and \phi = (\phi_{1}, \phi_{2}, \phi_{3})\in W_{h} . Rewrite system (2.12) as

    \begin{equation} \left\{ \begin{aligned} \frac{d u}{\partial t}& = A(t)u+F(t,u_{t}), &t > 0,\\ u_{0}& = \phi, \end{aligned} \right. \end{equation} (3.2)

    it then follows from references [36,Corollary 4 and Theorem 1] and [37,Corollary 8.1.3] that for any \phi \in W_{h} , a mild solution can be get as a continuous solution of the integral equation:

    \begin{equation} \left\{ \begin{aligned} u(t,\phi)& = T(t,0)\phi(0)+\int_{0}^{t}T(t,s)F(s,u_{s}) {\rm{d}}s, \; \; \; \forall t > 0,\\ u_{0}& = \phi. \end{aligned} \right. \end{equation} (3.3)

    Lemma 3.1. For any \phi\in W_{h} , system (2.12) has the unique solution, denoted by z(t, \cdot, \phi) on its maximal existence interval [0, t_{\phi}) with z_{0} = \phi , where t_{\phi}\leq \infty .Moreover, z(t, \cdot, \phi)\in \mathbb{Y}^{+}\times\mathbb{Y}^{+}\times\mathbb{Y}_{h} for any t\in [0, t_{\phi}) , andthis solution is a classical solution of system (2.12) for all t > \hat{\tau} .

    Before studying the global existence, go back to system (2.11) to make more information. We impose the compatibility conditions based on the biological meaning of \tau_{1}(t) and \tau_{2}(t) :

    \begin{equation} \begin{aligned} M(0,\cdot) = \int_{-\tau_{1}(0)}^{0}\Lambda_{v}(\xi_{1},\cdot)e^{-\int_{\xi_{1}}^{0}\mu_{m}(r_{2},\cdot) {\rm{d}}r_{2}} {\rm{d}}\xi_{1}, \end{aligned} \end{equation} (3.4)

    and

    \begin{equation} \begin{aligned} E_{v}(0,\cdot) = \int_{-\tau_{2}(0)}^{0}T_{1}(0,\xi_{1})\frac{\alpha\beta(\xi_{1},\cdot)pI_{h}(\xi_{1},\cdot)S_{v}(\xi_{1},\cdot)}{pI_{h}(\xi_{1},\cdot) +l[N(\cdot)-I_{h}(\xi_{1},\cdot)]} {\rm{d}}\xi_{1}. \end{aligned} \end{equation} (3.5)

    Now we define

    \begin{equation} \nonumber \begin{aligned} \mathcal{D}: = &\bigg\{\psi\in C\left([-\hat{\tau},0],C\left(\bar{\Omega},\mathbb{R}^{5}_{+}\right)\right):\psi_{5}(s,\cdot)\leq N(\cdot),\; \; \; \forall s\in [-\hat{\tau},0],\\ \; \; &\psi_{1}(0,\cdot) = \int_{-\tau_{1}(0)}^{0}\Lambda_{v}(\xi_{1},\cdot)e^{-\int_{\xi_{1}}^{0}\mu_{m}(r_{2},\cdot) {\rm{d}}r_{2}} {\rm{d}}\xi_{1},\\ \; \; &\psi_{3}(0,\cdot) = \int_{-\tau_{2}(0)}^{0}T_{1}(0,\xi_{1})\frac{\alpha\beta(\xi_{1},\cdot)pI_{h}(\xi_{1},\cdot)S_{v}(\xi_{1},\cdot)}{pI_{h}(\xi_{1},\cdot) +l[N(\cdot)-I_{h}(\xi_{1},\cdot)]} {\rm{d}}\xi_{1}\bigg\}. \end{aligned} \end{equation}

    Next, for any \psi\in \mathcal{D} , system (2.11) has the unique solution U(t, \cdot) = (M(t, \cdot), S_{v}(t, \cdot), E_{v}(t, \cdot), I_{v}(t, \cdot), I_{h}(t, \cdot)) satisfying U_{0} = \psi . Corollary 4 in reference [36] implies that S_{v}(t, x)\geq 0 , I_{v}(t, x)\geq 0 and I_{h}(t, x)\geq 0 on maximal existence interval. According to the uniqueness of the solutions and the compatibility conditions (3.4) and (3.5), one has M(t, \cdot) (see(2.6)) and

    \begin{equation} \begin{aligned} E_{v}(t,\cdot) = \int_{t-\tau_{2}(t)}^{t}T_{1}(t,\xi_{1})\frac{\alpha\beta(\xi_{1},\cdot)pI_{h}(\xi_{1},\cdot)S_{v}(\xi_{1},\cdot)}{pI_{h}(\xi_{1},\cdot) +l[N(\cdot)-I_{h}(\xi_{1},\cdot)]} {\rm{d}}\xi_{1}. \end{aligned} \end{equation} (3.6)

    Hence, M(t, x)\geq 0 and E_{v}(t, x)\geq 0 on the maximal interval of existence.

    Lemma 3.2. For any \phi\in W_{h} , system (2.12)has a unique solution u(t, \cdot, \phi) on [0, \infty) with u_{0} = \phi .Moreover, system (2.12) generates an \omega -periodic semiflow Q_{t}:W_{h}\rightarrow W_{h} , defined by Q_{t}: = u_{t}(\cdot) , \forall t\geq 0 , i.e., Q_{t}(\phi)(s, x) = u(t+s, x, \phi) , \forall\phi\in W_{h} , t\geq 0 , s\in [-\hat{\tau}, 0] , x\in \Omega , and Q: = Q_{\omega} has a strong global attractor in W_{h} .

    The proof of Lemma 3.2 is given in Appendix A.

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

    The current main work is to explore the extinction and persistence of malaria in a spatially heterogeneous environment. The basic reproduction number R_{0} is one of the most necessary concepts in epidemiology. It is usually defined as the mean number of secondary infections during the entire period when a typical infected person is introduced into a completely susceptible population. From a biological point of view, the basic reproduction number is a characteristic of the risk of infection. More meaningfully, it describes the threshold behavior of many infectious disease models.

    Let u_{2} = u_{3} = 0 , we have the equation of mosquitoes compartment

    \begin{equation} \begin{aligned} \frac{\partial u_{1}}{\partial t} = &D_{v}\Delta u_{1}+(1-\tau^{\prime}_{1}(t))\Lambda_{v}(t-\tau_{1}(t),x)e^{-\int_{t-\tau_{1}(t)}^{t}\mu_{m}(r_{2},x) {\rm{d}}r_{2}}-\mu_{v}(t,x)u_{1}, \end{aligned} \end{equation} (4.1)

    obeys the Neumann boundary condition. Based on reference [38,Lemma 2.1], it is easy to know that system (4.1) has a unique globally attractive positive \omega -periodic solution u_{1}^{*}(t, \cdot) in \mathbb{Y}^{+}\backslash \{0\} . Linearize system (2.12) at (u_{1}^{*}(t, \cdot), 0, 0) and consider the equations of infectious compartments,

    \begin{equation} \left\{ \begin{aligned} \frac{\partial w_{1}}{\partial t} = &D_{v}\Delta w_{1}-\mu_{v}(t,x)w_{1}+(1-\tau^{\prime}_{2}(t))\int_{\Omega} \Gamma(t,t-\tau_{2}(t),x,y)\cdot\\ &\frac{\alpha\beta(t-\tau_{2}(t),y)pw_{2}(t-\tau_{2}(t),y)u_{1}^{*}(t-\tau_{2}(t),y)}{lN(y)} {\rm{d}}y,&x\in \Omega,\\ \frac{\partial w_{2}}{\partial t} = &D_{h}\Delta w_{2}+c\beta(t,x)w_{1}-(\mu_{h}(x)+\varrho(x))w_{2}, &x\in \Omega,\\ \frac{\partial w_{1}}{\partial \nu} = &\frac{\partial w_{2}}{\partial \nu} = 0, &x\in \partial\Omega, \end{aligned} \right. \end{equation} (4.2)

    for t > 0 . Denote w: = (w_{1}, w_{2}) = (u_{2}, u_{3}) .

    Set \mathbb{X}_{2}: = C\left(\bar{\Omega}, \mathbb{R}^{2}\right) , \mathbb{X}_{2}^{+}: = C\left(\bar{\Omega}, \mathbb{R}_{+}^{2}\right) , and C_{\omega}(\mathbb{R}, \mathbb{X}_{2}) to be the Banach space consisting of all \omega -periodic and continuous functions from \mathbb{R} to \mathbb{X}_{2} , where \|\psi\|_{C_{\omega}(\mathbb{R}, \mathbb{X}_{2})}: = \max\limits_{\theta\in[0, \omega]}\|\psi(\theta)\|_{\mathbb{X}_{2}} for any \psi\in C_{\omega}(\mathbb{R}, \mathbb{X}_{2}) . Let \mathcal{X}_{2}: = C([-\hat{\tau}_{2}, 0], \mathbb{X}_{2}) and \mathcal{X}_{2}^{+}: = C([-\hat{\tau}_{2}, 0], \mathbb{X}_{2}^{+}) . Write \mathcal{F}(t):\mathcal{X}_{2} \rightarrow \mathbb{X}_{2} as

    \begin{equation} \nonumber \begin{aligned} &\mathcal{F}(t)\left(\begin{array}{c} \tilde{\psi}_{1}\\ \tilde{\psi}_{2} \end{array}\right) = \left(\begin{array}{c} (1-\tau^{\prime}_{2}(t))\int_{\Omega} \Gamma(t,t-\tau_{2}(t),\cdot,y) \frac{\alpha\beta(t-\tau_{2}(t),y)p\tilde{\psi}_{2}(-\tau_{2}(t),y)u_{1}^{*}(t-\tau_{2}(t),y) }{lN(y)} {\rm{d}}y \\ c\beta(t,\cdot)\tilde{\psi}_{1}(0,\cdot) \end{array}\right) \end{aligned} \end{equation}

    for any t\in \mathbb{R} , \tilde{\psi} = \left(\tilde{\psi}_{1}, \tilde{\psi}_{2}\right)\in \mathcal{X}_{2} and -V_{1}(t)w = \bar{D}\Delta w-W(t)w , where \bar{D} = {\rm{diag}}(D_{v}, D_{h}) and

    -[W(t)](x) = \left(\begin{array}{cc} -\mu_{v}(t,x)&0\\ 0 & -\left(\mu_{h}(x)+\varrho(x)\right) \end{array}\right), \; \; \; x \in \bar{\Omega}.

    Denote that \Phi(t, s) = {\rm{diag}}(T_{1}(t, s), T_{2}(t, s)) is an evolution operator related to the following system

    \begin{equation} \begin{aligned} \frac{ {\rm{d}}w}{ {\rm{d}}t} = -V_{1}(t)w, \end{aligned} \end{equation} (4.3)

    t\geq s , with the Neumann boundary condition. Obviously, \Phi(t, s) is a positive operator, that is, \Phi(t, s)\mathbb{X}_{2}^{+}\subseteq \mathbb{X}_{2}^{+} for all t\geq s . Then -V_{1}(t) is resolvent positive, based on reference [39,Theorem 3.12].

    The exponential growth bound of \Phi(t, s) can be written as

    \begin{equation} \nonumber \begin{aligned} \bar{\omega} = \inf\{\tilde{\omega}:\exists \mathcal{L}\geq 1\; {\rm{suth}}\; {\rm{that}}\; \|\Phi(t+s,s)\leq \mathcal{L}e^{\tilde{\omega}t}\|,\; \; \forall s\in \mathbb{R}, t\geq 0\}. \end{aligned} \end{equation}

    By [39,Proposition A.2], one has

    \begin{equation} \nonumber \begin{aligned} \bar{\omega}(\Phi) = \frac{\ln r(\Phi(\omega+\bar{s},\bar{s}))}{\omega}, \quad \bar{s}\in [0,\omega]. \end{aligned} \end{equation}

    According to Krein-Rutman theorem and [40,Lemma 14.2], one has

    \begin{equation} \nonumber \begin{aligned} 0 < r(\Phi(\omega,0)) = \max\{r(T_{1}(\omega,0)),r(T_{2}(\omega,0))\} < 1, \end{aligned} \end{equation}

    here, r(\Phi(\omega, 0)) denotes the spectral radius of \Phi(\omega, 0) . By [39,Proposition 5.6] with s = 0 , we obtain \bar{\omega}(\Phi) < 0 . Obviously, \mathcal{F}(t) and W(t) meet the assumptions: (1) \mathcal{F}(t): \mathcal{X}_{2} \rightarrow \mathbb{X}_{2} is positive, i.e., \mathcal{F}(t)\mathcal{X}_{2}^{+}\subseteq\mathbb{X}_{2}^{+} , for any t\geq 0 . (2) -W(t) is cooperative, and \bar{\omega}(\Phi) < 0 .

    To introduce R_{0} of system (2.12) based on references [41] and [42], it is assumed that both human and mosquito populations are close to (u_{1}^{*}(t, \cdot), 0, 0) . Suppose w\in C_{\omega}(\mathbb{R}, \mathbb{X}_{2}) , and w(t) = w(t) is the initial distribution of infected mosquitoes and humans at t\in\mathbb{R} . For any fixed s\geq 0 , the distribution of newly infected mosquitoes and individuals at time t-s (t\geq s) is \mathcal{F}(t-s)w_{t-s} , which is engendered by the infectious mosquitoes and humans introduced in [t-s-\hat{\tau}_{2}, t-s] . Thus, the distribution of infectious who were newly infected at time t-s and keep survive at time t , is \Phi(t, t-s)\mathcal{F}(t-s)w_{t-s} , t\geq s . Then, at time t , the distribution of the accumulative infectious humans and mosquitoes produced by all people introduced to t at all previous times is

    \begin{equation} \nonumber \begin{aligned} \int_{0}^{+\infty}\Phi(t,t-s)\mathcal{F}(t-s)w_{t-s} {\rm{d}}s = \int_{0}^{+\infty}\Phi(t,t-s)\mathcal{F}(t-s)w(t-s+\cdot) {\rm{d}}s. \end{aligned} \end{equation}

    Define the next generation operator L associated with system (2.12) as

    \begin{equation} \nonumber [Lw](t): = \int_{0}^{+\infty}\Phi(t,t-s)\mathcal{F}(t-s)w(t-s+\cdot) {\rm{d}}s, \; \; \; \; \forall t\in \mathbb{R}, \; \; \; \; w\in C_{\omega}(\mathbb{R}, \mathbb{X}_{2}). \end{equation}

    Inspired by the definition of next generation operators in references [39,42], the basic reproduction number of system (2.12) is defined by the spectral radius L , i.e.,

    \begin{equation} \begin{aligned} R_{0}: = r(L). \end{aligned} \end{equation} (4.4)

    For any \tilde{\psi} \in \mathcal{X}_{2} , set \bar{Q}_{t} to be solution map of system (4.2) on \mathcal{X}_{2} , i.e., \bar{Q}_{t}(\tilde{\psi}) = w_{t}(\tilde{\psi}) , t\geq 0 , where w_{t}(\tilde{\psi})(\theta)(x) = w(t+\theta, x, \tilde{\psi}) = (w_{1}(t+\theta, x, \tilde{\psi}), w_{2}(t+\theta, x, \tilde{\psi})) , and w(t, x, \tilde{\psi}) is the unique solution of system (4.2) with w(\theta, x) = \tilde{\psi}(\theta, x) for all \theta\in[-\hat{\tau}_{2}, 0] , x\in \bar{\Omega} . Thus \bar{Q}: = \bar{Q}_{\omega} is the Poincaré map related to system (4.2). r(\bar{Q}) represents the spectral radius of \bar{Q} . Based on reference [41,Theorem 3.7], the following result can be obtained.

    Lemma 4.1. R_{0}-1 has the same sign as r(\bar{Q})-1 .

    To explore the global dynamic behaviors of system (2.12) according to R_{0} , we demonstrate that system (4.2) engenders an eventually strong monotone periodic semiflow on the phase space:

    \begin{equation} \nonumber \begin{aligned} \mathcal{E}_{1}: = \mathbb{Y} \times C([-\tau_{2}(0),0],\mathbb{Y}), \; {\rm{and}}\; \mathcal{E}_{1}^{+}: = \mathbb{Y}\times C([-\tau_{2}(0),0],\mathbb{Y}^{+}). \end{aligned} \end{equation}

    Thus (\mathcal{E}_{1}, \mathcal{E}_{1}^{+}) is an ordered Banach space. For every \hat{\psi}\in \mathcal{E}_{1}^{+} , \tilde{w}(t, x, \hat{\psi}) = (\tilde{w}_{1}(t, x, \hat{\psi}), \tilde{w}_{2}(t, x, \hat{\psi})) denotes the unique solution of (4.2) with \tilde{w}_{0}(\hat{\psi})(\theta, x) for x\in \bar{\Omega} , \theta\in [-\tau_{2}(0), 0] , where

    \begin{equation} \nonumber \begin{aligned} \tilde{w}_{t}(\hat{\psi})(\theta,x) = \tilde{w}(t+\theta,x,\hat{\psi}) = (\tilde{w}_{1}(t,x),\tilde{w}_{2}(t+\theta_{2},x)), \forall \theta\in [-\tau_{2}(0),0], \forall t\geq 0. \end{aligned} \end{equation}

    Lemma 4.2. For any \hat{\psi}\in \mathcal{E}_{1}^{+} , system (4.2)has a unique nonnegative solution \tilde{w}(t, \cdot, \hat{\psi}) on [0, \infty) with \tilde{w}_{0} = \hat{\psi} .

    The proof of Lemma 4.2 is given in Appendix B.

    Remark 1. Through the uniqueness of solutions in Lemmas 3.2 and 4.2, then for any \tilde{\psi}\in \mathcal{X}_{2}^{+} and \hat{\psi} \in \mathcal{E}_{1}^{+} with \tilde{\psi}_{1}(\cdot) = \hat{\psi}_{1}(\cdot) and \tilde{\psi}_{2}(\theta, \cdot) = \hat{\psi}_{2}(\theta, \cdot) , \forall \theta\in[-\tau_{2}(0), 0] , one has w(t, \cdot, \tilde{\psi}) = \tilde{w}(t, \cdot, \hat{\psi}) , t\geq 0 , where w(t, \cdot, \tilde{\psi}) and \tilde{w}(t, \cdot, \hat{\psi}) are solutions of system (4.2) meeting w_{0} = \tilde{\psi} and \tilde{w}_{0} = \hat{\psi} , respectively.

    For every fixed t\geq 0 , let \tilde{Q}_{t} be the solution map of system (4.2) on the space \mathcal{E}_{1} , i.e., \tilde{Q}_{t}(\hat{\psi}) = \tilde{w}_{t}(\hat{\psi}) , for t\geq 0 and \forall \hat{\psi}\in\mathcal{E}_{1} . Thus, \tilde{Q}: = \tilde{Q}_{\omega} is the Poincaré map related to system (4.2) and r(\tilde{Q}) is the spectral radius of \tilde{Q} . The next lemma indicates that the periodic semiflow \tilde{Q}_{t} is eventually strongly positive.

    Lemma 4.3. For each \xi in \mathcal{E}_{1}^{+} , with \xi\not\equiv 0 , the solution \tilde{w}(t, \cdot, \xi) of (4.2) with \tilde{w}_{0} = \xi , meets \tilde{w}_{i}(t, \cdot) > 0 (i = 1, 2) for all t > 2\hat{\tau} , therefore, \tilde{Q}_{t}\xi \gg 0 for every t > 3\hat{\tau} .

    The proof of Lemma 4.3 is given in Appendix C.

    Fix an integer n_{0} satisfying n_{0}\omega > 3\hat{\tau} . Obviously, \tilde{Q}^{n_{0}} = \tilde{Q}_{n_{0}\omega} is strongly monotone, according to proof of Lemma 4.3. In addition, through similar arguments in reference [43,Lemma 2.6], one can get that \tilde{Q}^{n_{0}} is compact. Apply the Krein-Rutman theorem to the linear operator \tilde{Q}^{n_{0}} , and because r(\tilde{Q}^{n_{0}}) = (r(\tilde{Q}))^{n_{0}} , one has \bar{\lambda} = r(\tilde{Q}) > 0 , among them \bar{\lambda} is the simple eigenvalue of \tilde{Q} and there is a eigenvector \vartheta\in Int(\mathcal{E}_{1}^{+}) that is strongly positive. Based on reference [22,Lemma 3.8], one has the following lemma.

    Lemma 4.4. The spectral radii of the two Poincaré maps \bar{Q}:\mathcal{X}_{2}\rightarrow \mathcal{X}_{2} and \tilde{Q}:\mathcal{E}_{1}\rightarrow \mathcal{E}_{1} are the samei.e., r(\bar{Q}) = r(\tilde{Q}) .Furthermore, R_{0}-1 has the same sign as r(\tilde{Q})-1 .

    The fact that system has a special solution is given by the following lemma, which can be used to study long-term dynamics. The following argument is inspired by references [29,Lemma 5] and [44,Proposition 1.1].

    Lemma 4.5. There is a positive \omega -periodic function w^{*}(t, x) , such that e^{\tilde{\mu}t}w^{*}(t, x) is a solution of (4.2), where \tilde{\mu} = \frac{\ln r(\tilde{Q})}{\omega} .

    The proof of Lemma 4.5 is given in Appendix D.

    Next, we use R_{0} to give the threshold results of the global dynamics of system (2.12).

    Denote

    \begin{equation} \begin{aligned} \mathcal{E}_{2}: = C([-\tau_{2}(0),0],\mathbb{Y}^{+}) \times \mathbb{Y}^{+}\times C([-\tau_{2}(0),0],\mathbb{Y}_{h}). \end{aligned} \end{equation} (4.5)

    Through similar arguments in Lemma 4.2, it can be concluded that for any \phi\in W_{h} and \bar{\varphi}\in \mathcal{E}_{2} with \phi_{1}(\theta, \cdot) = \bar{\varphi}_{1}(\theta, \cdot) , \phi_{2}(0, \cdot) = \bar{\varphi}_{2}(\cdot) , and \phi_{3}(\theta, \cdot) = \bar{\varphi}_{3}(\theta, \cdot) , \forall\theta\in [-\tau_{2}(0), 0] , z(t, \cdot, \phi) = \bar{z}(t, \cdot, \bar{\varphi}) is established, t\geq 0 , where z(t, \cdot, \phi) and \bar{z}(t, \cdot, \bar{\varphi}) are solutions of system (2.12) meeting z_{0} = \phi and \bar{z}_{0} = \bar{\varphi} , respectively. Then we get that solutions of system (2.12) on \mathcal{E}_{2} exist globally on [0, +\infty) and ultimately bounded. In addition, it is easy to get the next result based on the discussion on references [22,Lemma 3.5], [43] and [45,Theorem 2.9].

    Lemma 4.6. Set \hat{Q}_{t} to be the solution map of system (2.12) on the space \mathcal{E}_{2} , i.e., \hat{Q}_{t}(\bar{\varphi}) = u_{t}(\bar{\varphi}) , for t\geq 0 and \forall \bar{\varphi}\in\mathcal{E}_{2} . Thus, \hat{Q}_{t} is an \omega -periodic semiflow on \mathcal{E}_{2} , in this sense, one has

    (i) \hat{Q}_{0} = I ;

    (ii) \hat{Q}_{t+\omega} = \hat{Q}_{t}\circ\hat{Q}_{\omega} , \forall t\geq 0 ;

    (iii) \hat{Q}_{t}(\bar{\varphi}) is continuous in (t, \bar{\varphi})\in[0, \infty)\times \mathcal{E}_{2} .

    Furthermore, \hat{Q}: = \hat{Q}_{\omega} admits a strong global attractor in \mathcal{E}_{2} .

    Next, we prove that the solution of system (2.12) is strictly positive.

    Lemma 4.7. Use u(t, \cdot, \bar{\varphi}) = (u_{1}(t, \cdot, \bar{\varphi}), u_{2}(t, \cdot, \bar{\varphi}), u_{3}(t, \cdot, \bar{\varphi})) to be the solution of (2.12) with the initial data u_{0} = \bar{\varphi}\in \mathcal{E}_{2} , if there is t_{1}\geq 0 such that u_{2}(t_{1}, \cdot, \bar{\varphi})\not\equiv 0 , and u_{3}(t_{1}, \cdot, \bar{\varphi})\not\equiv 0 , in that way, the solution of system (2.12) meets

    \begin{equation} \begin{aligned} u_{2}(t,x,\bar{\varphi}) > 0,\quad u_{3}(t,x,\bar{\varphi}) > 0,\quad \; \forall t > t_{1},\; \; x\in\bar{\Omega}. \end{aligned} \end{equation} (4.6)

    Besides, for \forall t > 0 , x\in \bar{\Omega} and any initial data \bar{\varphi}\in \mathcal{E}_{2} , one has u_{1}(t, x, \bar{\varphi}) > 0 , and

    \begin{equation} \begin{aligned} \liminf\limits_{t\rightarrow \infty} u_{1}(t,x,\bar{\varphi})\geq \varepsilon, \; \; unifromly\; for\; x\in \bar{\Omega}, \end{aligned} \end{equation} (4.7)

    where \varepsilon > 0 is a \bar{\varphi} -independent positive constant.

    Proof. According to the comparison principle for cooperative systems, we have u_{2}(t, x)\geq 0 and u_{3}(t, x)\geq 0 for t > 0 and x\in \bar{\Omega} . Moreover, for a fixed \bar{\varphi}\in \mathcal{E}_{2} , it is easy to know that u_{2}(t, x, \bar{\varphi}) and u_{3}(t, x, \bar{\varphi}) satisfy

    \begin{equation} \nonumber \left\{ \begin{aligned} &\frac{\partial u_{2}}{\partial t}\geq D_{v}\Delta u_{2}-\hat{\mu}_{v}u_{2}, \\ &\frac{\partial u_{3}}{\partial t}\geq D_{h}\Delta u_{3}-(\hat{\mu}_{h}+\hat{\varrho})u_{3}, \\ &\frac{\partial u_{2}}{\partial \nu} = \frac{\partial u_{3}}{\partial \nu} = 0, \quad x\in \partial\Omega, \end{aligned} \right. \end{equation}

    where \hat{\mu}_{v} = \max\limits_{t\in[0, \omega], x\in\bar{\Omega}}\mu(t, x) , \hat{\mu}_{h} = \max\limits_{x\in\bar{\Omega}}\mu_{h}(x) and \hat{\varrho} = \max\limits_{x\in\bar{\Omega}}\varrho(x) . If there exists t_{1}\geq 0 , such that u_{2}(t_{1}, \cdot, \bar{\varphi})\neq 0 and u_{3}(t_{1}, \cdot, \bar{\varphi})\neq 0 , one has u_{2}(t, \cdot, \bar{\varphi}) > 0 and u_{3}(t, \cdot, \bar{\varphi}) > 0 , \forall t > t_{1} based on the strong maximum principle [40,Proposition 13.1]. For t > 0 , let \bar{n}(t, x, \bar{\varphi}_{2}) be the solution of

    \begin{equation} \left\{ \begin{aligned} \frac{\partial \bar{n}(t,x)}{\partial t} = & D_{v}\Delta \bar{n}(t,x)+(1-\tau^{\prime}_{1}(t))\Lambda_{v}(t-\tau_{1}(t),x)e^{-\int^{t}_{t-\tau_{1}(t)}\mu_{m}(r,x) {\rm{d}}r} -\left(\frac{\alpha\beta(t,x)p}{l}+\mu_{v}(t,x)\right) \bar{n}(t,x), &x\in \Omega, \\ \frac{\partial \bar{n}(t,x)}{\partial \nu} = &0, &x\in \partial\Omega. \end{aligned} \right. \end{equation} (4.8)

    Clearly, \Lambda_{v}(t, x)\not \equiv 0 is non-negative non-trivial and Hölder continuous on \mathbb{R}\times \bar{\Omega} . Using the comparison principle, one has

    \begin{equation} \nonumber \begin{aligned} u_{1}(t,x)\geq \bar{n}(t,x) > 0, \; t > 0, \; x\in \bar{\Omega}. \end{aligned} \end{equation}

    In addition, from reference [38,Lemma 2.1], system (4.8) allows a unique globally attractive periodic solution \bar{n}^{*}(t, \cdot) . Then

    \begin{equation} \nonumber \begin{aligned} \liminf\limits_{t\rightarrow \infty} u_{1}(t,x,\bar{\varphi}_{1})\geq \varepsilon: = \min\limits_{t\in [0,\omega],x\in\bar{\Omega}}\bar{n}^{\ast}(t,x),\quad {\rm{uniformly}}\; {\rm{for}}\; x\in \bar{\Omega}. \end{aligned} \end{equation}

    This proof is complete.

    For every fixed \bar{\varphi}\in \mathcal{E}_{2} , u(t, x, \bar{\varphi}) is the unique solution of system (2.12) with u_{0} = \bar{\varphi} . The following conclusion indicates that R_{0} is a threshold for the transmission of disease.

    Theorem 4.8. Let R_{0} be defined in (4.4).We have the following statements:

    (i) When R_{0} < 1 , the disease-free \omega -periodic solution (u_{1}^{*}(t, x), 0, 0) of system (2.12) is globally attractive.

    (ii) When R_{0} > 1 , system (2.12) has at least one positive \omega -periodic solution (\bar{z}_{1}^{*}(t, x), \bar{z}_{2}^{*}(t, x), \bar{z}_{3}^{*}(t, x)) , and there exists a constant \bar{\sigma} > 0 such that for every \bar{\varphi}\in \mathcal{E}_{2} with \bar{\varphi}_{2}(\cdot)\not\equiv 0 , \bar{\varphi}_{3}(0, \cdot)\not\equiv 0 , on has

    \begin{equation} \begin{aligned} \liminf\limits_{t\rightarrow \infty}\min\limits_{x\in\bar{\Omega}}u_{i}(t,x,\bar{\varphi})&\geq \bar{\sigma}, \quad i = 1,2,3. \end{aligned} \end{equation} (4.9)

    Proof. Since \hat{Q}_{t} is the solution map of system (2.12) on \mathcal{E}_{2} , based on reference [22,Lemma 3.5], \hat{Q}_{t} is an \omega -periodic semiflow of system (2.12). Thence, \hat{Q}: = \hat{Q}_{\omega} is the Poincaré map, and \{\hat{Q}^{n}\}_{n\geq 0} is a discrete-time system on \mathcal{E}_{2} . For every fixed \bar{\varphi} \in \mathcal{E}_{2} , the omega limit set of the orbit \{\hat{Q}^{n}(\bar{\varphi})\}_{n\geq 0} is \omega(\bar{\varphi}) . Then \omega(\bar{\varphi}) is an internally chain transitive set for \{\hat{Q}^{n}\} on \mathcal{E}_{2} .

    (i) As R_{0} < 1 , according to Lemma 4.4, one has r(\tilde{Q}) < 1 , and therefore \tilde{\mu} = \frac{\ln r(\tilde{Q})}{\omega} < 0 . For t > 0 , discuss the following system with parameter \epsilon > 0 ,

    \begin{equation} \left\{ \begin{aligned} \frac{\partial w^{\epsilon}_{1}}{\partial t} = &D_{v}\Delta w^{\epsilon}_{1}-\mu_{v}(t,x)w^{\epsilon}_{1}+(1-\tau^{\prime}_{2}(t))\int_{\Omega}\Gamma\left(t,t-\tau_{2}(t),x,y\right) \\ &\frac{\alpha \beta(t-\tau_{2}(t),y)p[u_{1}^{*}(t-\tau_{2}(t),y)+\epsilon]w^{\epsilon}_{2}(t-\tau_{2}(t),y)}{lN(y)} {\rm{d}}y, &x\in \Omega,\\ \frac{\partial w^{\epsilon}_{2}}{\partial t} = &D_{h}\Delta w_{2}^{\epsilon}+c\beta(t,x)w^{\epsilon}_{1}-(\mu_{h}(x)+\varrho(x))w^{\epsilon}_{2}, & x\in \Omega,\\ \frac{\partial w^{\epsilon}_{1}}{\partial \nu} = &\frac{\partial w^{\epsilon}_{2}}{\partial \nu} = 0, & x\in \partial\Omega. \end{aligned} \right. \end{equation} (4.10)

    Denote w^{\epsilon}(t, x, \hat{\varphi}) = (w^{\epsilon}_{1}(t, x, \hat{\varphi}), w^{\epsilon}_{2}(t, x, \hat{\varphi})) to be the unique solution of system (4.10) for any \hat{\varphi}\in \mathcal{E}_{1} , with w^{\epsilon}_{0}(\hat{\varphi})(\theta, x) = (\hat{\varphi}(0, x), \hat{\varphi}(\theta, x)) for all \theta\in [-\tau_{2}(0), 0] , x\in \bar{\Omega} , where

    \begin{equation} \nonumber \begin{aligned} w^{\epsilon}_{t}(\hat{\varphi})(\theta,x) = &w^{\epsilon}(t+\theta,x,\hat{\varphi}) = (w^{\epsilon}_{1}(t,x,\hat{\varphi}), w^{\epsilon}_{2}(t+\theta,x,\hat{\varphi})),\; \theta\in [-\tau_{2}(0),0], \forall t\geq 0, x\in \bar{\Omega}. \end{aligned} \end{equation}

    Let \tilde{Q}_{t}^{\epsilon}: = \mathcal{E}_{1}\rightarrow \mathcal{E}_{1} be the solution map of system (4.10) and \tilde{Q}^{\epsilon}: = \tilde{Q}_{\omega}^{\epsilon} be the corresponding Poincaré map, i.e., \tilde{Q}^{\epsilon}(\hat{\varphi}) = w^{\epsilon}_{\omega}(\hat{\varphi}) , for \forall \hat{\varphi}\in \mathcal{E}_{1} , and make r(\tilde{Q}^{\epsilon}) be the spectral radius of \tilde{Q}^{\epsilon} . Then the eventual strong monotonicity and compactness of \tilde{Q}^{\epsilon} on \mathcal{E}_{1} are easy to prove, for t\geq 3\hat{\tau}_{2} . From the continuity of the principle eigenvalue to parameters, one has \lim\limits_{\epsilon\rightarrow 0}r(\tilde{Q}^{\epsilon}) = r(\tilde{Q}) . Allow \epsilon > 0 to be small enough, one has r(\tilde{Q}^{\epsilon}) < 1 . From Lemma 4.5, a positive \omega -periodic function w^{*}_{\epsilon}(t, x) exists such that w^{\epsilon}(t, x) = e^{\tilde{\mu}_{\epsilon } t}w^{*}_{\epsilon}(t, x) is a solution of system (4.10), where \tilde{\mu}_{\epsilon} = \frac{\ln r(\tilde{Q}^{\epsilon})}{\omega} < 0 . By system (4.1), u_{1}^{*}(t, \cdot) is globally attractive and the comparison principle, there must be a large enough integer n_{1} > 0 such that n_{1}\omega > \hat{\tau}_{2} and

    \begin{equation} \nonumber \begin{aligned} u_{1}(t,x)\leq u_{1}^{*}(t,x)+\epsilon,\qquad \forall t\geq n_{1}\omega-\hat{\tau}_{2}, \quad x\in \bar{\Omega}. \end{aligned} \end{equation}

    Thus, as t\geq n_{1}\omega , one has

    \begin{equation} \left\{ \begin{aligned} \frac{\partial u_{2}}{\partial t}\leq &D_{v}\Delta u_{2}-\mu_{v}(t,x)u_{2}+(1-\tau^{\prime}_{2}(t))\int_{\Omega}\Gamma\left(t,t-\tau_{2}(t),x,y\right)\cdot\\ &\frac{\alpha \beta(t-\tau_{2}(t),y)p[u_{1}^{*}(t-\tau_{2}(t),y)+\epsilon]u_{3}(t-\tau_{2}(t),y)}{lN(y)} {\rm{d}}y, &x\in \Omega,\\ \frac{\partial u_{3}}{\partial t}\leq &D_{h}\Delta u_{3}+c\beta(t,x)u_{2}-(\mu_{h}(x)+\varrho(x))u_{3}, &x\in \Omega,\\ \frac{\partial u_{2}}{\partial \nu} = &\frac{\partial u_{3}}{\partial \nu} = 0, &x\in \partial\Omega, \end{aligned} \right. \end{equation} (4.11)

    there exists some C_{1} > 0 , for any given \bar{\varphi}\in \mathcal{E}_{2} , such that

    \begin{equation} \nonumber \begin{aligned} (u_{2}(t,x,\bar{\varphi}),u_{3}(t,x,\bar{\varphi}))\leq C_{1}(w_{1}^{\epsilon}(t,x),w_{2}^{\epsilon}(t,x)), \quad t\in [n_{1}\omega-\hat{\tau}_{2}, n_{1}\omega],\; x\in \bar{\Omega}. \end{aligned} \end{equation}

    Using systems (4.10) and (4.11), and the comparison theorem [36,Proposition 1], we have

    \begin{equation} \nonumber \begin{aligned} \left(u_{2}(t,x,\bar{\varphi}),u_{3}(t,x,\bar{\varphi})\right)&\leq C_{1}e^{\tilde{\mu}_{\epsilon}t} w^{*}_{\epsilon}(t,x), \; \; \; \forall t\geq n_{1}\omega, \; \; \; x\in \bar{\Omega}. \end{aligned} \end{equation}

    Accordingly, \lim\limits_{t\rightarrow \infty}\left(u_{2}(t, x, \bar{\varphi}), u_{3}(t, x, \bar{\varphi})\right) = (0, 0) uniformly for x\in \bar{\Omega} . Using the theory of internally chain transitive sets [32,Section 1.2], we can prove that \lim\limits_{t\rightarrow \infty}(u_{1}(t, x, \bar{\varphi})-u_{1}^{*}(t, x)) = 0 uniformly for x\in\bar{\Omega} , where u_{1}^{*}(t, x) is a globally attractive solution of system (4.1). Based on the above argument, u_{1}(t, \cdot, \bar{\varphi}) meets a non-autonomous system, and this system tends to periodic system (4.1) asymptotically. It is known that system (4.1) can produce a solution semiflow \hat{P}_{t} , t\geq 0 on C([-\tau_{2}(0), 0], \mathbb{Y}^{+}) . Then \hat{P}: = \hat{P}_{\omega} is a Poincaré map related to system (4.1). Obviously, \hat{P} has a global attractor in C([-\tau_{2}(0), 0], \mathbb{Y}^{+}) .

    Let \mathcal{J} = \omega(\bar{\varphi}) be the omega limit set of \bar{\varphi} = (\bar{\varphi}_{1}, \bar{\varphi}_{2}, \bar{\varphi}_{3})\in \mathcal{E}_{2} for \hat{Q}_{t} . Owing to \lim\limits_{t\rightarrow \infty}u_{2}(t, x, \bar{\varphi}) = 0 and \lim\limits_{t\rightarrow \infty}u_{3}(t, x, \bar{\varphi}) = 0 uniformly for x\in\bar{ \Omega} , one has \mathcal{J} = J\times\{0\}\times \{\hat{0}\} . Here, \hat{0} represents function which is identical to 0 , i.e., \hat{0}(\theta, \cdot) = 0 , \forall \theta\in [-\tau_{2}(0), 0] . Following Lemma 4.7, we have \hat{0}\not\in J . Based on [32,Lemma 1.2.1], we know that \mathcal{J} is an internally chain transitive set for \hat{Q} , then J is an internally transitive chain set for \hat{P} . Define u^{0}\in C([-\tau_{2}(0), 0], \mathbb{Y}^{+}) by u^{0}(\theta, \cdot) = u_{1}^{*}(\theta, \cdot) for \theta\in [-\tau_{2}(0), 0] . Since J\neq \{\hat{0}\} and u^{0} is globally attractive in C\left([-\tau_{2}(0), 0], \mathbb{Y}^{+}\right)\backslash \{\hat{0}\} , one has J\cap W^{s}(u^{0})\neq \emptyset , where W^{s}(u^{0}) is the stable set of u^{0} . According to [32,Lemma 1.2.1], one has J = \{u^{0}\} . \mathcal{J} = \{u^{0}, 0, \hat{0}\} is proved, thence

    \begin{equation} \nonumber \begin{aligned} \lim\limits_{t\rightarrow \infty}\| (u_{1}(t,x,\bar{\varphi}),u_{2}(t,x,\bar{\varphi}),u_{3}(t,x,\bar{\varphi}))-(u_{1}^{*}(t,\cdot),0,0)\| = 0. \end{aligned} \end{equation}

    (ii) When R_{0} > 1 , Lemmas 4.1, 4.4 and 4.5 indicate r(\tilde{Q}) > 1 , consequently, \tilde{\mu} = \frac{\ln r(\tilde{Q})}{\omega} > 0 . Let

    \begin{equation} \nonumber \begin{aligned} \mathbb{P}_{0} = \{\bar{\varphi}\in \mathcal{E}_{2}:\bar{\varphi}_{2}(\cdot)\not \equiv 0\; {\rm{and}}\; \bar{\varphi}_{3}(0,\cdot)\not \equiv0\}, \end{aligned} \end{equation}

    and

    \begin{equation} \nonumber \begin{aligned} \partial\mathbb{P}_{0} = \mathcal{E}_{2}\backslash \mathbb{P}_{0} = \{\bar{\varphi}\in \mathcal{E}_{2}:\bar{\varphi}_{2}(\cdot)\equiv0\; {\rm{or}}\; \bar{\varphi}_{3}(0,\cdot)\equiv 0\}. \end{aligned} \end{equation}

    For every \bar{\varphi}\in \mathbb{P}_{0} , Lemma 4.7 reveals that u_{2}(t, x, \bar{\varphi}) > 0 and u_{3}(t, x, \bar{\varphi}) > 0 , \forall t > 0 , x\in \bar{\Omega} . Thus, \hat{Q}^{n}(\mathbb{P}_{0})\subset \mathbb{P}_{0} , \forall n\in \mathbb{N} . Lemma 4.6 evidences that \hat{Q} has a strong global attractor in \mathcal{E}_{2} .

    Set

    \begin{equation} \nonumber \begin{aligned} M_{\partial}: = \{\bar{\varphi}\in \partial\mathbb{P}_{0}: \hat{Q}^{n}(\bar{\varphi})\in \partial\mathbb{P}_{0}, \forall n\in \mathbb{N}\}, \end{aligned} \end{equation}

    and the omega limit set of the orbit \gamma^{+}: = \{\hat{Q}^{n}(\bar{\varphi}):\forall n\in \mathbb{N}\} is \omega(\bar{\varphi}) . Denote H = (u_{1}^{*}(t, x), 0, \hat{0}) . Next, we prove that H cannot form a cycle for \hat{Q} in \partial\mathbb{P}_{0} .

    Claim 1. For every \tilde{\phi}\in M_{\partial} , the omega limit set \omega(\tilde{\phi}) = H .

    For every fixed \tilde{\phi}\in M_{\partial} , \hat{Q}^{n}(\tilde{\phi})\in \partial\mathbb{P}_{0} , \forall n\in \mathbb{N} . Then either u_{2}(n\omega, \cdot, \tilde{\phi})\equiv 0 or u_{3}(n\omega, \cdot, \tilde{\phi})\equiv 0 , for each n\in \mathbb{N} . Furthermore, by contradiction and Lemma 4.7, it is obvious that for each t\geq0 , either u_{2}(t, \cdot, \tilde{\phi})\equiv 0 or u_{3}(t, \cdot, \tilde{\phi})\equiv 0 . If u_{2}(t, \cdot, \tilde{\phi})\equiv 0 for t\geq 0 , then \lim\limits_{t\rightarrow \infty}\left(u_{1}(t, x, \tilde{\phi})-u_{1}^{*}(t, x)\right) = 0 uniformly for x\in \bar{\Omega} . Thereby, the u_{3} equation in system (2.12) meets

    \begin{equation} \nonumber \begin{aligned} \frac{\partial u_{3}}{\partial t}\leq D_{h}\Delta u_{3}-(\bar{\mu}_{h}+\bar{\varrho})u_{3}, \end{aligned} \end{equation}

    where \bar{\varrho} = \min\limits_{x\in \bar{\Omega}}\varrho(x) and \bar{\mu}_{h} = \min\limits_{x\in \bar{\Omega}}\mu_{h}(x) . According to the comparison principle, one has \lim\limits_{t\rightarrow \infty}u_{3}(t, x, \tilde{\phi}) = 0 uniformly for x\in\bar{\Omega} . Suppose that u_{2}(t_{2}, \cdot, \tilde{\phi})\neq 0 for some t_{2}\geq 0 , in the light of Lemma 4.7, we obtain u_{2}(t, \cdot, \tilde{\phi}) > 0 , \forall t\geq t_{2} . Accordingly, u_{3}(t, \cdot, \tilde{\phi})\equiv 0 . From the u_{2} equation in system (2.12), one has \lim\limits_{t\rightarrow \infty}u_{2}(t, x, \tilde{\phi}) = 0 uniformly for x\in \bar{\Omega} . Consequently, the u_{1} equation in system (2.12) abides by a nonautonomous system, which is asymptomatic to the periodic system (4.1). The theory of asymptomatically periodic system [32,Section 3.2] implies that \lim\limits_{t\rightarrow \infty}(u_{1}(t, x, \tilde{\phi})-u_{1}^{*}(t, x)) = 0 uniformly for x\in \bar{\Omega} . Therefore, \omega(\tilde{\phi}) = H for any \tilde{\phi}\in M_{\partial} , and H can not form a cycle for \hat{Q} in \partial\mathbb{P}_{0} .

    For t > 0 , study the time-periodic parabolic system with parameter \bar{\delta} > 0

    \begin{equation} \left\{ \begin{aligned} \frac{\partial w^{\bar{\delta}}_{1}}{\partial t} = &D_{v}\Delta w^{\bar{\delta}}_{1}-\mu_{v}(t,x)w^{\bar{\delta}}_{1}+(1-\tau^{\prime}_{2}(t))\int_{\Omega} \Gamma(t,t-\tau_{2}(t),x,y) \cdot\\ &\frac{\alpha\beta(t-\tau_{2}(t),y)p\left[u_{1}^{*}(t-\tau_{2}(t),y)-\bar{\delta}\right]w^{\bar{\delta}}_{2}(t-\tau_{2}(t),y)}{p\bar{\delta} +lN(y)} {\rm{d}}y, &x\in \Omega,\\ \frac{\partial w^{\bar{\delta}}_{2}}{\partial t} = &D_{h}\Delta w^{\bar{\delta}}_{2}+\frac{c\beta(t,x)l\left[N(x)-\bar{\delta}\right]w^{\bar{\delta}}_{1}}{p\bar{\delta}+lN(x)}-(\mu_{h}(x)+\varrho(x))w^{\bar{\delta}}_{2},&x\in \Omega,\\ \frac{\partial w^{\bar{\delta}}_{1}}{\partial \nu} = &\frac{\partial w^{\bar{\delta}}_{2}}{\partial \nu} = 0, &x\in \partial\Omega. \end{aligned} \right. \end{equation} (4.12)

    For every \hat{\varphi}\in \mathcal{E}_{1} , set w^{\bar{\delta}}(t, x, \hat{\varphi}) = (w^{\bar{\delta}}_{1}(t, x, \hat{\varphi}), w^{\bar{\delta}}_{2}(t, x, \hat{\varphi})) to be the unique solution of system (4.12) with w^{\bar{\delta}}_{0}(\hat{\varphi})(\theta, x) = (\hat{\varphi}_{1}(0, x), \hat{\varphi}_{2}(\theta, x)) for all \theta\in [-\tau_{2}(0), 0] , x\in \bar{\Omega} , where

    \begin{equation} \nonumber \begin{aligned} w^{\delta}_{t}(\hat{\varphi})(\theta,x) = w^{\bar{\delta}}(t+\theta,x,\hat{\varphi}) = \left(w^{\bar{\delta}}_{1}(t,x,\hat{\varphi}),w^{\bar{\delta}}_{2}(t+\theta,x,\hat{\varphi})\right), \end{aligned} \end{equation}

    for any t\geq 0 and x\in \bar{\Omega} . Let \tilde{Q}_{t}^{\bar{\delta}}:\mathcal{E}_{1}\rightarrow \mathcal{E}_{1} be the solution map of system (4.12), i.e., \tilde{Q}_{t}^{\bar{\delta}}(\hat{\varphi}) = w^{\bar{\delta}}_{t}(\hat{\varphi}) , and \tilde{Q}^{\bar{\delta}}: = \tilde{Q}_{\omega}^{\bar{\delta}} be the corresponding Poincaré map, for \forall \hat{\varphi}\in \mathcal{E}_{1} . Let r(\tilde{Q}^{\bar{\delta}}) be the spectral radius of \tilde{Q}^{\bar{\delta}} . Because \lim\limits_{\bar{\delta}\rightarrow 0}r(\tilde{Q}^{\bar{\delta}}) = r(\tilde{Q}) > 1 , a fully small \bar{\delta} > 0 can be selected, such that

    \begin{equation} \nonumber \begin{aligned} \bar{\delta} < \min\left\{\min _{t\in [0,\omega],x\in \bar{\Omega}}u_{1}^{*}(t,x),\; \; \min\limits_{x\in \bar{\Omega}}N(x)\right\}, \quad {\rm{and}} \quad r(\tilde{Q}^{\bar{\delta}}) > 1. \end{aligned} \end{equation}

    For the certain \bar{\delta} > 0 , according to the continuous dependence of the solution on the initial data, there must be a constant \delta^{*} > 0 such that for all \bar{\varphi} with \|\bar{\varphi} -H\| < \delta^{*} , which leads to \|\hat{Q}_{t}(\bar{\varphi}) -\hat{Q}_{t}(H)\| < \bar{\delta} for all t\in [0, \omega] .

    Claim 2. For any \bar{\varphi}\in \mathbb{P}_{0} , then \limsup\limits_{n\rightarrow \infty}\|\hat{Q}^{n}(\varphi)-H\|\geq \delta^{*} .

    Assuming for some \bar{\varphi}_{0}\in \mathbb{P}_{0} , one has \limsup\limits_{n\rightarrow \infty}\|\hat{Q}^{n}(\bar{\varphi}_{0})-H\| < \delta^{*} , by contradiction. So there is n_{2}\geq 1 such that \|\hat{Q}^{n}(\bar{\varphi}_{0})-H\| < \delta^{*} for all n\geq n_{2} . For each t\geq n_{2}\omega , set t = n\omega+t^{\prime} with n = [t/ \omega] and t^{\prime}\in [0, \omega) , one has

    \begin{equation} \begin{aligned} \|\hat{Q}_{t}(\bar{\varphi}_{0})-\hat{Q}_{t}(H)\| = \|\hat{Q}_{t^{\prime}}(\hat{Q}^{n}(\bar{\varphi}_{0})) -\hat{Q}_{t^{\prime}}(H)\| < \bar{ \delta}. \end{aligned} \end{equation} (4.13)

    Following system (4.13) and Lemma 4.7, we can receive

    \begin{equation} \nonumber \begin{aligned} u_{1}^{*}(t,x)-\bar{\delta} < u_{1}(t,x,\bar{\varphi}_{0})& < u_{1}^{*}(t,x)+\bar{\delta},\quad 0 < u_{2}(t,x,\bar{\varphi}_{0}) < \bar{\delta},\quad 0 < u_{3}(t,x,\bar{\varphi}_{0}) < \bar{\delta}, \end{aligned} \end{equation}

    for each t > n_{2}\omega-\hat{\tau}_{2} , and x\in\bar{\Omega} . As t\geq n_{2}\omega , u_{2}(t, x, \bar{\varphi}_{0}) and u_{3}(t, x, \bar{\varphi_{0}}) meet

    \begin{equation} \left\{ \begin{aligned} \frac{\partial u_{2}}{\partial t}\geq& D_{v}\Delta u_{2}-\mu_{v}(t,x) u_{2}+(1-\tau^{\prime}_{2}(t)) \int_{\Omega}\Gamma\left(t,t-\tau_{2}(t),x,y\right)\cdot\\ & \frac{\alpha\beta(t-\tau_{2}(t),y)p\left[u_{1}^{*}(t-\tau_{2}(t),y)-\bar{\delta}\right]u_{3}(t-\tau_{2}(t),y)} {p\bar{\delta}+lN(y)} {\rm{d}}y, &x\in \Omega, \\ \frac{\partial u_{3}}{\partial t}\geq& D_{h}\Delta u_{3}+\frac{c\beta(t,x)l\left[N(x)-\bar{\delta}\right]u_{2}}{p\bar{\delta}+lN(x)}-(\mu_{h}(x)+\varrho(x))u_{3}, & x\in \Omega,\\ \frac{\partial u_{2}}{\partial \nu} = &\frac{\partial u_{3}}{\partial \nu} = 0, &x\in \partial\Omega. \end{aligned} \right. \end{equation} (4.14)

    Based on u(t, x, \bar{\varphi}_{0})\gg 0 for every t\geq 0 and x\in \bar{\Omega} , there must be a constant C_{2} > 0 , such that

    \begin{equation} \nonumber \begin{aligned} (u_{2}(t,x,\bar{\varphi}_{0}),u_{3}(t,x,\bar{\varphi}_{0}))\geq C_{2}e^{\tilde{\mu}_{\bar{\delta}}t}u^{*}_{\bar{\delta}}(t,x), \forall t\in [n_{2}\omega-\hat{\tau}_{2},n_{2}\omega], x\in \bar{\Omega}, \end{aligned} \end{equation}

    where u^{*}_{\bar{\delta}}(t, x) is a positive \omega -periodic function and e^{\tilde{\mu}_{\bar{\delta}}t}u^{*}_{\bar{\delta}}(t, x) is a solution of system (4.12), then \tilde{\mu}_{\bar{\delta}} = \frac{\ln r(\tilde{Q}^{\bar{\delta}})}{\omega} > 1 . By system (4.14) and the comparison theorem

    \begin{equation} \nonumber \begin{aligned} (u_{2}(t,x,\bar{\varphi}_{0}),u_{3}(t,x,\bar{\varphi}_{0}))\geq C_{2}e^{\tilde{\mu}_{\bar{\delta}}t}u^{*}_{\bar{\delta}}(t,x), \quad \forall t\geq n_{2}\omega,\quad x\in\bar{\Omega}. \end{aligned} \end{equation}

    Obviously, \tilde{\mu}_{\bar{\delta}} > 0 , then u_{2}(t, \cdot, \bar{\varphi}_{0})\rightarrow +\infty and u_{3}(t, \cdot, \bar{\varphi}_{0})\rightarrow +\infty as t\rightarrow +\infty , which leads to a contradiction.

    With the above claim, it is easy to obtain that H is an isolated invariant set for \hat{Q} in \mathcal{E}_{2} , and W^{s}(H)\cap \mathbb{P}_{0} = \emptyset , where W^{s}(H) is the stable set of H for \hat{Q} . For every integer n and n\omega > \hat{\tau}_{2} , \hat{Q}^{n}: = \hat{Q}_{n\omega} is compact, then \hat{Q} is compact, and \hat{Q} is asymptotically smooth on \mathcal{E}_{2} . Moreover, similar to Lemma 4.6, we get that \hat{Q} has a global attractor on \mathcal{E}_{2} . According to reference [45,Theorem 3.7], \hat{Q} has a global attractor A_{0} in \mathbb{P}_{0} . Based on the acyclicity theorem on uniform persistence for map [32,Theorem 1.3.1 and Remark 1.3.1], we obtain that \hat{Q} is uniformly persistent with respect to (\mathbb{P}_{0}, \partial \mathbb{P}_{0}) , in this sense, there exists an \gamma_{1} > 0 , such that

    \begin{equation} \begin{aligned} \liminf\limits_{n\rightarrow \infty} {\rm{d}}(\hat{Q}^{n}(\bar{\varphi}),\partial \mathbb{P}_{0})\geq \gamma_{1},\quad \forall \bar{\varphi} \in \mathbb{P}_{0}. \end{aligned} \end{equation} (4.15)

    Now, the practical persistence are derived. Following A_{0} = \hat{Q}(A_{0}) , we obtain that \bar{\varphi}_{2}(\cdot) > 0 and \bar{\varphi}_{3}(0, \cdot) > 0 for any \bar{\varphi}\in A_{0} . Set \mathcal{B}: = \bigcup\limits_{t\in[0, \omega]}\hat{Q}_{t}(A_{0}) . Obviously, we have \mathcal{B}\subset\mathbb{P}_{0} and \lim\limits_{t\rightarrow \infty} {\rm{d}}(\hat{Q}_{t}(\bar{\varphi}), \mathcal{B}) = 0 , \forall \bar{\varphi}\in \mathbb{P}_{0} . The continuous function q:\mathcal{E}_{2}\rightarrow \mathbb{R}_{+} is defined by

    \begin{equation} \nonumber \begin{aligned} q(\bar{\varphi}): = \min\left\{\min\limits_{x\in\bar{\Omega}}\bar{\varphi}_{2}(0,x),\min\limits_{x\in \bar{\Omega}}\bar{\varphi}_{3}(0,x)\right\},\quad \forall \bar{\varphi}\in (\bar{\varphi}_{1},\bar{\varphi}_{2},\bar{\varphi}_{3})\in \mathcal{E}_{2}. \end{aligned} \end{equation}

    Because \mathcal{B} is a compact subset of \mathbb{P}_{0} , we obtain that \inf\limits_{\bar{\varphi} \in \mathcal{B}}q(\bar{\varphi}) = \min\limits_{\bar{\varphi}\in \mathcal{B}}q(\bar{\varphi}) > 0 . Consequently, there is an \gamma_{2} > 0 such that

    \begin{equation} \nonumber \begin{aligned} \liminf\limits_{t\rightarrow \infty} q(\hat{Q}_{t}(\bar{\varphi})) = \liminf\limits_{t\rightarrow \infty}\min\left(\min\limits_{x\in \bar{\Omega}} u_{2}(t,x,\bar{\varphi}), \min\limits_{x\in\bar{\Omega}} u_{3}(t,x,\bar{\varphi})\right)\geq \gamma_{2}, \; \forall \bar{\varphi} \in \mathbb{P}_{0}. \end{aligned} \end{equation}

    Furthermore, according to Lemma 4.7, there must be an \gamma_{3}\in (0, \gamma_{2}) such that

    \begin{equation} \nonumber \begin{aligned} \liminf\limits_{t\rightarrow \infty}\min\limits_{x\in\bar{\Omega}}u_{j}(t,x, \bar{\varphi}) \geq \gamma_{3}, \quad \forall\bar{ \varphi} \in \mathbb{P}_{0}\; \; (j = 1,2,3). \end{aligned} \end{equation}

    By references [29,Lemma 8] and [32,Theorem 3.5.1], for every t > 0 , the solution map Q_{t}: = W_{h}\rightarrow W_{h} of system (2.12), defined in Lemma 3.2, is a \kappa -contraction about an equivalent norm on W_{h} , where \kappa is Kuratowski measure of noncompactness. Define

    \begin{equation} \nonumber \begin{aligned} \mathbb{W}_{0} = \left\{\phi\in W_{h}:\phi_{2}(0,\cdot)\not\equiv 0 \; \; {\rm{and}} \; \; \phi_{3}(0,\cdot)\not\equiv 0\right\}, \end{aligned} \end{equation}

    and

    \begin{equation} \nonumber \begin{aligned} \partial\mathbb{W}_{0} = W_{h} /\mathbb{W}_{0} = \left\{\phi\in W_{h}: \phi_{2}(0,\cdot) = 0 \; \; {\rm{or}} \; \; \phi_{3}(0,\cdot) = 0\right\}. \end{aligned} \end{equation}

    For any \phi\in \mathbb{W}_{0} , Q is \rho -uniformly persistent with \rho(\phi) = {\rm{d}}(\phi, \partial \mathbb{W}_{0}) is easily attainable. It then follows from [45,Theorem 4.5], applicable to Q , that system (2.12) has an \omega -periodic solution (z_{1}^{*}(t, \cdot), z_{2}^{*}(t, \cdot), z_{3}^{*}(t, \cdot)) with (z_{1t}^{*}, z_{2t}^{*}, z_{3t}^{*})\in \mathbb{W}_{0} . Let \bar{z}_{1}^{*}(\theta, \cdot) = z_{1}^{*}(\theta, \cdot) , \bar{z}_{2}^{*}(0, \cdot) = \bar{z}_{2}^{*}(0, \cdot) and \bar{z}_{3}^{*}(\theta, \cdot) = z_{3}^{*}(\theta, \cdot) where \theta\in [-\tau_{2}(0), 0] . Combing the uniqueness of solutions and Lemma 4.7, we have that (\bar{z}_{1}^{*}(t, \cdot), \bar{z}_{2}^{*}(t, \cdot), \bar{z}_{3}^{*}(t, \cdot)) is periodic solution system (2.12) and it is also strictly positive.

    Next, we prove the asymptomatic behaviour of M and E_{v} in system (2.11). By Eq (2.6), and \Lambda_{v}(t, x)\not\equiv 0 is nonnegative and Hölder continuous function for t > 0 and x\in\bar{\Omega} . Then there is a unique positive \omega -periodic solution M^{*}(t, x) which is globally attractive. From Eq (3.6) and Theorem 4.8(i), when R_{0} < 1 , one has \lim\limits_{t\rightarrow \infty}E_{v}(t, \cdot) = 0 . With the arguments similar to those above, we can prove that if R_{0} < 1 , then \lim\limits_{t\rightarrow \infty}(E_{v}(t, x), I_{v}(t, x), I_{h}(t, x)) = (0, 0, 0) and \lim\limits_{t\rightarrow \infty}\left((M(t, \cdot), S_{v}(t, x))-(M^{*}(t, x), u_{1}^{*}(t, x))\right) = 0 , uniformly for x\in \bar{\Omega} . By Eq (2.12), Theorem (4.8)(ii) and Eqs (2.6) and (3.6), then if R_{0} > 1 , for any \bar{\varphi}\in \mathcal{E}_{2} with \bar{\varphi}_{2}(\cdot)\not\equiv 0 and \bar{\varphi}_{3}(0, \cdot)\not\equiv 0 , one has

    \begin{equation} \nonumber \begin{aligned} \liminf\limits_{t\rightarrow \infty}S_{v}(t,x)\geq \gamma_{3},\quad \liminf\limits_{t\rightarrow \infty}I_{v}(t,x)\geq \gamma_{3},\quad \liminf\limits_{t\rightarrow \infty}I_{h}(t,x)\geq \gamma_{3}, \end{aligned} \end{equation}

    uniformly for x\in\bar{\Omega} . Since \Lambda_{v}(t, x) is Hölder continuous and non-negative function for t > 0 and x\in \Omega from Eq (2.6), then there exists \gamma_{4} > 0 such that

    \begin{equation} \nonumber \begin{aligned} \liminf\limits_{t\rightarrow \infty}M(t,x)\geq \gamma_{4}, \end{aligned} \end{equation}

    uniformly for x\in\bar{\Omega} with M(0, x) satisfies compatibility condition (3.4). By the integral form Eq (3.6), then exists \gamma_{5} > 0 such that

    \begin{equation} \nonumber \begin{aligned} \liminf\limits_{t\rightarrow \infty}\min\limits_{x\in\bar{\Omega}}E_{v}(t,x) \geq \gamma_{5}, \end{aligned} \end{equation}

    with E_{v}(0, x) satisfies compatibility condition (3.5). In addition, if (S_{v}(t, \cdot), I_{v}(t, \cdot), I_{h}(t, \cdot)) is \omega -periodic in t , then E_{v}(t, \cdot) is \omega -periodic is easy to be checked. Therefore, the following result is obtained.

    Theorem 4.9. Let R_{0} be defined as in (4.4).Set U(t, x, \psi) to be the solution of system (2.11) with U_{0} = \psi , where (U_{1}(t, x), U_{2}(t, x), U_{3}(t, x), U_{4}(t, x), U_{5}(t, x)) = (M(t, x), S_{v}(t, x), E_{v}(t, x), I_{v}(t, x), I_{h}(t, x)) .One has:

    (i) When R_{0} < 1 , the disease-free \omega -periodic solution (M^{*}(t, \cdot), u_{1}^{*}(t, \cdot), 0, 0, 0) of system (2.11) is globally attractive.

    (ii) When R_{0} > 1 , system (2.11) has at least one positive \omega -periodic solution \bar{U}(t, x) , where \bar{U}_{2}(t, x) = \bar{z}_{1}^{\ast}(t, x) , \bar{U}_{4}(t, x) = \bar{z}_{2}^{\ast}(t, x) , \bar{U}_{5}(t, x) = \bar{z}_{3}^{\ast}(t, x) and there is a constant \gamma_{6} = \max\{\gamma_{3}, \gamma_{4}, \gamma_{5}\} > 0 that makes for every \psi\in \mathcal{D} with \psi_{4}(\cdot)\not\equiv 0 , \psi_{5}(0, \cdot)\not\equiv 0 , we have

    \begin{equation} \begin{aligned} \liminf\limits_{t\rightarrow \infty}\min\limits_{x\in\bar{\Omega}}\bar{U}_{i}(t,x,\psi)\geq \gamma_{6}, \quad i = 1,2,3,4,5,6. \end{aligned} \end{equation} (4.16)

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

    Here, some simulations to support the analysis results built in Section 4 are provided. More precisely, the potential dynamics outcomes of system (2.11) are given. For convenience, we will focus our attention on the domain \Omega , which is one-dimensional, denoted by [0, \pi] to simulate the long-time behavior of system (2.11). The basic reproduction number R_{0} often (but not always) controls the extinction and persistence of the disease. Specifically, system (2.11) applies to the transmission of malaria in Maputo Province, Mozambique, a sub-Saharan African country where malaria is at risk of spreading in the country. The WHO report shows that the number of confirmed malaria cases has been increasing in recent years, and the condition deteriorated severely in 2014.

    According to reference [24], suppose the total population density N(x) = \bar{N} = 53\; ({\rm{km}}^{2})^{-1} . The unit of time is month. Set the periodic \omega = 1\; {\rm{year}} = 12 months. We set \varrho = a_{1}(2.05-\cos(2x))\; {\rm{Month}}^{-1} , where a_{1} = 0.0187 means that people living in urban areas (near the center area) can get better medical services (number of doctors and hospitals, advanced medical equipment, supply of medicines) than people in rural areas. This results in a higher recovery rate near the center. See Appendix E for the values and definitions of other parameters.

    To compute R_{0} , we use the numerical method proposed in reference [41,Lemma 2.5 and Remark 3.2]. Using the above parameters set, the basic reproduction number can be numerically calculate as R_{0} = 1.2233 > 1 , which implies that in mosquito and human populations, the disease persist. Figure 3 shows numerical plots of humans and mosquitoes compartments with initial data

    \hat{u}(\theta,x) = \left(\begin{array}{c} 500-5\cos2x\\ 300-5\cos2x\\ 32-5\cos2x\\ 40-5\cos2x\\ 5-2\cos2x \end{array}\right), \; \; \forall \theta \in [-\hat{\tau},0],\; \; \; x \in [0,\pi],

    where \hat{u}(t, x): = (M, S_{v}, E_{v}, I_{v}, I_{h}) . This conforms to the conclusion of Theorem 4.9(ii). The time interval can be truncated [20,30] to prove that the existence of the positive periodic solution.

    Figure 3.  The long-term behavior of the solutions of system (2.11) when R_{0} > 1 .

    Reduce the bite rate to 0.5\beta(t) , and increase the death rate of mosquito to 1.5\mu_{v}(t) , which can be achieved by using insecticide-treated nets, then R_{0} = 0.2983 < 1 . Thus, immature mosquitoes and susceptible adult mosquitoes show periodic fluctuations. The exposed adult mosquitoes, infectious adult mosquitoes and infectious humans converge to 0 (Figure 4). This is coincident with the conclusion of Theorem 4.9(i).

    Figure 4.  The long-term behavior of system (2.11) solution, when R_{0} < 1 .

    We are interested in how system parameters affect the disease risk R_{0} . In this section, we choose the parameters that are consistent with those in Section 5.1.

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

    Figure 5.  The relationship between R_{0} and parameters. There are two curves. The green one represents the time-average maturity delay [\tau_{1}] , and the red one represents \tau_{1}(t) affected by the periodic temperature.

    Looking at Figure 5, we let l/p vary in [0, 1] , and other parameters are same as those in Figure 3 and explore the impact of the vector-bias effect, l/p are used to represent the relative attractiveness of susceptible to infected humans. From Figure 5(a), it can be easily seen that R_{0} is an increasing function of l/p . Therefore, ignoring the vector-bias effect will ultimately underestimate the risk of disease transmission. And this conclusion is valid regardless of whether the periodic or the time-average maturation delay is considered.

    In this subsection, we examine the impact of interventions on R_{0} . Note that in Figure 3(e), the population density of infected individuals in urban (central) areas is less than in rural (border) areas, which prompts us to explore the relationship between human recovery rate \varrho and disease risk R_{0} in a heterogeneous space. Since medical resources (advanced medical equipment, the number of hospitals and doctors, medicine supply) are abundant in urban areas than in rural areas, humans can get better medical services in urban areas with a higher recovery rate.

    In this subsection, we analyze the impact of medical resources on the spread of malaria, including two aspects: (i) the impact of increasing medical resources; (ii) for fixed medical resources, the impact of different allocation measures on disease transmission. In Figure 5(b), we let a_{1} in \varrho vary in [0.02, 0.5] , and the other parameters are same as those in Figure 3. It proposes a strategy, i.e., in the current allocation of the medical resources, the recovery rate can be improved through efforts to develop new drugs. Figure 5(b) shows that increasing medical resources can reduce R_{0} to less than 1, which means that the disease can be controlled. In this paper, we use the recovery rate \int_{0}^{\pi}\varrho(x) {\rm{d}}x of the entire region to represent medical resources.

    In the case of fixed medical resources, we explore the impact of different allocation on the spread of disease. Choose a_{1} = 0.058 and a new parameter \delta_{1} is added to \varrho , i.e., \varrho(x) = 0.058\times (1 -(1-\delta_{1})\cos(2x)) , where \delta_{1}\in [0, 1] , which measures the difference in the allocation of urban and rural medical resources. Other parameters are the same as those in Figure 3. When \delta_{1} = 0 , the gap between rural and urban medical resources is the largest, and the recovery rate in urban area is the highest (near the center of space area, i.e., x = \frac{\pi}{2} ). As \delta_{1} changes from 0 to 1, medical resources are allocated to nearby rural area (near x = 0 and x = \pi ), and the gap in medical conditions between rural and urban areas gradually narrows, and is eventually evenly distributed in space. Through calculation, it can be found that if \delta_{1} = 0 , in order to make R_{0} less than 1, a_{1}\geq0.0637 is needed, but if \delta_{1} = 1 , the medical resources only need 0.0576, which can make R_{0} < 1 . From Figure 5(c), we come to an interesting conclusion that in the case of limited medical resources, the spatial heterogeneity of the distribution of medical resources increases the risk of disease transmission. This result indicates that the distribution of medical resources dose play an important role in the design of malaria control programs.

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

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

    Figure 6.  The evolution of infection compartments of humans and mosquitoes with/without diffusion.

    As we all know, in the study of malaria transmission, the combined effects of temperature-dependent EIP and maturation delay, vector-bias, spatial structure and seasonal variation are worth studying. Here, we establish and analyze a reaction-diffusion model of malaria, including periodic delays and vector-bias in a heterogeneous environment.

    Using the theory proposed by Liang et al. [41] and Zhao [42], the basic reproduction number R_{0} can be introduced. We first prove that in a defined phase space, linear system (4.2) generates a periodic semiflow which is eventually strongly monotonic. Then, we describe R_{0} as a threshold for the extinction and persistence of malaria. It is proved that if R_{0} < 1 , then malaria will be eliminated. When R_{0} > 1 , there admits at least one positive \omega -periodic solution.

    In the numerical part, data for exploring the transmission of malaria in Maputo Provence, Mozambique have been published. The numerical simulation is divided into three parts. First, we give some simulation conclusions to support the theoretical analysis. Second, we discuss the sensitivity of R_{0} to the parameters in the model. (I) Our numerical results show that if we only consider the time-average maturation delay, we can underestimate the risk of malaria transmission assessed by R_{0} . (II) It is easy to know that R_{0} is decreasing with respect to l/p . Thence, if the impact of vector-bias effect on malaria transmission is ignored, the risk of transmission will be underestimated. (III) We analyze the impact of intervention measures on R_{0} from two aspects: (i) the impact of increasing medical resources. This case is mainly due to the uneven distribution of medical resources in rural and urban areas, which increases the impact of medical resources on R_{0} ; (ii) for fixed medical resources, the impact of different allocation measures on disease transmission. Through the analysis of this case, we find that in the case of limited medical resources, narrowing the gap in the allocation of medical resources between rural and urban areas can better reduce the risk of malaria transmission. Third, we explore the impact of diffusion on R_{0} in a spatially heterogeneous environment, assuming that the distribution of medical resources and the initial data of mosquitoes and humans depend on space. Biologically speaking, we understand this spatial heterogeneity as the difference between rural and urban areas. In addition, we also study the impact of diffusion on the final size of malaria. We find that under the set parameters, if the diffusion of mosquitoes and humans is not considered, the disease will become extinct in urban areas, and persist in rural areas. But under this set of parameters, the disease will become extinct if the diffusion is not considered.

    We are very grateful to anonymous referees and editor for careful reading and helpful suggestions which led to an improvement of our original manuscript. The authors are grateful for the support and useful comments provided by Professor Xiaoqiang Zhao (Memorial University of Newfoundland). This research is supposed by National Natural Science Foundation of China (grant number: 11971013) and Postgraduate Research & Practice Innovation Program of Jiangsu Province (KYCX 21\_0176 ) and Nanjing University of Aeronautics and Astronautics PhD short-term visiting scholar project (ZDGB2021027).

    The authors declare there is no conflict of interest.

    Proof. Firstly, the local existence of the unique mild solution is proved. Clearly, F is locally Lipschitz continuous. It is necessary to show

    \begin{equation} \begin{aligned} \lim\limits_{\hat{\theta}\rightarrow 0} {\rm{dist}}\left(\phi(0,\cdot)+\hat{\theta}F(t,\phi),\mathbb{X}_{1}^{+}\right) = 0, \forall (t,\phi)\in[0,+\infty)\times W_{h}. \end{aligned} \end{equation} (A.1)

    For all (t, \phi)\in[0, +\infty)\times W_{h} and \hat{\theta}\geq 0 , because p\geq l > 0 , one has

    \begin{array}{l} \phi(0,x)+\hat{\theta} F(t,\phi)(x)\\ = \left(\begin{array}{c} \phi_{1}(0,x)+\hat{\theta}\bigg( (1-\tau_{1}^{\prime}(t))\Lambda_{v}(t-\tau_{1}(t))e^{-\int_{t-\tau_{1}(t)}^{t}\mu_{m}(r,x) {\rm{d}}r} -\frac{\alpha \beta(t,x)p\phi_{3}(0,x)\phi_{1}(0,x)}{p\phi_{3}(0,x)+l[N(x)-\phi_{3}(0,x)]}\bigg)\\ \phi_{2}(0,x)+\hat{\theta}(1-\tau_{2}^{\prime}(t))\int_{\Omega}\Gamma(t,t-\tau_{2}(t),x,y) \frac{\alpha \beta(t-\tau_{2}(t),y)p\phi_{3}(-\tau_{2}(t),y)\phi_{1}(-\tau_{2}(t),y)}{p\phi_{3}(-\tau_{2}(t),y)+l[N(y)-\phi_{3}(-\tau_{2}(t),y)]} {\rm{d}}y\\ \phi_{3}(0,x)+\hat{\theta}\frac{c\beta(t,x)l[N(x)-\phi_{3}(0,x)]\phi_{2}(0,x)}{p\phi_{3}(0,x)+l[N(x)-\phi_{3}(0,x)]} \end{array}\right)\\ \geq\left(\begin{array}{c} \phi_{1}(0,x)+\hat{\theta}\bigg( (1-\tau_{1}^{\prime}(t))\Lambda_{v}(t-\tau_{1}(t))e^{-\int_{t-\tau_{1}(t)}^{t}\mu_{m}(r,x) {\rm{d}}r} -\frac{\alpha \beta(t,x)p\phi_{3}(0,x)\phi_{1}(0,x)}{lN(x)}\bigg) \\ \phi_{2}(0,x)+\hat{\theta}(1-\tau_{2}^{\prime}(t))\int_{\Omega}\Gamma(t,t-\tau_{2}(t),x,y) \frac{\alpha \beta(t-\tau_{2}(t),y)p\phi_{3}(-\tau_{2}(t),y)\phi_{1}(-\tau_{2}(t),y)}{lN(y)} {\rm{d}}y \\ \phi_{3}(0,x)+\hat{\theta}\frac{c\beta(t,x)l[N(x)-\phi_{3}(0,x)]\phi_{2}(0,x)}{pN(x)} \end{array}\right) \\ \geq\left(\begin{array}{c} \phi_{1}(0,x)\left(1-\hat{\theta}\frac{\alpha \hat{\beta}p\phi_{3}(0,x)}{l\tilde{N}} \right) \\ \phi_{2}(0,x) \\ \phi_{3}(0,x)\left(1-\hat{\theta}\frac{c\hat{\beta}l\phi_{2}(0,x)}{p\tilde{N}}\right) \end{array}\right), \end{array}

    for t\geq 0 and x\in \bar{\Omega} , where \hat{\beta} = \max\limits_{t\in[0, \omega], x\in\bar{\Omega}}\beta(t, x) and \tilde{N} = \min\limits_{x\in\bar{\Omega}}N(x) , and

    \begin{equation} \nonumber \begin{aligned} N(x)-[\phi_{3}(0,x)+\hat{\theta}F_{3}(t,\phi)(x)] = &[N(x)-\phi_{3}(0,x)] \left[1-\frac{\hat{\theta}c\beta(t,x)l\phi_{2}(0,x)}{p\phi_{3}(0,x)+l[N(x)-\phi_{3}(0,x)]}\right]\\ \geq& [N(x)-\phi_{3}(0,x)]\left[1-\frac{\hat{\theta}c\beta(t,x)\phi_{2}(0,x)}{N(x)}\right]. \end{aligned} \end{equation}

    This implies that

    \begin{equation} \nonumber \begin{aligned} \lim\limits_{\hat{\theta}\rightarrow 0} {\rm{dist}}\left(\phi(0,\cdot)+\hat{\theta}F(t,\phi),\mathbb{Y}^{+}\times\mathbb{Y}^{+}\times\mathbb{Y}_{h}\right) = 0, \quad \forall (t,\phi)\in[0,+\infty)\times W_{h}. \end{aligned} \end{equation}

    Consequently, by reference [36,Corollary 4] with K = \mathbb{Y}^{+}\times\mathbb{Y}^{+}\times\mathbb{Y}_{h} and S(t, s) = T(t, s) , system (2.12) has a unique mild solution u(t, \cdot, \phi) on its maximum existence interval t\in[0, t_{\phi}) and u_{0} = \phi where t_{\phi}\leq \infty , u(t, \cdot, \phi)\in \mathbb{Y}^{+}\times\mathbb{Y}^{+}\times\mathbb{Y}_{h} . In addition, when t > \hat{\tau} , u(t, \cdot, \phi) is a classical solution, based on the analyticity of T(t, s) , s, t\in \mathbb{R} and t > s .

    Analyze the following time-periodic reaction-diffusion equation

    \begin{equation} \left\{ \begin{aligned} \frac{\partial \hat{w}(t,x)}{\partial t} = &D_{v}\Delta \hat{w}(t,x)-\mu_{v}(t,x)\hat{w}(t,x)+(1-\tau^{\prime}_{1}(t))\Lambda_{v}(t-\tau_{1}(t),x) e^{-\int_{t-\tau_{1}(t)}^{t}\mu_{m}(r,x) {\rm{d}}r},&x\in\Omega, \\ \frac{\partial \hat{w}(t,x)}{\partial \nu} = &0, &x\in\partial\Omega. \end{aligned} \right. \end{equation} (A.2)

    By reference [38,Lemma 2.1], system (A.2) has a unique globally attractive positive \omega -periodic solution u_{1}^{*}(t, x) in \mathbb{Y}^{+} . Since system (A.2) can dominate the first equation of (2.12), there is a B_{1} > 0 such that for any \phi \in W_{h} , a positive integer l_{1} = l_{1}(\phi) > 0 exists and meets u_{1}(t, x, \phi)\leq B_{1} for all t\geq l_{1}\omega and x\in \bar{\Omega} .

    Next, similar to reference [46,Theorem 2.1], the ultimate boundedness of solutions can be proved. Set

    \begin{equation} \nonumber \begin{aligned} \tilde{\Lambda}_{v} = \max\limits_{t\in [0,\omega]}\int_{\Omega}\Lambda_{v}(t,x) {\rm{d}}x, \end{aligned} \end{equation}

    for any \phi\in W_{h} , let (u_{1}(t, x), u_{2}(t, x), u_{3}(t, x)): = (u_{1}(t, \phi)(x), u_{2}(t, \phi)(x), u_{3}(t, \phi)(x)) , t\geq 0 , x\in \bar{\Omega} , and \bar{u}_{i}(t) = \int_{\Omega}u_{i}(t, x) {\rm{d}}x where i = 1, 2, 3 . Using the Green's formula and integrating the first equation of system (2.12), one has

    \begin{equation} \nonumber \begin{aligned} \frac{ {\rm{d}} \bar{u}_{1}(t)}{ {\rm{d}}t} = & \int_{\Omega}(1-\tau_{1}^{\prime}(t)) \Lambda_{v}\left(t-\tau_{1}(t),x\right)e^{-\int_{t-\tau_{1}(t)}^{t}\mu_{m}(r,x) {\rm{d}}r} {\rm{d}}x \\ &- \int_{\Omega}\frac{\alpha\beta(t,x)pu_{3}u_{1}}{pu_{3} +l[N(x)-u_{3}]} {\rm{d}}x-\int_{\Omega}\mu_{v}(t,x)u_{1} {\rm{d}}x\\ \leq & \tilde{\Lambda}_{v}-\bar{\mu}_{v}\bar{u}_{1}(t)-\int_{\Omega}\frac{\alpha\beta(t,x)pu_{3}u_{1}}{pu_{3} +l[N(x)-u_{3}]} {\rm{d}}x, \end{aligned} \end{equation}

    that is,

    \begin{equation} \begin{aligned} \int_{\Omega}\frac{\alpha\beta(t,x)pu_{3}u_{1}}{pu_{3} +l[N(x)-u_{3}]} {\rm{d}}x\leq \tilde{\Lambda}_{v}-\bar{\mu}_{v}\bar{u}_{1}(t)-\frac{ {\rm{d}} \bar{u}_{1}(t)}{ {\rm{d}}t},\quad t > 0. \end{aligned} \end{equation} (A.3)

    By system (A.3), the property of the fundamental solutions [47] and Green's formula, the second equation of (2.12) are integrated to get

    \begin{equation} \begin{aligned} \frac{ {\rm{d}}\bar{u}_{2}(t)}{ {\rm{d}}t} = &-\int_{\Omega}\mu_{v}(t,x)u_{2} {\rm{d}}x+ (1-\tau_{2}^{\prime}(t))\int_{\Omega}\int_{\Omega}\Gamma(t,t-\tau_{2}(t),x,y)\cdot\\ &\frac{\alpha\beta(t-\tau_{2}(t),y)pu_{3}(t-\tau_{2}(t),y) u_{1}(t-\tau_{2}(t),y)}{pu_{3}(t-\tau_{2}(t),y) +l[N(y)-u_{3}(t-\tau_{2}(t),y)]} {\rm{d}}y {\rm{d}}x\\ \leq & -\bar{\mu}_{v}\bar{u}_{2}(t)-k_{1}\bar{u}_{2}(t-\tau_{2}(t))-k_{2}\frac{ {\rm{d}}\bar{u}_{1}(t-\tau_{2}(t))}{ {\rm{d}}t}+k_{3}, \quad\forall t\geq l_{1}\omega+\hat{\tau}, \end{aligned} \end{equation} (A.4)

    where k_{1} , k_{2} and k_{3} are positive numbers independent of \phi . Select k_{1}\leq \bar{\mu}_{v}k_{2} in system (15.4), so that

    \begin{equation} \nonumber \begin{aligned} \frac{ {\rm{d}}}{ {\rm{d}}t}[\bar{u}_{2}(t)+k_{2}\bar{u}_{1}(t-\tau_{2}(t))]&\leq -\bar{\mu}_{v}\bar{u}_{2}(t)-k_{1}\bar{u}_{1}(t-\tau_{2}(t))+k_{3}\leq -\frac{k_{1}}{k_{2}}\bar{u}_{2}(t)-k_{1}\bar{u}_{1}(t-\tau_{2}(t))+k_{3}\\ &\leq-\frac{k_{1}}{k_{2}}[\bar{u}_{2}(t)+k_{2}\bar{u}_{1}(t-\tau_{2}(t))]+k_{3},\quad\forall t\geq l_{1}\omega+\hat{\tau}, \end{aligned} \end{equation}

    which yields \bar{u}_{2}(t)+k_{2}\bar{u}_{1}(t-\tau_{2}(t))\leq \frac{k_{2}k_{3}}{k_{1}}+1 for t\geq l^{\prime}_{1}\omega+\hat{\tau} , where l^{\prime}_{1} > l_{1} is some integer. Hence,

    \begin{equation} \nonumber \begin{aligned} \bar{u}_{2}(t) = \|u_{2}(t,\cdot)\|_{L^{1}(\Omega)}\leq \frac{k_{2}k_{3}}{k_{1}}+1, \quad t\geq l^{\prime}_{1}\omega. \end{aligned} \end{equation}

    Integrating the third equation on model (2.12), and applying Green's formula, we have

    \begin{equation} \nonumber \begin{aligned} \frac{ {\rm{d}}\bar{u}_{3}(t)}{ {\rm{d}}t} = &\int_{\Omega}\frac{c\beta(t,x)l[N(x)-u_{3}(t,x)]u_{2}(t,x)}{pu_{3}(t,x)+l[N(x)-u_{3}(t,x)]} {\rm{d}}x -\int_{\Omega}(\mu_{h}(x)+\varrho(x))u_{3}(t,x) {\rm{d}}x \leq c\hat{\beta}\bar{u}_{2}(t)-(\bar{\mu}_{h}+\bar{\varrho})\bar{u}_{3}(t). \end{aligned} \end{equation}

    Thus,

    \begin{equation} \nonumber \begin{aligned} \bar{u}_{3}(t)\leq\frac{\left(\frac{k_{2}k_{3}}{k_{1}}+1\right)c\hat{\beta}}{\bar{\mu}_{h}+\bar{\varrho}}, \end{aligned} \end{equation}

    for t\geq l_{2}\omega+\hat{\tau} , (l_{2}\geq l_{1}) and x\in \bar{\Omega} . Therefore, we have t_{\phi} = \infty for each \phi \in W_{h} .

    Let \{Q_{t}\}_{t\geq 0} be a family of operators on W_{h} and Q_{t}(\phi)(s, x) = u_{t}(s, x, \phi) = u(t+s, x, \phi) for t\geq 0 , s\in [-\hat{\tau}, 0] , x\in \Omega and \phi \in W_{h} . Similar to the proof of Lemma 2.1 in reference [38], we can show that \{Q_{t}\}_{t\geq 0} is an \omega -periodic semiflow on W_{h} . From the above proofs, we conclude that Q_{t} : W_{h} \rightarrow W_{h} is point dissipative. Moreover Q: = Q_{\omega} is \kappa -contraction, where \kappa is Kuratowski measure of noncompactness, and hence, Q is asymptotically smooth. By references [48,Theorem 2.4.7] and [32,Theorem 1.1.2], it follows that Q has a global compact attractor.

    Proof. Let \bar{\tau}_{2} = \min\limits_{t\in[0, \omega]}\tau_{2}(t) . For all t\in [0, \bar{\tau}_{2}] , since t-\tau_{2}(t) strictly increase with respect to t , one has

    \begin{equation} \nonumber \begin{aligned} -\tau_{2}(0) = 0-\tau_{2}(0)\leq t-\tau_{2}(t)\leq \bar{\tau}_{2}-\tau_{2}(\bar{\tau}_{2})\leq \bar{\tau}_{2}-\bar{\tau}_{2} = 0, \end{aligned} \end{equation}

    and therefore, \tilde{w}_{2}(t-\tau_{2}(t), \cdot) = \hat{\psi}_{2}(t-\tau_{2}(t), \cdot) . Thus, for any t\in [0, \bar{\tau}_{2}] , there holds

    \begin{equation} \nonumber \left\{ \begin{aligned} \frac{\partial \tilde{w}_{1}}{\partial t} = &D_{v}\Delta \tilde{w}_{1}-\mu_{v}(t,x)\tilde{w}_{1}+(1-\tau^{\prime}_{2}(t)) \int_{\Omega} \Gamma(t,t-\tau_{2}(t),x,y)\cdot\\ &\frac{\alpha\beta(t-\tau_{2}(t),y)p\hat{\psi}_{2}(t-\tau_{2}(t),y)u_{1}^{*}(t-\tau_{2}(t),y)}{lN(y)} {\rm{d}}y,&x\in \Omega,\\ \frac{\partial \tilde{w}_{2}}{\partial t} = &D_{h}\Delta \tilde{w}_{2}+c\beta(t,x)\tilde{w}_{1}-(\mu_{h}(x)+\varrho(x))\tilde{w}_{2}, &x\in \Omega,\\ \frac{\partial \tilde{w}_{1}}{\partial \nu} = &\frac{\partial \tilde{w}_{2}}{\partial \nu} = 0, &x\in \partial\Omega. \end{aligned} \right. \end{equation}

    Fix \hat{\psi}\in \mathcal{E}_{1}^{+} , for t\in [0, \bar{\tau}_{2}] , the solution (\tilde{w}_{1}(t, \cdot), \tilde{w}_{1}(t, \cdot)) of the above linear system exists uniquely. To put this another way, one has z_{1}(\theta_{3}, \cdot) = \tilde{w}_{1}(\theta_{3}, \cdot) , where \theta_{3}\in [0, \bar{\tau}_{2}] and z_{2}(\theta_{4}, \cdot) = \tilde{w}_{2}(\theta_{4}, \cdot) , \theta_{4}\in [-\tau_{2}(0), \bar{\tau}_{2}] .

    By repeating the above step to [n\bar{\tau}_{2}, (n+1)\bar{\tau}_{2}] , we see that the solution with initial date \hat{\psi}\in \mathcal{E}_{1}^{+} exists uniquely for all t\geq 0 .

    Proof. Similar to the proof of Lemma 4.2, on each interval [n\bar{\tau}_{2}, (n+1)\bar{\tau}_{2}] , n\in \mathbb{N} , one easily obtain that \tilde{w}_{i}(t, \cdot) \geq 0 , ( i = 1, 2 ) for all t\geq 0 . Select a K > \max\{\hat{\mu}_{v}, \hat{\mu}_{h}+\hat{\varrho}\} , such that for every t\in \mathbb{R} , g_{1}(t, \cdot, \tilde{w}_{1}): = -\mu_{v}(t, \cdot)\tilde{w}_{1}+K\tilde{w}_{1} is strictly increasing with respect to \tilde{w}_{1} , and g_{2}(t, \cdot, \tilde{w}_{2}): = -(\mu_{h}(\cdot)+\varrho(\cdot))\tilde{w}_{2}+K\tilde{w}_{2} is increasing in \tilde{w}_{2} . Then, for t > 0 , both \tilde{w}_{1} and \tilde{w}_{2} meet the following system

    \begin{equation} \nonumber \left\{ \begin{aligned} \frac{\partial \tilde{w}_{1}}{\partial t} = &D_{v}\Delta \tilde{w}_{1}-K \tilde{w}_{1}+g_{1}(t,\cdot, \tilde{w}_{1})+(1-\tau^{\prime}_{2}(t)) \int_{\Omega} \Gamma(t,t-\tau_{2}(t),x,y)\cdot\\ &\frac{\alpha\beta(t-\tau_{2}(t),y)p \tilde{w}_{2}(t-\tau_{2}(t),y)u_{1}^{*}(t-\tau_{2}(t),y)}{lN(y)} {\rm{d}}y, &x\in \Omega,\\ \frac{\partial \tilde{w}_{2}}{\partial t} = &D_{h}\Delta \tilde{w}_{2}-K \tilde{w}_{2}+g_{2}(t,\cdot, \tilde{w}_{2})+c\beta(t,x) \tilde{w}_{1},&x\in \Omega,\\ \frac{\partial \tilde{w}_{1}}{\partial \nu} = &\frac{\partial \tilde{w}_{2}}{\partial \nu} = 0,&x\in\partial\Omega. \end{aligned} \right. \end{equation}

    Therefore, for a fixed \bar{\xi}\in \mathcal{E}_{1}^{+} , one has

    \begin{equation} \left\{ \begin{aligned} \tilde{w}_{1}(t,\cdot,\bar{\xi}) = &\tilde{T}_{1}(t,0)\bar{\xi}_{1}(0)+\int_{0}^{t}\tilde{T}_{1}(t,s) g_{1}(s,\cdot, \tilde{w}_{1}(s,\cdot)) {\rm{d}}s+\int_{0}^{t}\tilde{T}_{1}(t,s)(1-\tau^{\prime}_{2}(t)) \cdot\\ & \int_{\Omega}\Gamma(s,s-\tau_{2}(s),\cdot,y)\frac{\alpha\beta(s-\tau_{2}(s),y)p \tilde{w}_{2}(s-\tau_{2}(s),y)u_{1}^{*}(s-\tau_{2}(s),y)}{lN(y)} {\rm{d}}y {\rm{d}}s,\\ \tilde{w}_{2}(t,\cdot,\bar{\xi}) = &\tilde{T}_{2}(t,0)\bar{\xi}_{2}(0)+\int_{0}^{t}\tilde{T}_{2}(t,s)g_{2}(s,\cdot, \tilde{w}_{2}(s,\cdot)) {\rm{d}}s +\int_{0}^{t}\tilde{T}_{2}(t,s)c\beta(s,\cdot) \tilde{w}_{1}(s,\cdot) {\rm{d}}s, \end{aligned} \right. \end{equation} (C.1)

    where \tilde{T}_{1}(t, s), \tilde{T}_{2}(t, s): = \mathbb{Y}\rightarrow \mathbb{Y} are evolution operators related to \frac{\partial \tilde{w}_{1}}{\partial t} = D_{v}\Delta \tilde{w}_{1}-K \tilde{w}_{1} and \frac{\partial \tilde{w}_{2}}{\partial t} = D_{h}\Delta \tilde{w}_{2}-K \tilde{w}_{2} with Nuemann boundary condition. Since m(t): = t-\tau_{2}(t) are increasing in t\in \mathbb{R} , it easily follows that [-\tau_{2}(0), 0]\subseteq m([0, \hat{\tau}_{2}]) . Generally, we assume that \xi_{2} > 0 . Then there exists an (\theta, x_{0})\in [-\tau_{2}(0), 0]\times \Omega such that \tilde{w}_{2}(\theta, x_{0}) > 0 . According to the equation of (17.1), we have \tilde{w}_{1}(t, \cdot, \xi) > 0 for all t > \hat{\tau}_{2} . If s > 2\hat{\tau}_{2} , then s-\tau_{2}(s) > 2\hat{\tau}_{2}-\hat{\tau}_{2} = \hat{\tau}_{2} . It is easy to know that \forall t\geq 0 , \tilde{w}_{2}(t, \cdot, \xi) > 0 for all t > 2\hat{\tau}_{2} according to the second equation of (17.1). This means that there is \tilde{w}_{i}(t, \cdot)\geq 0 for all t > 2\hat{\tau}_{2} , i = 1, 2 , so the solution map \tilde{Q}_{t} is strongly positive for all t > 3\hat{\tau}_{2} .

    Proof. For every \hat{\psi}\in \mathcal{E}_{1}^{+} , let \tilde{w}(t, x, \hat{\psi}) = (\tilde{w}_{1}(t, x, \hat{\psi}), \tilde{w}_{2}(t, x, \hat{\psi})) be the solution of (4.2) with \tilde{w}_{0}(\hat{\psi}) = \hat{\psi} . Since \hat{\psi}\gg 0 , we get \tilde{w}_{t}(\hat{\psi})\gg 0 for all t\geq 0 will be a piece of cake. Denote

    \begin{equation} \nonumber \begin{aligned} w_{1}^{*}(t,x)& = e^{-\tilde{\mu}t}\tilde{w}_{1}(t,x,\hat{\psi}),\; \; \; t\geq 0,\; \; \; \qquad x\in \bar{\Omega}, \\ w_{2}^{*}(t,x)& = e^{-\tilde{\mu}t}\tilde{w}_{2}(t,x,\hat{\psi}),\; \; \; t\geq -\tau_{2}(0),\; \; x\in \bar{\Omega}. \end{aligned} \end{equation}

    Then w^{*}(t, x) = (w^{*}_{1}(t, x), w^{*}_{2}(t, x))\gg 0 for t\geq 0 , x\in\bar{\Omega} , and w^{*}(t, x) meets the linear-periodic system with \tilde{\mu}:

    \begin{equation} \left\{ \begin{aligned} \frac{\partial w_{1}^{*}(t,x)}{\partial t} = &D_{v}\Delta w_{1}^{*}(t,x)-(\mu_{v}(t,x)+\tilde{\mu})w_{1}^{*}(t,x)+(1-\tau^{\prime}_{1}(t)) \int_{\Omega}\Gamma\left(t,t-\tau_{2}(t),x,y\right)\cdot\\ &\frac{\alpha \beta(t-\tau_{2}(t),y)pu_{1}^{*}(t-\tau_{2}(t),y)e^{-\tilde{\mu}\tau_{2}(t)}w_{2}^{*}(t-\tau_{2}(t),y)}{lN(y)} {\rm{d}}y, &x\in \Omega,\\ \frac{\partial w_{2}^{*}(t,x)}{\partial t} = &D_{h}\Delta w_{2}^{*}(t,x)+c\beta(t,x)w_{1}^{*}(t,x)-(\mu_{h}(x)+\varrho(x)+\tilde{\mu})w_{2}^{*}(t,x), &x\in \Omega,\\ \frac{\partial w_{1}^{*}(t,x)}{\partial \nu} & = \frac{\partial w_{2}^{*}(t,x)}{\partial \nu} = 0, &x\in \partial\Omega, \end{aligned} \right. \end{equation} (D.1)

    for t > 0 . Hence, for all \theta\in[-\tau_{2}(0), 0] , x\in\bar{\Omega} , w^{*}(t, x) is the solution of the \omega -periodic system (D.1) with w_{0}^{*}(\theta, x) = (w^{*}_{1}(0, x), w^{*}_{2}(\theta, x)) = \left(\hat{\psi}_{1}(x), e^{-\tilde{\mu}\theta}\hat{\psi}_{2}(\theta, x)\right) , where w_{t}^{*}(\cdot, \cdot) = \left(w^{*}_{1}(t, \cdot), w^{*}_{2t}(\cdot, \cdot)\right) for each t\geq 0 with

    \begin{equation} \nonumber \begin{aligned} &w_{1}^{*}(t,x) = e^{-\tilde{\mu}t}\tilde{w}_{1}(t,x,\hat{\psi}),\qquad \qquad \qquad \qquad \qquad \quad\forall x\in\bar{\Omega},\\ &w^{*}_{2t}(\theta,x) = \tilde{w}_{2}^{*}(t+\theta,x) = e^{-\tilde{\mu}(t+\theta)}\tilde{w}_{2}(t+\theta,x,\hat{\psi}),\; \; \; \forall (\theta,x)\in [-\tau_{2}(0),0]\times\bar{\Omega}. \end{aligned} \end{equation}

    For every \theta\in[-\tau_{2}(0), 0] , x\in\bar{\Omega} , one has

    \begin{equation} \nonumber \begin{aligned} w_{1}^{*}(\omega,x) = &e^{-\tilde{\mu}\omega}(\tilde{Q}(\hat{\psi}))_{1}(\theta,x) = e^{-\tilde{\mu}\omega}r(\tilde{Q})\hat{\psi}_{1}(x) = \hat{\psi}_{1}(x) = w_{1}^{*}(0,x),\\ w_{2}^{*}(\omega+\theta,x) = &e^{-\tilde{\mu}(\omega+\theta)}(\tilde{Q}(\hat{\psi}))_{2}(\theta,x) = e^{-\tilde{\mu}(\omega+\theta)}r(\tilde{Q})\hat{\psi}_{2}(\theta,x) = e^{-\tilde{\mu}\theta}\hat{\psi}_{2}(\theta,x) = w_{2}^{*}(\theta,x). \end{aligned} \end{equation}

    Consequently, w^{*}_{0}(\theta, \cdot) = w^{*}_{\omega}(\theta, \cdot) for all \theta\in[-\tau_{2}(0), 0] , and by the existence and uniqueness of system (D.1), one has

    \begin{equation} \nonumber \begin{aligned} w_{1}^{*}(t,x)& = w_{1}^{*}(t+\omega,x),\; \; \; \forall t\geq 0,\; \; \; \qquad x\in \bar{\Omega},\\ w_{2}^{*}(t,x)& = w_{2}^{*}(t+\omega,x),\; \; \; \forall t\geq -\tau_{2}(0),\; \; x\in \bar{\Omega}. \end{aligned} \end{equation}

    Thus, w^{*}(t, x) is an \omega -periodic solution of (D.1) and e^{\tilde{\mu}t}w^{*}(t, x) is a solution of (4.2).

    The authors in reference [20] discussed the impact of seasonality on the spread of malaria including assessing the seasonal related bite rate \beta(t) , periodic maturity delay \tau_{1}(t) , periodic EIP \tau_{2}(t) , periodic death \mu_{v}(t) and recruitment rate \Lambda_{v}(t) of mosquitoes, where

    \begin{equation} \begin{aligned} \beta(t) = &6.983-0.0709\sin(\pi t/2)-1.993\cos(\pi t/6)-0.128\cos(\pi t/2) -0.007642\sin(\pi t/3)\\ &-0.4247\cos(\pi t/3)+0.05452\sin(2\pi t/3) -1.459\sin(\pi t/6)-0.04095\cos(2\pi t/3)\\ &+0.0005486\cos(5\pi t/6)-0.06235 \sin(5\pi t/3)\; \; {\rm{Month}}^{-1}, \end{aligned} \end{equation} (E.1)
    \begin{equation} \begin{aligned} \tau_{1}(t) = &0.5605-0.03366 \cos (\pi t/6)+0.01178\sin (\pi t/6)-0.0008135\cos(\pi t/2)-0.00266\cos(2\pi t/3) \\ &-0.0001281\cos(5\pi t/6)+0.005958\cos(\pi t/3)+0.01178\sin (\pi t/6)-0.001072\sin (2\pi t/3)\\ &-0.001936\sin (\pi t/2)-0.0001059\sin(\pi t/3)+0.005026\sin (5\pi t/6) \rm { Month}, \end{aligned} \end{equation} (E.2)
    \begin{equation} \begin{aligned} \tau_{2}(t) = &1/30.4(17.25+4.806\sin(\pi t/6)-0.3257\cos(\pi t)+8.369\cos(\pi t/6)+2.857\sin(\pi t/3) \\ &+1.197\cos(\pi t/2) +0.03578\cos(2\pi t/3)+0.6354\sin(5\pi t/6)+1.035\sin(2\pi t/3)\\ &-0.3505\cos(5\pi t/6)+1.963\sin(\pi t/2)+3.27\cos(\pi t/3)+0\sin(\pi t)) \; {\rm{Month}}, \end{aligned} \end{equation} (E.3)
    \begin{equation} \begin{aligned} \mu_{v}(t) = &3.086+0.01942\cos(\pi t/3)+0.0007665\cos(2\pi t/3)+0.007133\cos(\pi t/2)\\ &+0.04788\cos(\pi t/6)-0.001459\cos(5\pi t/6)+0.005687\sin(2\pi t/3)+0.01819\sin(\pi t/3)\\ &+0.01135\sin(\pi t/2)+0.02655\sin(\pi t/6)+0.003198 \sin(5\pi t/3)\; \; {\rm{Month}}^{-1}, \end{aligned} \end{equation} (E.4)
    \begin{equation} \begin{aligned} \Lambda_{v}(t) = \hat{k}\times \beta(t)\; ( {\rm{km}}^{2} {\rm{Month}})^{-1}, \; {\rm{where}} \; \hat{k} = 53.13\times 5. \end{aligned} \end{equation} (E.5)
    \begin{equation} \begin{aligned} \mu_{m}(t) = 9.0288\left(1+0.08\cos\left(\frac{\pi}{6}t\right)\right)\; \; {\rm{Month}}^{-1}. \end{aligned} \end{equation} (E.6)

    Other parameter values are defined as

    \begin{equation} \nonumber \begin{aligned} &D_{v} = 0.0125\; {\rm{km}}^{2} {\rm{Month}}^{-1}, \quad D_{h} = 0.1\; {\rm{km}}^{2} {\rm{Month}}^{-1},\quad \alpha = 0.64,\\ &c = 0.21,\quad \mu_{h} = 0.001574,\quad p = 0.8, \quad l = 0.6. \end{aligned} \end{equation}
    Table A1.  Parameter definitions and values.
    Parameter Definition Value(range) References
    N_{h} Total human population number 53 ({\rm{km}}^{2})^{-1} [24]
    \mu_{h} Natural mortality rate of humans 0.00157 {\rm{Month}}^{-1} [20]
    \varrho Constant recovery rate of humans (0.04256, 0.5168) {\rm{Month}}^{-1} [24]
    l/p Vector-bias parameter (0, 1) [13]
    c Transmission probability (0.01, 0.27) [50]
    from infectious mosquitoes to humans
    D_{h} Diffusion rate of humans 0.1 {\rm{km}}^{2} {\rm{Month}}^{-1} [24]
    D_{v} Diffusion rate of mosquitoes 0.0125 {\rm{km}}^{2} {\rm{Month}}^{-1} [24]
    \alpha Transmission probability from (0.072, 0.64) [50]
    infectious humans to mosquitoes per bite
    \beta Biting rate of mosquitoes (E.1) [20]
    \mu_{v} Death rate of mosquitoes (E.4) [20]
    \tau_{1} Maturation delay of (E.2) [51]
    mosquitoes in the aquatic stage
    \tau_{2} EIP of mosquitoes (E.3) [24]
    \Lambda_{v} Recruitment rate of mosquitoes (E.5) [24]
    \mu_{m} Natural death rate for the aquatic stage (E.6) [52]

     | Show Table
    DownLoad: CSV


    [1] V. A. Ambartsumian, On the fluctuation of the brightness of the milky way, Dokl. Akad. Nauk. USSR, 44 (1994), 223–226.
    [2] J. Patade, S. Bhalekar, On analytical solution of Ambartsumian equation, Natl. Acad. Sci. Lett., 40 (2017), 291–293. https://doi.org/10.1007/s40009-017-0565-2 doi: 10.1007/s40009-017-0565-2
    [3] H. O. Bakodah, A. Ebaid, Exact solution of Ambartsumian delay differential equation and comparison with Daftardar-Gejji and Jafari approximate method, Mathematics, 6 (2018), 331. https://doi.org/10.3390/math6120331 doi: 10.3390/math6120331
    [4] A. Ebaid, A. A. Enazi, B. Z. Albalawi, M. D. Aljoufi, Accurate approximate solution of Ambartsumian delay differential equation via decomposition method, Math. Comput. Appl., 24 (2019), 7. https://doi.org/10.3390/mca24010007 doi: 10.3390/mca24010007
    [5] E. A. Algehyne, E. R. E. Zahar, F. M. Alharbi, A. Ebaid, Development of analytical solution for a generalized Ambartsumian equation, AIMS Math., 5 (2019), 249–258. https://doi.org/10.3934/math.2020016 doi: 10.3934/math.2020016
    [6] D. Kumar, J. Singh, D. Baleanu, S. Rathore, Analysis of a fractional model of the Ambartsumian equation, Eur. Phys. J. Plus, 133 (2018), 259. https://doi.org/10.1140/epjp/i2018-12081-3 doi: 10.1140/epjp/i2018-12081-3
    [7] M. Sezera, A. A. Dascıoglu, A Taylor method for numerical solution of generalized pantograph equations with linear functional argument, J. Comput. Appl. Math., 200 (2007), 217–225. https://doi.org/10.1016/j.cam.2005.12.015 doi: 10.1016/j.cam.2005.12.015
    [8] S. Sedaghat, Y. Ordokhani, M. Dehghan, Numerical solution of the delay differential equations of pantograph type via Chebyshev polynomials, Commun. Nonlinear Sci., 17 (2012), 4815–4830. https://doi.org/10.1016/j.cnsns.2012.05.009 doi: 10.1016/j.cnsns.2012.05.009
    [9] E. Tohidi, A. H. Bhrawy, K. Erfani, A collocation method based on Bernoulli operational matrix for numerical solution of generalized pantograph equation, Appl. Math. Model., 37 (2013), 4283–4294. https://doi.org/10.1016/j.apm.2012.09.032 doi: 10.1016/j.apm.2012.09.032
    [10] S. Javadi, E. Babolian, Z. Taheri, Solving generalized pantograph equations by shifted orthonormal Bernstein polynomials, J. Comput. Appl. Math., 303 (2016), 1–14. https://doi.org/10.1016/j.cam.2016.02.025 doi: 10.1016/j.cam.2016.02.025
    [11] J. Shen, T. Tang, L. Wang, Spectral methods: Algorithms, analysis and applications, Springer, 41 2011. https://doi.org/10.1007/978-3-540-71041-7
    [12] S. S. E. Eldien, On solving systems of multi-pantograph equations via spectral tau method, Appl. Math. Comput., 321 (2018), 63–73. https://doi.org/10.1016/j.amc.2017.10.014 doi: 10.1016/j.amc.2017.10.014
    [13] H. Jafari, M. Mahmoudi, M. H. N. Skandari, A new numerical method to solve pantograph delay differential equations with convergence analysis, Adv. Differ. Equ., 2021,129. https://doi.org/10.1186/s13662-021-03293-0
    [14] O. R. Isik, T. Turkoglu, A rational approximate solution for generalized pantograph-delay differential equations, Math. Method. Appl. Sci., 39 (2016), 2011–2024. https://doi.org/10.1002/mma.3616 doi: 10.1002/mma.3616
    [15] A. Ebaid, A. M. Wazwaz, E. Alali, B. Masaedeh, Hypergeometric series solution to a class of second-order boundary value problems via Laplace transform with applications to nanofuids, Commun. Theor. Phys., 67 (2017), 231. https://doi.org/10.1088/0253-6102/67/3/231 doi: 10.1088/0253-6102/67/3/231
    [16] P. V. Pavani, U. L. Priya, B. A. Reddy, Solving differential equations by using Laplace transforms, Int. J. Res. Anal. Rev., 5 (2018), 1796–1799.
    [17] H. S. Ali, E. Alali, A. Ebaid, F. M. Alharbi, Analytic solution of a class of singular second-order boundary value problems with applications, Mathematics, 7 (2019), 172. https://doi.org/10.3390/math7020172 doi: 10.3390/math7020172
    [18] N. Dogan, Solution of the system of ordinary differential equations by combined Laplace transform-Adomian decomposition method, Math. Comput. Appl., 17 (2012), 203–211. https://doi.org/10.3390/mca17030203 doi: 10.3390/mca17030203
    [19] A. Atangana, B. S. Alkaltani, A novel double integral transform and its applications, J. Nonlinear Sci. Appl., 9 (2016), 424–434. https://doi.org/10.22436/jnsa.009.02.08 doi: 10.22436/jnsa.009.02.08
    [20] K. Abbaoui, Y. Cherruault, Convergence of Adomian's method applied to nonlinear equations, Math. Comput. Model., 20 (1994), 69–73. https://doi.org/10.1016/0895-7177(94)00163-4 doi: 10.1016/0895-7177(94)00163-4
    [21] Y. Cherruault, G. Adomian, Decomposition methods: A new proof of convergence, Math. Comput. Model., 18 (1993), 103–106. https://doi.org/10.1016/0895-7177(93)90233-O doi: 10.1016/0895-7177(93)90233-O
    [22] A. Alshaery, A. Ebaid, Accurate analytical periodic solution of the elliptical Kepler equation using the Adomian decomposition method, Acta Astronaut., 140 (2017), 27–33. https://doi.org/10.1016/j.actaastro.2017.07.034 doi: 10.1016/j.actaastro.2017.07.034
    [23] Z. Ayati, J. Biazar, On the convergence of homotopy perturbation method, J. Egypt. Math. Soc., 23 (2015), 424–428. https://doi.org/10.1016/j.joems.2014.06.015 doi: 10.1016/j.joems.2014.06.015
    [24] S. A. Pasha, Y. Nawaz, M. S. Arif, The modified homotopy perturbation method with an auxiliary term for the nonlinear oscillator with discontinuity, J. Low Freq. Noise V. A., 38 (2019), 1363–1373. https://doi.org/10.1177/0962144X18820454 doi: 10.1177/0962144X18820454
    [25] M. Bayat, I. Pakar, M. Bayat, Approximate analytical solution of nonlinear systems using homotopy perturbation method, P. I. Mech. Eng. E-J. Pro., 230 (2016), 10–17. https://doi.org/10.1177/0954408914533104 doi: 10.1177/0954408914533104
    [26] S. Ahmad, A. Ullah, A. Akgül, M. D. L. Sen, A novel homotopy perturbation method with applications to nonlinear fractional order KdV and Burger equation with exponential-decay kernel, J. Funct. Space., 2021, 8770488. https://doi.org/10.1155/2021/8770488
    [27] A. Arikoglu, I. Ozkol, Solution of fractional differential equations by using differential transform method, Chaos Soliton. Fract., 34 (2007), 1473–1481. https://doi.org/10.1016/j.chaos.2006.09.004 doi: 10.1016/j.chaos.2006.09.004
    [28] S. H. Chang, I. L. Chang, A new algorithm for calculating one-dimensional differential transform of nonlinear functions, Appl. Math. Comput., 195 (2008), 799–808. https://doi.org/10.1016/j.amc.2007.05.026 doi: 10.1016/j.amc.2007.05.026
    [29] S. J. Liao, The proposed homotopy analysis technique for the solution of nonlinear problems, Ph. D. Thesis, Shanghai Jiao Tong University, Shanghai, China, 1992.
    [30] S. J. Liao, Homotopy analysis method: A new analytical technique for nonlinear problems, Commun. Nonlinear Sci., 2 (1997), 95–100. https://doi.org/10.1016/S1007-5704(97)90047-2 doi: 10.1016/S1007-5704(97)90047-2
    [31] S. J. Liao, An explicit, totally analytic approximation of Blasius' viscous flow problems, Int. J. Nonlin. Mech., 34 (1999), 759–778. https://doi.org/10.1016/S0020-7462(98)00056-0 doi: 10.1016/S0020-7462(98)00056-0
    [32] S. J. Liao, Beyond perturbation: Introduction to the homotopy analysis method, Chapman & Hall/CRC Press: Boca Raton, FL, USA, 2003.
    [33] S. J. Liao, On the homotopy analysis method for nonlinear problems, Appl. Math. Comput., 147 (2004), 499–513. https://doi.org/10.1016/S0096-3003(02)00790-7 doi: 10.1016/S0096-3003(02)00790-7
    [34] O. Nave, Modification of semi-analytical method applied system of ODE, Modern Appl. Sci. CCSE, 14 (2020), 75–81. https://doi.org/10.5539/mas.v14n6p75 doi: 10.5539/mas.v14n6p75
    [35] A. H. S. A. Enazy, A. Ebaid, E. A. Algehyne, H. K. A. Jeaid, Advanced study on the delay differential equation y'(t) = ay(t)+by(ct), Mathematics, 10 (2022), 4302. https://doi.org/10.3390/math10224302 doi: 10.3390/math10224302
    [36] S. M. Khaled, Applications of standard methods for solving the electric train mathematical model with proportional delay, Int. J. Anal. Appl., 20 (2022), 27. https://doi.org/10.28924/2291-8639-20-2022-27 doi: 10.28924/2291-8639-20-2022-27
    [37] A. B. Albidah, N. E. Kanaan, A. Ebaid, H. K. A. Jeaid, Exact and numerical analysis of the pantograph delay differential equation via the homotopy perturbation method, Mathematics, 11 (2023), 944. https://doi.org/10.3390/math11040944 doi: 10.3390/math11040944
    [38] E. R. E. Zahar, A. Ebaid, Analytical and numerical simulations of a delay model: The pantograph delay equation, Axioms, 11 (2022), 741. https://doi.org/10.3390/axioms11120741 doi: 10.3390/axioms11120741
    [39] R. Alrebdi, H. K. A. Jeaid, Accurate solution for the pantograph delay differential equation via Laplace transform, Mathematics, 11 (2023), 2031. https://doi.org/10.3390/math11092031 doi: 10.3390/math11092031
    [40] Y. Plusquellec, J. L. Steimer, An analytical solution for a class of delay-differential equations in pharmacokinetic compartment modeling, Math. Biosci., 70 (1984), 39–56. https://doi.org/10.1016/0025-5564(84)90045-2. doi: 10.1016/0025-5564(84)90045-2
    [41] F. V. Difonzo, P. Przybyłowicz, Y. Wu, Existence, uniqueness and approximation of solutions to Carathéodory delay differential equations, J. Comput. Appl. Math., 436 (2024), 115411. https://doi.org/10.1016/j.cam.2023.115411 doi: 10.1016/j.cam.2023.115411
  • This article has been cited by:

    1. Hocine Makheloufi, Stability Analysis of Viscoelastic Swelling Porous Elastic Soils With Nonlinear Kelvin–Voigt Dampings, 2025, 0170-4214, 10.1002/mma.10889
  • Reader Comments
  • © 2024 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(423) PDF downloads(21) Cited by(0)

Figures and Tables

Figures(16)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog