Processing math: 40%
Research article Special Issues

Transition graph decomposition for complex balanced reaction networks with non-mass-action kinetics


  • Reaction networks are widely used models to describe biochemical processes. Stochastic fluctuations in the counts of biological macromolecules have amplified consequences due to their small population sizes. This makes it necessary to favor stochastic, discrete population, continuous time models. The stationary distributions provide snapshots of the model behavior at the stationary regime, and as such finding their expression in terms of the model parameters is of great interest. The aim of the present paper is to describe when the stationary distributions of the original model, whose state space is potentially infinite, coincide exactly with the stationary distributions of the process truncated to finite subsets of states, up to a normalizing constant. The finite subsets of states we identify are called copies and are inspired by the modular topology of reaction network models. With such a choice we prove a novel graphical characterization of the concept of complex balancing for stochastic models of reaction networks. The results of the paper hold for the commonly used mass-action kinetics but are not restricted to it, and are in fact stated for more general setting.

    Citation: Daniele Cappelletti, Badal Joshi. Transition graph decomposition for complex balanced reaction networks with non-mass-action kinetics[J]. Mathematical Biosciences and Engineering, 2022, 19(8): 7649-7668. doi: 10.3934/mbe.2022359

    Related Papers:

    [1] Xia Wang, Yuming Chen . An age-structured vector-borne disease model with horizontal transmission in the host. Mathematical Biosciences and Engineering, 2018, 15(5): 1099-1116. doi: 10.3934/mbe.2018049
    [2] Maria Guadalupe Vazquez-Peña, Cruz Vargas-De-León, Jorge Velázquez-Castro . Global stability for a mosquito-borne disease model with continuous-time age structure in the susceptible and relapsed host classes. Mathematical Biosciences and Engineering, 2024, 21(11): 7582-7600. doi: 10.3934/mbe.2024333
    [3] Yanxia Dang, Zhipeng Qiu, Xuezhi Li . Competitive exclusion in an infection-age structured vector-host epidemic model. Mathematical Biosciences and Engineering, 2017, 14(4): 901-931. doi: 10.3934/mbe.2017048
    [4] Rundong Zhao, Qiming Liu, Huazong Zhang . Dynamical behaviors of a vector-borne diseases model with two time delays on bipartite networks. Mathematical Biosciences and Engineering, 2021, 18(4): 3073-3091. doi: 10.3934/mbe.2021154
    [5] Mostafa Adimy, Abdennasser Chekroun, Claudia Pio Ferreira . Global dynamics of a differential-difference system: a case of Kermack-McKendrick SIR model with age-structured protection phase. Mathematical Biosciences and Engineering, 2020, 17(2): 1329-1354. doi: 10.3934/mbe.2020067
    [6] Yicang Zhou, Zhien Ma . Global stability of a class of discrete age-structured SIS models with immigration. Mathematical Biosciences and Engineering, 2009, 6(2): 409-425. doi: 10.3934/mbe.2009.6.409
    [7] Rocio Caja Rivera, Shakir Bilal, Edwin Michael . The relation between host competence and vector-feeding preference in a multi-host model: Chagas and Cutaneous Leishmaniasis. Mathematical Biosciences and Engineering, 2020, 17(5): 5561-5583. doi: 10.3934/mbe.2020299
    [8] Jianxin Yang, Zhipeng Qiu, Xue-Zhi Li . Global stability of an age-structured cholera model. Mathematical Biosciences and Engineering, 2014, 11(3): 641-665. doi: 10.3934/mbe.2014.11.641
    [9] Jinliang Wang, Ran Zhang, Toshikazu Kuniya . A note on dynamics of an age-of-infection cholera model. Mathematical Biosciences and Engineering, 2016, 13(1): 227-247. doi: 10.3934/mbe.2016.13.227
    [10] Qiuyi Su, Jianhong Wu . Impact of variability of reproductive ageing and rate on childhood infectious disease prevention and control: insights from stage-structured population models. Mathematical Biosciences and Engineering, 2020, 17(6): 7671-7691. doi: 10.3934/mbe.2020390
  • Reaction networks are widely used models to describe biochemical processes. Stochastic fluctuations in the counts of biological macromolecules have amplified consequences due to their small population sizes. This makes it necessary to favor stochastic, discrete population, continuous time models. The stationary distributions provide snapshots of the model behavior at the stationary regime, and as such finding their expression in terms of the model parameters is of great interest. The aim of the present paper is to describe when the stationary distributions of the original model, whose state space is potentially infinite, coincide exactly with the stationary distributions of the process truncated to finite subsets of states, up to a normalizing constant. The finite subsets of states we identify are called copies and are inspired by the modular topology of reaction network models. With such a choice we prove a novel graphical characterization of the concept of complex balancing for stochastic models of reaction networks. The results of the paper hold for the commonly used mass-action kinetics but are not restricted to it, and are in fact stated for more general setting.



    1. Introduction

    Vector-borne diseases are infectious diseases caused by pathogens and parasites in human populations that are transmitted to people by blood-sucking arthropods, such as mosquitoes, ticks and fleas. They include some of the world's most destructive diseases, for instance, malaria, schistosomiasis, plague, and dengue fever. According to WHO [1], vector-borne diseases account for more than 17 % of all infectious diseases, causing more than 1 million deaths annually. In the past two decades, some vector-borne diseases, such as malaria and schistosomiasis, have continued to threaten human health. Furthermore, other vector-borne diseases have reemerged and broken out in different parts of the world, such as the 2014 Guangzhou outbreak of dengue fever and the outbreak of West Nile virus in North America since 1999. Any outbreak of the vector-borne diseases causes great harm to public health. As far as the 2014 Guangzhou outbreak of dengue fever is concerned, the total number of dengue fever cases reached 36,889 as of October 21st, 2014 [2], according to the provincial health and family planning commission. Due to the great harm to the public health caused by the vector diseases, it is imperative to understand the transmission dynamics of the vector-borne diseases firstly, and then discuss strategies to prevent and contain their outbreaks.

    Mathematical modeling has contributed significantly to our understanding of the epidemiology of infectious diseases [3,5]. Over the past two decades, there have been many published mathematical models focused on understanding the transmission dynamics of the vector-borne diseases ([4,12,23,26,27,30] and references therein). These models provided useful insights into the transmission dynamics of the vector-borne diseases. Almost all of the above models are described by ordinary differential equations (ODEs); therefore, some of the assumptions implicitly made in the formulation of these models [28] include: (1) infectious individuals are equally infectious during their infectious period; (2) the stage durations of the latent and infectious periods are exponentially distributed. Although in many cases these simplifying assumptions may provide a reasonable approximation to the biological process being modeled, it is important to examine how the model results may be influenced by these assumptions, which calls for an investigation of models that use more realistic assumptions [28].

    In this paper, we develop an age-structured model to study how transmission dynamics of the vector-borne diseases are affected by the incubation and infectious ages. The model studied in the paper incorporates both incubation age of the exposed hosts and infection age of the infectious hosts. Incubation age of the exposed hosts describes the different removal rates in the latent period, and infection age of the infectious hosts describes the variable infectiousness in the infectious period. Several recent studies [0,16,19,20,24,25,29] on age structured models have shown that age of infection may play an important role in the transmission dynamics of infectious diseases. Thieme and Castillo-Chavez [29] studied the effect of infection-age-independent infectivity on dynamics of HIV transmission, and showed that undamped oscillations may occur in particular if the variable infectivity is highly concentrated at certain parts of the incubation period. Lloyd [19,20] studied the epidemic model with the inclusion of non-exponential distributions of infectious periods. The results indicated that the inclusion of more realistic description of the recovery process may cause a significant destabilization of the model, and less dispersed distributions are seen to have two important epidemiological consequences: (1) less stable behavior is seen within the model; (2) disease persistence is diminished.

    Epidemic models with age of infection are usually described by first order partial differential equations, whose complexity makes them more difficult to theoretically analyze, particularly, their global behavior. Most existing studies on age-structured models focus only on the existence of non-trivial steady states [17,11] or give local stability results [32]. The stability analysis of nonlinear dynamical systems has always been a topic of both theoretical and practical importance since global stability is one of the most important issues related to their dynamic behaviors. However, proving the global stability is a very challenging task, especially for nonlinear systems described by PDEs due to the lack of generically applicable tools. The global stability results for the age-structured epidemic models were first obtained in [7,8,9]. The method of Lyapunov functions is the most common tool used to prove the global stability, especially for ODE models [14,15,18]. In recent years, Lyapunov function has been also used to study the global stability of epidemic models with age of infection [21,22,31].

    In this paper, we also use Lyapunov functions to study the global dynamics of a vector-borne disease model with incubation age of the exposed hosts and infection age of the infectious hosts. By using a class of Lyapunov functions we show that the global dynamics of the system is completely determined by the basic reproduction number R0: if R0<1 the disease-free equilibrium is globally asymptotically stable; if R0>1, a unique endemic equilibrium is globally asymptotically stable.

    This paper is organized as follows. In the next section we formulate a vector-borne epidemic model with incubation age of exposed hosts and infection age of infectious hosts. The two infection ages describe the different removal rates in the latent stage and the variable infectiousness in the infectious stage, respectively. We obtain an explicit formula for the basic reproduction number of system. Then we discuss the trivial and non-trivial equilibria and their stabilities. In Section 3, the global stability of the infection-free equilibrium of the system is analyzed by constructing a Lyapunov function. In Section 4, we show uniform strong persistence of the vector-borne disease if R0>1. In Section 5, we again use a Lyapunov function to derive the global stability of the epidemic equilibrium. Finally, a brief discussion is given in Section 6.


    2. The vector-borne disease model with two ages of infection and the local stabilities

    To introduce the model, we divide the host population under consideration into four groups: susceptible hosts at time t, denoted by Sh(t), infected but not infectious individuals Eh(τ,t), infected and infectious individuals Ih(a,t), and the number of recovered or immune individuals, denoted by Rh(t). The vector population, on the other hand, is divided into three compartmental classes: susceptible vector at time t, denoted by Sv(t), the number of recovered or immune vectors, denoted by Rv(t), and infected/infectious vectors Iv(t).

    With the above notation, we study the following infection-age-structured vector-borne epidemic model:

    {Sv(t)=ΛvSv(t)0βv(a)Ih(a,t)daμvSv(t),Iv(t)=Sv(t)0βv(a)Ih(a,t)da(μv+αv)Iv(t),Rv(t)=αvIv(t)μvRv(t),Sh(t)=ΛhβhSh(t)Iv(t)μhSh(t),Eh(τ,t)τ+Eh(τ,t)t=(μh+m(τ))Eh(τ,t),Eh(0,t)=βhSh(t)Iv(t),Ih(a,t)a+Ih(a,t)t=(μh+αh(a)+rh(a))Ih(a,t),Ih(0,t)=0m(τ)Eh(τ,t)dτ,Rh(t)=0rh(a)Ih(a,t)daμhRh(t). (1)

    In equation (1), Λh is the birth /recruitment rate of the host population. Let μv,μh be the natural death rate of the vectors and the host, respectively. m(τ) denotes the removal rate of the infected hosts with age of incubation τ from the latent period. αh(a) gives the additional disease induced death rate due to vector-borne disease at age of infection a. αv denotes the recovery rate of the infected vectors. rh(a) denotes the recovery rate of the infected hosts with age of infection a. Furthermore, βv(a) is the transmission coefficient of the infected host individuals at age of infection a, and βh is the transmission coefficient from infected vectors to healthy host individuals.

    To understand the model, notice that susceptible host individuals are recruited at a rate Λh. Susceptible host individuals can become infected by a bite of an infected mosquito with disease. Upon infection through biting by infected mosquitoes, the newly infected individuals move to the latent class, then progress into the infectious class with the progression rate m(τ). The non-infectious and infectious individuals infected by disease with age-since-infection equal to zero move to the boundary condition. The number total recovery rate from the infected class Ih(a,t) is given by the integral over all ages-since-infection. The susceptible vectors are recruited at a rate Λv. Susceptible mosquitos can become infected through biting on an infected individual of any age-since-infection at a specific age-infection transmission rate. As a consequence, the force of infection of susceptible vectors is given by the integral over all ages-since-infection. The total recovery rate from the infected vector class Iv(t) is given by αvIv(t).

    We notice that the equations for the recovered individuals and the recovered vectors are decoupled from the system and the analysis of system (1) is equivalent to the analysis of the system

    {Sv(t)=ΛvSv(t)0βv(a)Ih(a,t)daμvSv(t),Iv(t)=Sv(t)0βv(a)Ih(a,t)da(μv+αv)Iv(t),Sh(t)=ΛhβhSh(t)Iv(t)μhSh(t),Eh(τ,t)τ+Eh(τ,t)t=(μh+m(τ))Eh(τ,t),Eh(0,t)=βhSh(t)Iv(t),Ih(a,t)a+Ih(a,t)t=(μh+αh(a)+rh(a))Ih(a,t),Ih(0,t)=0m(τ)Eh(τ,t)dτ. (2)

    Model (2) is equipped with the following initial conditions:

    Sv(0)=Sv0,Iv(0)=Iv0,Sh(0)=Sh0,Eh(τ,0)=φ(τ),Ih(a,0)=ψ(a).

    All parameters are nonnegative, Λv>0, Λh>0, and βh>0, μv>0, μh>0. We make the following assumptions on the parameter-functions.

    Assumption 2.1 The parameter-functions satisfy the following.

    1. The functions βv(a) is bounded and uniformly continuous. When βv(a) is of compact support, the support has non-zero Lebesgue measure;

    2. The functions m(τ), αh(a), rh(a) belong to L(0,);

    3. The functions φ(τ), ψ(a) are integrable.

    Define the space of functions

    X=R×R×R×(L1(0,))×(L1(0,)).

    It can be verified that solutions of (2) with nonnegative initial conditions belong to the positive cone for t0. Furthermore, adding the first and the second equations we have

    ddt(Sv(t)+Iv(t))Λvμv(Sv(t)+Iv(t)).

    Hence,

    lim supt(Sv(t)+Iv(t))Λvμv.

    The number of the hosts can be bounded as follows:

    ddt(Sh(t)+0Eh(τ,t)dτ+0Ih(a,t)da)Λhμh(Sh(t)+0Eh(τ,t)dτ+0Ih(a,t)da).

    Hence,

    lim supt(Sh(t)+0Eh(τ,t)dτ+0Ih(a,t)da)Λhμh.

    Therefore, the following set is positively invariant for system

    Ω={(Sv,Iv,Sh,Eh,Ih)X+|(Sv(t)+Iv(t))Λvμv,(Sh(t)+0Eh(τ,t)dτ+0Ih(a,t)da)Λhμh}.

    Finally, since the exit rate of exposed host individuals from the incubation compartment is given by μh+m(τ), then the probability of still being latent after τ time units is given by

    π1(τ)=eμhτeτ0m(σ)dσ. (3)

    The exit rate of infected individuals from the infective compartment is given by μh+αh(a)+rh(a), thus the probability of still being infectious after a time units is given by

    π2(a)=eμhaea0(αh(σ)+rh(σ))dσ. (4)

    The reproduction number of disease in system (2) is given by the following expression

    R0=βhΛvΛhμvμh(μv+αv)0m(τ)π1(τ)dτ0βv(a)π2(a)da. (5)

    The reproduction number of disease gives the number of secondary infections produced in an entirely susceptible population by a typical infected individual during its entire infectious period. R0 gives the strength of vector-borne disease to invade when rare and alone. In particular, we notice that the reproduction number for vector-borne diseases is a product of the reproduction numbers of the two transmission processes: human-to-vector Rh and vector-to-human Rv,

    Rh=Λvμv0βv(a)π2(a)da,Rv=βhΛhμh(μv+αv)0m(τ)π1(τ)dτ,

    that is R0=RvRh. In the next section we compute explicit expressions for the equilibria and establish their local stability.

    System (2) always has a unique disease-free equilibrium E0, which is given by

    E0=(Sv0, 0, Sh0, 0, 0),

    where

    Sv0=Λvμv,Sh0=Λhμh.

    In addition, for Dengue virus there is a corresponding endemic equilibrium E1 given by

    E1=(Sv, Iv, Sh, Eh(τ), Ih(a)).

    We denote by

    Λ=βhΛhΛvμhμv(μv+αv),b=0m(τ)π1(τ)dτ0βv(a)π2(a)da,b(λ)=0m(τ)eλτπ1(τ)dτ0βv(a)eλaπ2(a)da. (6)

    The non-zero components of the equilibrium E1 are given by

    Iv=μvμh(R01)βh(Λhb+μv),Sv=Λv(μv+αv)Ivμv,Sh=ΛhβhIv+μh,Eh(τ)=Eh(0)π1(τ),Eh(0)=βhShIv,Ih(a)=Ih(0)π2(a),Ih(0)=Eh(0)0m(τ)π1(τ)dτ. (7)

    Next, we turn to the linearized equations for the disease-free equilibrium. To introduce the linearization at the disease-free equilibrium E0, we let Sv(t)=Sv0+xv(t), Iv(t)=yv(t), Sh(t)=Sh0+xh(t), Eh(τ,t)=zh(τ,t), Ih(a,t)=yh(a,t). The linearized system becomes

    {dxv(t)dt=Sv00βv(a)yh(a,t)daμvxv(t),dyv(t)dt=Sv00βv(a)yh(a,t)da(μv+αv)yv(t),dxh(t)dt=βhSh0yv(t)μhxh(t),zh(τ,t)τ+zh(τ,t)t=(μh+m(τ))zh(τ,t),zh(0,t)=βhSh0yv(t),yh(a,t)a+yh(a,t)t=(μh+αh(a)+rh(a))yh(a,t),yh(0,t)=0m(τ)zh(τ,t)dτ. (8)

    To study system (2), we look for solutions of the form xv(t)=ˉxveλt, yv(t)=ˉyveλt, xh(t)=ˉxheλt, zh(τ,t)=ˉzh(τ)eλt and yh(a,t)=ˉyh(a)eλt. We obtain the following eigenvalue problem

    {λˉxv=Sv00βv(a)ˉyh(a)daμvˉxv,λˉyv=Sv00βv(a)ˉyh(a)da(μv+αv)ˉyv,λˉxh=βhSh0ˉyvμhˉxh,dˉzh(τ)dτ=(λ+μh+m(τ))ˉzh(τ),ˉzh(0)=βhSh0ˉyv,dˉyh(a)da=(λ+μh+αh(a)+rh(a))ˉyh(a),ˉyh(0)=0m(τ)ˉzh(τ)dτ. (9)

    We notice that the two equations for ˉxv and ˉxh are decoupled from the equation for ˉyv, ˉzh, ˉyh. Hence, the equations for ˉxv and ˉxh are independent from the equations for ˉyv, ˉzh, ˉyh. Solving the differential equations for ˉzh, ˉyh, we have

    ˉzh(τ)=ˉzh(0) eλτπ1(τ)=βhSh0ˉyv eλτπ1(τ),ˉyh(a)=ˉyh(0) eλaπ2(a)=βhSh0ˉyv eλaπ2(a)0m(τ) eλτπ1(τ)dτ.  (10)

    Substituting for ˉyh(a) in the second equation of (9), we can obtain the following equation

    λ+μv+αv=βhSv0Sh00m(τ)eλτπ1(τ)dτ0βv(a)eλaπ2(a)da. (11)

    Now we are ready to establish the following result.

    Proposition 1. If

    R0<1,

    then the disease-free equilibrium is locally asymptotically stable. If R0>1, it is unstable.

    Proof. Assume

    R0<1.

    We set

    LHSdef=λ+μv+αv,RHSdef=G1(λ)=βhSv0Sh00m(τ)eλτπ1(τ)dτ0βv(a)eλaπ2(a)da. (12)

    Consider λ with λ0. For such λ, following from (12), we have that

    |LHS|μv+αv,|RHS|G1(λ)G1(0)=βhSv0Sh00m(τ)π1(τ)dτ0βv(a)π2(a)da=βhΛvΛhμvμh0m(τ)π1(τ)dτ0βv(a)π2(a)da=R0(μv+αv)<|LHS|.

    This gives a contradiction. Hence, we have shown that equation (11) cannot have any roots with non-negative real parts. Therefore, the disease-free equilibrium E0 depends on the eigenvalues of the equations for xv and xh. It is evident that λ=μv and λ=μh, so the disease-free equilibrium E0 is locally asymptotically stable if R0<1.

    Now assume

    R0>1.

    We rewrite the characteristic equation (11) in the form

    (λ+μv+αv)βhSv0Sh00m(τ)eλτπ1(τ)dτ0βv(a)eλaπ2(a)da=0. (13)

    We denote

    G2(λ)=(λ+μv+αv)βhSv0Sh00m(τ)eλτπ1(τ)dτ0βv(a)eλaπ2(a)da. (14)

    Thus equation (13) has turned into the following characteristic equation

    G2(λ)=0. (15)

    For λ real we have

    G2(0)=(μv+αv)βhSv0Sh00m(τ)π1(τ)dτ0βv(a)π2(a)da=(μv+αv)(1R0)<0.

    Furthermore, limλG2(λ)=+. Hence, the characteristic equation (15) has a real positive root. Therefore, the endemic equilibrium E0 is unstable. This concludes the proof.

    Now we turn to the local stability of the endemic equilibrium E1 if R0>1. The result on local stability of the equilibrium E1 is summarized below

    Proposition 2. Assume R0>1, then the endemic equilibrium E1 is locally asymptotically stable.

    Proof. We study the linearized equation around the endemic equilibrium E1. We introduce the following notation for the perturbations Sv(t)=Sv+xv(t), Iv(t)=Iv+yv(t),  Sh(t)=Sh+xh(t), Eh(τ,t)=Eh(τ)+zh(τ,t), Ih(a,t)=Ih(a)+yh(a,t). The system for the perturbations becomes (16)

    {dxv(t)dt=Sv0βv(a)yh(a,t)daxv(t)0βv(a)Ih(a)daμvxv(t),dyv(t)dt=Sv0βv(a)yh(a,t)da+xv(t)0βv(a)Ih(a)da(μv+αv)yv(t),dxh(t)dt=βhShyv(t)βhxh(t)Ivμhxh(t),dzh(τ)dτ=(λ+μh+m(τ))zh(τ,t),zh(0,t)=βhShyv(t)+βhxh(t)Iv,dyh(a)da=(λ+μh+αh(a)+rh(a))yh(a,t),yh(0,t)=0m(τ)zh(τ,t)dτ. (16)

    An approach similar to [8] (see Appendix B in [8]) can show that the linear stability of the system is in fact determined by the eigenvalues of the linearized system (16). To investigate the point spectrum, we look for exponential solutions (see the case of the disease-free equilibrium) and obtain a linear eigenvalue problem.

    {λxv=Sv0βv(a)yh(a)daxv0βv(a)Ih(a)daμvxv,λyv=Sv0βv(a)yh(a)da+xv0βv(a)Ih(a)da(μv+αv)yv,λxh=zh(0)μhxh,dzh(τ)dτ=(λ+μh+m(τ))zh(τ),zh(0)=βhShyv+βhIvxh,dyh(a)da=(λ+μh+αh(a)+rh(a))yh(a),yh(0)=0m(τ)zh(τ)dτ. (17)

    Solving the differential equation, we have

    zh(τ)=zh(0) eλτπ1(τ),yh(a)=yh(0) eλaπ2(a)=zh(0) eλaπ2(a)0m(τ) eλτπ1(τ)dτ.

    Substituting for yh in the second equation of (17), we can obtain the following equation

    {(λ+μv+0βv(a)Ih(a)da)xv+Svb(λ)zh(0)=0,xv0βv(a)Ih(a)da+(λ+μv+αv)yvSvb(λ)zh(0)=0,(λ+μh)xh+zh(0)=0,βhIvxhβhShyv+zh(0)=0. (18)

    By direct calculation, we obtain the following characteristic equation:

    (λ+μv+0βv(a)Ih(a)da)(λ+μv+αv)(λ+μh+βhIv)=βhShSvb(λ)(λ+μv)(λ+μh). (19)

    We divide both sides by (λ+μv)(λ+μh), then we introduce the following notation.

    G3(λ)=(λ+μv+0βv(a)Ih(a)da)(λ+μv+αv)(λ+μh+βhIv)(λ+μv)(λ+μh),G4(λ)=βhShSvb(λ)=βhShSv0m(τ)eλτπ1(τ)dτ0βv(a)eλaπ2(a)da. (20)

    Thus (19) can be expressed as the the equation

    G3(λ)=G4(λ). (21)

    If λ is a root with λ0, it follows from equation (20) that

    |G3(λ)|>|λ+μv+αv|μv+αv. (22)

    From system (2), we have

    βhSvSh0m(τ)π1(τ)dτ0βv(a)π2(a)da=μv+αv.

    Hence,

    |G4(λ)||G4(λ)|G4(0)=βhSvSh0m(τ)π1(τ)dτ0βv(a)π2(a)da=μv+αv<|G3(λ)|. (23)

    This leads to contradiction. Hence, for λ0, (21) has no solutions. Thus, the characteristic equation (19) has only solutions with negative real parts. Therefore, the endemic equilibrium E1 is locally asymptotically stable if R0>1. This concludes the proof.


    3. Global stability of the disease-free equilibrium

    In the previous section, we have established that equilibria are locally stable, that is, given the conditions on the parameters, if the initial conditions are close enough to the equilibrium, the solution will converge to that equilibrium. In this section our objective is to extend these results to global results. That is, given the conditions on the parameters, convergence to the equilibrium occurs independently of the initial conditions.

    As a first step, we establish the global stability of the disease-free equilibrium. We will use a Lyapunov function to approach the problem. We need to integrate the differential equation along the characteristic lines. Denote the initial condition by BE(t), BI(t):

    BE(t)=Eh(0,t),BI(t)=Ih(0,t).

    Integrating along the characteristic lines, we obtain

    Eh(τ,t)={BE(tτ)π1(τ), t>τ,φ(τt)π1(τ)π1(τt), t<τ,Ih(a,t)={BI(ta)π2(a), t>a,ψ(at)π2(a)π2(at), t<a. (24)

    Theorem 3.1. Assume

    R01.

    Then the disease-free equilibrium E0 is globally asymptotically stable.

    Proof. We will use a Lyapunov function. We adopt the Volterra-type function used in [7,10,13]. Define

    f(x)=x1lnx.

    We note that f(x)0 for all x>0. f(x) achieves its global minimum at one, with f(1)=0. Let

    q(a)=aβv(s)esa(μh+αh(σ)+rh(σ))dσds,p(τ)=βhΛhΛvμhμv(μv+αv)q(0)τm(s)esτ(μh+m(σ))dσds. (25)

    We notice that

    p(0)=R0.

    Differentiating (25) first, we obtain

    q(a)=βv(a)+(μh+αh(a)+rh(a))q(a),p(τ)=βhΛhΛvμhμv(μv+αv)q(0)m(τ)+(μh+m(τ))p(τ). (26)

    According to (26), we have Λ=βhΛhΛvμhμv(μv+αv). So we define the following Lyapunov function:

    U1(t)=U11(t)+U12(t)+U13(t)+U14(t)+U15(t), (27)

    where

    U11(t)=Λf(SvSv0),U12(t)=ΛSv0Iv(t),U13(t)=Sh0f(ShSh0),U14(t)=0p(τ)Eh(τ,t)dτ,U15(t)=Λ0q(a)Ih(a,t)da..

    Because of the complexity of the expressions, we take the derivative of each component of the Lyapunov function separately

    U11(t)=ΛSv0(1Sv0Sv)(ΛvSv0βv(a)Ih(a,t)daμvSv)=ΛSv0(1Sv0Sv)(μvSv0μvSvSv0βv(a)Ih(a,t)da)=Λμv(SvSv0)2SvSv0ΛSv0Sv0βv(a)Ih(a,t)da+Λ0βv(a)Ih(a,t)da. (28)
    U12(t)=ΛSv0[Sv0βv(a)Ih(a,t)da(μv+αv)Iv]=ΛSv0Sv0βv(a)Ih(a,t)daβhSh0Iv. (29)

    Noting that Eh(0,t)=βhShIv, we have

    U13(t)=(1Sh0Sh)(ΛhβhShIvμhSh)=(1Sh0Sh)(μhSh0μhShβhShIv)=μh(ShSh0)2ShEh(0,t)+βhSh0Iv. (30)
    U14(t)=0p(τ)Eh(τ,t)tdτ=0p(τ)[Eh(τ,t)τ+(μh+m(τ))Eh(τ,t)]dτ=[0p(τ)dEh(τ,t)+0(μh+m(τ))p(τ)Eh(τ,t)dτ]=[p(τ)Eh(τ,t)|00Eh(τ,t)dp(τ)+0(μh+m(τ))p(τ)Eh(τ,t)dτ]=p(0)Eh(0,t)Λq(0)0m(τ)Eh(τ,t)dτ=R0Eh(0,t)Λq(0)Ih(0,t). (31)

    Similarly to (31), we obtain

    U15(t)=Λ0q(a)[Ih(a,t)a+(μh+αh(a)+rh(a))Ih(a,t)]da=Λq(0)Ih(0,t)Λ0βv(a)Ih(a,t)da. (32)

    Now differentiating (27) we have

    U1(t)=Λμv(SvSv0)2SvSv0ΛSv0Sv0βv(a)Ih(a,t)da+Λ0βv(a)Ih(a,t)da+ΛSv0Sv0βv(a)Ih(a,t)daβhSh0Ivμh(ShSh0)2ShEh(0,t)+βhSh0Iv+R0Eh(0,t)Λq(0)Ih(0,t)+Λq(0)Ih(0,t)Λ0βv(a)Ih(a,t)da. (33)

    Canceling all terms that cancel, we simplify the above expression:

    U1(t)=Λμv(SvSv0)2SvSv0μh(ShSh0)2Sh+(R01)Eh(0,t). (34)

    The last inequality follows from the fact that R01. Notice that U1 equals zero implies that Sv=Sv0, Sh=Sh0, Eh(0,t)=0. We define a set

    Θ1={(Sv,Iv,Sh,Eh,Ih)Ω|U1(t)=0}.

    LaSalle's Invariance Principle [9] implies that the bounded solutions of (2) converge to the largest compact invariant set of Θ1. We will show that this largest compact invariant set is the singleton given by the disease-free equilibrium. First, we notice that equality in (34) occurs if and only if Sv=Sv0, Sh=Sh0, Eh(0,t)=0. Thus, from the solution for the equation along the characteristic line (24), we have that Eh(τ,t)=Eh(0,tτ)π1(τ)=0 for all t>τ. Hence, limtEh(τ,t)=0 for t>τ. Noting that

    Ih(0,t)=0m(τ)E(τ,t)dτ.

    So we have limtIh(0,t)=0. Thus, we have

    limtIh(a,t)=0, t>a.

    Therefore, we conclude that the disease-free equilibrium is globally stable. This completes the proof.

    Our next step is to show the global asymptotic stability of the epidemic equilibrium in system (2)


    4. The uniform strong persistence of the vector-borne disease

    In the previous section, we saw that if the reproduction number is less or equal to one, The vector-borne disease dies out. In this section, we assume that for R0>1, we will show that the vector-borne disease persists.

    From Proposition 2 we know that under the specified conditions the equilibrium E1 is locally asymptotically stable. It remains to be established that E1 is globally stable. We expect to show this result using a Lyapunov function, similar to the one used in [7,10,13]. With f(x)=x1lnx, we define the following Lyapunov function

    U2(t)=U21(t)+U22(t)+U23(t)+U24(t)+U25(t)+U26(t)+U27(t)+U28(t), (35)

    where

    {U21(t)=1q(0)0m(τ)π1(τ)dτf(SvSv),U22(t)=1Svq(0)0m(τ)π1(τ)dτIvf(IvIv),U23(t)=Shf(ShSh),U24(t)=1R00p(τ)Eh(τ)f(Eh(τ,t)Eh(τ))dτ,U25(t)=1q(0)0m(τ)π1(τ)dτ0q(a)Ih(a)f(Ih(a,t)Ih(a))da,U26(t)=tShSh(s)Eh(0,s)ds,U27(t)=tSh(s)Sh(Eh(0))2Eh(0,s)ds,U28(t)=2Eh(0)t. (36)

    One difficulty with the Lyapunov function U2 above is that the component U21 is not defined if Sv=0, the component U22 is not defined if  Iv=0, the component U23,U26 is not defined if Sh=0, the component U24 is not defined if Eh(τ,t)=0, and the component U25 is not defined if Ih(a,t)=0. To show that the Lyapunov function above is valid, we need to show that the vector borne disease persists both in the hosts and in the vectors. For this to be the case, we need to guarantee that the initial conditions we start from are non-trivial, that is, the initial conditions are such that they lead to new infections of individuals and vectors either initially or at some future point. Mathematically speaking this means that the support of the initial density of latent individuals φ(τ) intersects the support of m(\tau) or the support of the initial density of infectious individuals \psi(a) intersects the support of \beta_v(a) either initially or at some future point. Thus, we define the following set

    \displaystyle \hat{\Omega }_1=\bigg\{\varphi\in L^1_+(0, \infty)\bigg|\exists s\geq0:\ \int^\infty_0 m(\tau+s)\varphi(\tau)d\tau>0\bigg\},
    \displaystyle \hat{\Omega }_2=\bigg\{\psi\in L^1_+(0, \infty)\bigg|\exists s\geq0:\ \int^\infty_0 \beta_v(a+s)\psi(a)da>0\bigg\}.

    Define

    \displaystyle \Omega_0=\mathbb{R}_+\times \mathbb{R}_+\times \mathbb{R}_+\times \hat{\Omega }_1\times \hat{\Omega }_2.

    Finally, define X_0=\Omega\cap\Omega_0. We notice that X_0 is forward invariant. It is not hard to see that \Omega is a forward invariant. To see that \hat{\Omega }_2 is forward invariant, let us assume that the inequality holds for the initial condition. The inequality says that the condition is such that if the support of \beta_v(a) is transferred s units to the right, it will intersect the support of the initial condition. But if that happens for the initial time, it will happen for any other time since the support of the initial condition only moves to the right. Similarly, \hat{\Omega }_1 is also forward invariant.

    We want to formulate the persistence result for the vector-borne disease which on one side will justify the use of the Lyapunov functional U_2(t), and on the other, will show that when \mathcal{R}_0>1 the disease persists in the form of the endemic equilibrium. Consequently, we identify conditions which lead to the prevalence in individuals and vectors being bounded away from zero. There are many different types of persistence [21]. We identify here the two that we will be working with.

    Definition 4.1. We call the vector-borne disease uniformly weakly persistent if there exists some \gamma>0 independent of the initial conditions such that

    \displaystyle \limsup\limits_{t\rightarrow\infty}\int^\infty_0E_h(\tau, t)d\tau > \gamma \; \text{whenever}\; \int^\infty_0 \varphi(\tau)d\tau > 0,
    \displaystyle \limsup\limits_{t\rightarrow\infty}\int^\infty_0I_h(a,t)da > \gamma \; \text{whenever}\; \int^\infty_0\psi(a)da> 0,

    and

    \displaystyle \limsup\limits_{t\rightarrow\infty}I_v(t) > \gamma \; \text{whenever}\; I_{v_0}> 0.

    for all solutions of model (2).

    One of the important implications of uniform weak persistence of the disease is that the disease-free equilibrium is unstable.

    Definition 4.2. We call the vector borne diease uniformly strongly persistent if there exists some \gamma>0 independent of the initial conditions such that

    \displaystyle \liminf\limits_{t\rightarrow\infty}\int^\infty_0E_h(\tau, t)d\tau > \gamma \; \text{whenever}\; \int^\infty_0\varphi(\tau)d\tau> 0,
    \displaystyle \liminf\limits_{t\rightarrow\infty}\int^\infty_0I_h(a,t)da > \gamma \; \text{whenever}\; \int^\infty_0\psi(a)da > 0,

    and

    \displaystyle \liminf\limits_{t\rightarrow\infty}I_v(t) > \gamma \; \text{whenever} \; I_{v_0} > 0.

    for all solutions of model (2).

    It is evident from the definitions that, if the disease is uniformly strongly persistent, it is also uniformly weakly persistent. To show uniform strong persistence for the vector-borne disease, we need to show two components.

    1. We have to show that the vector-borne disease is uniformly weakly persistent.

    2. We need to show that the solution semiflow of system (2.2) has a global compact attractor \mathfrak{T}.

    First, we show uniform weak persistence of the vector-borne disease. The following proposition states that result.

    Proposition 3. Assume \mathcal{R}_0 >1. Then, for all initial conditions that belong to X_0, the vector-borne disease is uniformly weakly persistent, that is, there exists \gamma>0 such that

    \displaystyle \limsup\limits_t\beta_hI_v(t)\geq\gamma, \limsup\limits_t\int^\infty_0 m(\tau)E_h(\tau, t)d\tau\geq\gamma,
    \limsup\limits_t\int^\infty_0 \beta_v(a)I_h(a, t)da\geq\gamma.

    Proof. We argue by contradiction. Assume that the vector-borne disease dies out. In particular, assume that for every \varepsilon>0 and an initial condition in X_0 we have

    \displaystyle \limsup\limits_t\beta_hI_v(t)<\varepsilon, \limsup\limits_t\int^\infty_0 m(\tau)E_h(\tau, t)d\tau<\varepsilon, \limsup\limits_t\int^\infty_0 \beta_v(a)I_h(a, t)da<\varepsilon.

    Hence, there exist T>0 such that for all t>T, we have

    \displaystyle \beta_hI_v(t)<\varepsilon, \int^\infty_0 m(\tau)E_h(\tau, t)d\tau<\varepsilon, \int^\infty_0 \beta_v(a)I_h(a, t)da<\varepsilon.

    By shifting the dynamical system we may assume that the above inequality holds for all t\geq0. From the first equation in (2), and taking into account the above inequality, we have

    \displaystyle S_v'(t) \geq\Lambda_v-\varepsilon S_v-\mu_v S_v, S_h'(t) \geq\Lambda_h-\varepsilon S_h-\mu_h S_h.

    Therefore,

    \begin{array}{l} \mathop {\lim \sup }\limits_t {\beta _h}{I_v}(t) < \varepsilon,\quad \mathop {\lim \sup }\limits_t \int\limits_0^\infty m (\tau ){E_h}(\tau,t)d\tau < \varepsilon,\\ \quad \begin{array}{*{20}{l}} {\mathop {\quad \quad \quad \lim {{\sup }_t}}\limits_{} \int\limits_0^\infty {{\beta _v}} (a){I_h}(a,t)da < \varepsilon .} \end{array} \end{array}

    Recall that we are using the following notation B_E(t)=E_h(0, t),\ B_I(t)=I_h(0, t) . Using the inequality above we obtain

    \left\{\begin{array}{ll} \displaystyle B_E(t)=E_h(0, t)=\beta_hS_hI_v\geq\beta_h\frac{\Lambda_h}{\varepsilon+\mu_h}I_v,\\ \displaystyle \frac{dI_v(t)}{dt}\geq \frac{\Lambda_v}{\varepsilon+\mu_v}\int^\infty_0 \beta_v(a)I_h(a, t)da-(\mu_v+\alpha_v)I_v. \end{array}\right. (37)

    Now, we apply expression (24) to obtain the following system of inequalities in B_E(t),\ B_I(t) and I_v(t):

    \left\{ \begin{array}{ll} \displaystyle B_E(t)\geq\beta_h\frac{\Lambda_h}{\varepsilon+\mu_h}I_v,\\ \displaystyle B_I(t)=\int^\infty_0 m(\tau)E_h(\tau, t)d\tau\geq\int^t_0 m(\tau)B_E( t-\tau)\pi_1(\tau)d\tau,\\ \displaystyle \frac{dI_v(t)}{dt}\geq \frac{\Lambda_v}{\varepsilon+\mu_v}\int^t_0 \beta_v(a)B_I( t-a)\pi_2(a)da-(\mu_v+\alpha_v)I_v. \end{array}\right. (38)

    We will take the Laplace transform of both sides of inequalities (38). Since all functions above are bounded, their Laplace transform exists for \lambda>0. We denote by \hat{B}_E(\lambda) the Laplace transform of B_E(t), by \hat{B}_I(\lambda) the Laplace transform of B_I(t), and by \hat{I}_v(\lambda) the Laplace transform of I_v(t). Furthermore,

    \displaystyle \hat{K}_1(\lambda)=\int^\infty_0m(\tau)\pi_1(\tau)e^{-\lambda\tau}d\tau, \displaystyle \hat{K}_2(\lambda)=\int^\infty_0\beta_v(a)\pi_2(a)e^{-\lambda a}da. (39)

    Taking the Laplace transform of inequalities (38) and using the convolution property of the Laplace transform, we obtain the following system of inequalities for \hat{B}_E(\lambda),\ \hat{B}_I(\lambda) and \hat{I}_v(\lambda).

    \left\{ \begin{array}{ll} \displaystyle\hat{ B}_E(\lambda)\geq\beta_h\frac{\Lambda_h}{\varepsilon+\mu_h}\hat{I}_v(\lambda),\\ \displaystyle \hat{B}_I(\lambda)\geq\hat{K}_1(\lambda)\hat{ B}_E(\lambda),\\ \displaystyle \lambda \hat{I}_v(\lambda)-I_v(0)\geq \frac{\Lambda_v}{\varepsilon+\mu_v}\hat{K}_2(\lambda)\hat{B}_I(\lambda)-(\mu_v+\alpha_v)\hat{I}_v(\lambda). \end{array}\right. (40)

    Eliminating \hat{B}_I(\lambda) and \hat{I}_v(\lambda) from the system above, we obtain

    \begin{array}{ll} \displaystyle \hat{B}_E(\lambda)\geq\frac{\beta_h\Lambda_v\Lambda_h\hat{K}_1(\lambda)\hat{K}_2(\lambda)} {(\varepsilon+\mu_v)(\varepsilon+\mu_h)(\lambda+\mu_v+\alpha_v)}\hat{B}_E(\lambda)\displaystyle+\frac{\beta_h\Lambda_h} {(\varepsilon+\mu_h)(\lambda+\mu_v+\alpha_v)}I_v(0). \end{array}

    This last inequality should hold for the given \varepsilon\approx0 and for any \lambda>0. But this is impossible since for \varepsilon\approx0 and \lambda\approx0, the coefficient in front \hat{B}_E(\lambda) on the right hand side is approximately \mathcal{R}_0>1, that is,

    \displaystyle \frac{\beta_h\Lambda_v\Lambda_h\hat{K}_1(\lambda)\hat{K}_2(\lambda)} {(\varepsilon+\mu_v)(\varepsilon+\mu_h)(\lambda+\mu_v+\alpha_v)}\approx\mathcal{R}_0>1.

    In addition, there is another positive term on the right side of this equality. This is a contradiction with our assumption that

    \begin{align*} \displaystyle \limsup\limits_t\beta_hI_v(t)<\varepsilon, \limsup\limits_t\int^\infty_0 m(\tau)E_h(\tau, t)d\tau<\varepsilon, \end{align*}
    \begin{array}{ll} \limsup_t\int^\infty_0 \beta_v(a)I_h(a, t)da<\varepsilon. \end{array} (41)

    Therefore, there exists at least one limit supremum which is bounded below by \gamma for any initial condition in X_0 and some \gamma>0.

    Note that

    \begin{array}{ll} \displaystyle E_h(0, t)&\displaystyle= S_h \beta_hI_v(t)\leq\frac{\Lambda_h}{\mu_h} \beta_hI_v(t)\\[2ex] \displaystyle I_h(0, t)&\displaystyle=\int^\infty_0m(\tau)E_h(\tau, t)d\tau\\[2ex] &\displaystyle= \int^t_0m(\tau)E_h(0, t-\tau)\pi_1(\tau)d\tau+\int^\infty_t m(\tau)\varphi(\tau-t)\frac{\pi_1(\tau)}{\pi_1(\tau-t)}d\tau\\[2ex] \displaystyle \frac{dI_v(t)}{dt}&\displaystyle=\int^\infty_0\beta_v(a)I_h(a, t)da-(\mu_v+\alpha_v)I_v(t)\\[2ex] &\displaystyle= \int^t_0\beta_v(a)I_h(0, t-a)\pi_2(a)da+\int^\infty_t \beta_v(a)\psi(a-t)\frac{\pi_2(a)}{\pi_2(a-t)}da\\[2ex] &\displaystyle-(\mu_v+\alpha_v)I_v(t). \end{array} (42)

    Following (42), we get

    \begin{array}{ll} \displaystyle \limsup\limits_tE_h(0, t)&\displaystyle\leq\frac{\Lambda_h}{\mu_h} \limsup\limits_t\beta_hI_v(t)\\[2ex] \displaystyle \limsup\limits_tI_h(0, t)&\displaystyle\leq\int^\infty_0m(\tau)\pi_1(\tau)d\tau\limsup\limits_tE_h(0, t)\\[2ex] &\displaystyle \leq \bar{m}\int^\infty_0e^{-\mu_h\tau}d\tau\limsup\limits_tE_h(0, t)\displaystyle = \frac{\bar{m}}{\mu_h}\limsup\limits_tE_h(0, t)\\[2ex] \displaystyle \limsup\limits_t\frac{dI_v(t)}{dt}&\displaystyle\leq\int^\infty_0\beta_v(a)\pi_2(a)da\limsup\limits_tI_h(0, t)-(\mu_v+\alpha_v)\limsup\limits_tI_v(t)\\[2ex] &\displaystyle \leq m_0\int^\infty_0e^{-\mu_h a}da\limsup\limits_tI_h(0, t)-(\mu_v+\alpha_v)\limsup\limits_tI_v(t)\\[2ex] &\displaystyle=\frac{m_0}{\mu_h}\limsup\limits_tI_h(0, t)-(\mu_v+\alpha_v)\limsup\limits_tI_v(t), \end{array} (43)

    where \bar{m}=\sup_\tau\{m(\tau)\},~m_0=\sup_a\{\beta_v(a)\}. The last inequality means that

    \displaystyle \limsup\limits_tI_v(t)\leq \frac{m_0}{\mu_h(\mu_v+\alpha_v)}\limsup\limits_tI_h(0, t).

    Thus we obtain that if any inequality in (41) holds, all the three inequalities are less than a constant \times \varepsilon. There is another contradiction with the above result that there exists at least one limit supremum which is bounded below by \gamma. As a result, there exists \gamma>0 such that for any initial condition in X_0, we have

    \displaystyle \limsup\limits_t\beta_hI_v(t)\geq\gamma, \limsup\limits_t\int^\infty_0 m(\tau)E_h(\tau, t)d\tau\geq\gamma, \limsup\limits_t\int^\infty_0 \beta_v(a)I_h(a, t)da\geq\gamma.

    In addition, the differential equation for I_v can be rewritten in the form

    \displaystyle \frac{dI_v}{dt}\geq\frac{\Lambda_v\gamma}{\gamma+\mu_v}-(\mu_v+\alpha_v)I_v,

    which in turn, implies a lower bound for I_v. This concludes the proof.

    Our next goal is to prove that system (2) has a global compact attractor \mathfrak{T}. As a first step, we define the semiflow \Psi of the solutions of system (2)

    \begin{array}{ll} \displaystyle \Psi\bigg(t: S_{v_0}, I_{v_0}, S_{h_0}, \varphi(\cdot), \psi(\cdot)\bigg)=\bigg( S_v(t), I_v( t), S_h(t), E_h(\tau, t), I_h(a, t)\bigg). \end{array}

    Definition 4.3. The semiflow is a mapping \Psi: [0, \infty)\times X_0\rightarrow X_0. A set \mathfrak{T} in X_0 is called a global compact attractor for \Psi, if \mathfrak{T} is a maximal compact invariant set and if for all open sets \mathfrak{U} containing \mathfrak{T} and all bounded sets \mathcal{B} of X_0 there exists some T>0 such that \Psi(t, \mathcal{B})\subseteq\mathfrak{U}, for all t>T.

    The following proposition establishes the presence of a global compact attractor.

    Proposition 4. Assume \mathcal{R}_0>1. Then, there exists \mathfrak{T}, a compact subset of X_0, which is a global attractor for the solution semiflow \Psi of (2) in X_0. Moreover, \mathfrak{T} is invariant under the solution semiflow, that is

    \displaystyle \Psi(t, x^0)\subseteq\mathfrak{T} \; for every\; x^0\in\mathfrak{T},\ \forall t\geq 0.

    Proof. To establish this result, we will apply Lemma 3.1.3 and Theorem 3.4.6 in [22]. To show the assumptions of Lemma 3.1.3 and Theorem 3.4.6 in [22], we split the solution semiflow into two components. For an initial condition x^0\in X_0 we have \Psi(t, x^0)=\hat{\Psi}(t, x^0)+\tilde{\Psi}(t, x^0). The splitting is done in such a way that \hat{\Psi}(t, x^0)\rightarrow 0 as t\rightarrow\infty for every x^0\in X_0, and for a fixed t and any bounded set \mathcal{B} in X_0, the set \{\tilde{\Psi}(t, x^0):\ x^0\in\mathcal{B}\} is precompact. The two components of the semiflow are defined as follows:

    \begin{array}{ll} \displaystyle \hat{\Psi}\bigg(t:S_{v_0}, I_{v_0}, S_{h_0}, \varphi(\cdot), \psi(\cdot)\bigg)=\bigg( 0, 0, 0,\hat{ E}_h(\cdot, t), \hat{ I}_h(\cdot, t)\bigg)\\[2ex] \displaystyle \tilde{\Psi}\bigg(t: S_{v_0}, I_{v_0}, S_{h_0}, \varphi(\cdot), \psi(\cdot)\bigg)=\bigg( S_v(t), I_v( t), S_h(t), \tilde{E}_h(\cdot, t), \tilde{I}_h(\cdot, t)\bigg), \end{array} (44)

    where \displaystyle E_h(\tau, t)=\hat{E}_h(\tau, t)+\tilde{E}_h(\tau, t),\ I_h(a, t)=\hat{I}_h(a, t)+\tilde{I}_h(a, t) and \hat{E}_h(\tau, t),\hat{I}_h(a, t),

    \tilde{E}_h(\tau, t),\tilde{I}_h(a, t) are the solutions of the following equations (the remaining equations are as in system (2)

    \left\{ \begin{array}{ll} \displaystyle \frac{\partial \hat{E}_h}{\partial t}+\frac{\partial \hat{E}_h}{\partial \tau}=-(\mu_h+m(\tau))\hat{E}_h(\tau, t),\\[2ex] \displaystyle \hat{E}_h(0, t)=0,\\[2ex] \displaystyle \hat{E}_h(\tau, 0)=\varphi(\tau), \end{array}\right. (45)
    \left\{ \begin{array}{ll} \displaystyle \frac{\partial \hat{I}_h}{\partial t}+\frac{\partial \hat{I}_h}{\partial a}=-(\mu_h+\alpha_h(a)+r_h(a))\hat{I}_h(\tau, t),\\[2ex] \displaystyle \hat{I}_h(0, t)=0,\\[2ex] \displaystyle \hat{I}_h(a, 0)=\psi(a), \end{array}\right. (46)

    and

    \left\{ \begin{array}{ll} \displaystyle \frac{\partial \tilde{E}_h}{\partial t}+\frac{\partial \tilde{E}_h}{\partial \tau}=-(\mu+m(\tau))\tilde{E}_h(\tau, t),\\[2ex] \displaystyle \tilde{E}_h(0, t)=\beta_hS_hI_v,\\[2ex] \displaystyle \tilde{E}_h(\tau, 0)=0, \end{array}\right. (47)
    \left\{ \begin{array}{ll} \displaystyle \frac{\partial \tilde{I}_h}{\partial t}+\frac{\partial \tilde{I}_h}{\partial a}=-(\mu_h+\alpha_h(a)+r_h(a))\tilde{I}_h(\tau, t),\\[2ex] \displaystyle \tilde{I}_h(0, t)=\int^\infty_0m(\tau)\tilde{E}_h(\tau, t)d\tau,\\[2ex] \displaystyle \tilde{I}_h(\tau, 0)=0. \end{array}\right. (48)

    System (45) is decoupled from the remaining equations. Using the formula (24) to integrate along the characteristic lines, we obtain

    \begin{array}{ll} \displaystyle \hat{E}_h(\tau, t)=\left\{\begin{array}{ll} \displaystyle 0,\mbox{ } &t>\tau,\\ \displaystyle \varphi(\tau-t)\frac{\pi_1(\tau)}{\pi_1(\tau-t)},\mbox{ } &t<\tau, \end{array}\right.\\ \end{array} (49)
    \begin{array}{ll} \displaystyle \hat{I}_h(a, t)=\left\{\begin{array}{ll} \displaystyle 0,\mbox{ } &t>a,\\ \displaystyle \psi(a-t)\frac{\pi_2(a)}{\pi_2(a-t)},\mbox{ } &t<a. \end{array}\right.\\ \end{array} (50)

    Integrating \hat{E}_h with respect to \tau, we obtain:

    \displaystyle \int^\infty_t \varphi(\tau-t)\frac{\pi_1(\tau)}{\pi_1(\tau-t)}d\tau=\int^\infty_0 \varphi(\tau)\frac{\pi_1(t+\tau)}{\pi_1(\tau)}d\tau\leq e^{-\mu_h t}\int^\infty_0 \varphi(\tau)d\tau\rightarrow0,

    as t\rightarrow+\infty. Integrating \hat{I}_h with respect to a, we obtain:

    \displaystyle \int^\infty_t \psi(a-t)\frac{\pi_2(a)}{\pi_2(a-t)}da=\int^\infty_0 \psi(a)\frac{\pi_2(t+a)}{\pi_2(a)}da\leq e^{-\mu_h t}\int^\infty_0 \psi(a)da\rightarrow 0,

    as t\rightarrow+\infty. This shows the first claim, that is, it shows that \hat{\Psi}(t, x^0)\rightarrow0 as t\rightarrow+\infty uniformly for every x^0\in\mathcal{B}\subseteq X_0, where \mathcal{B} is a ball of a given radius.

    To show the second claim, we need to show compactness. We fix t and let x^0\in X_0. Note that X_0 is bounded. We have to show that for that fixed t the family of functions defined by

    \displaystyle \tilde{\Psi}(t, x^0)= \bigg( S_v(t), I_v( t), S_h(t), \tilde{E}_h(\tau, t), \tilde{I}_h(a, t)\bigg),

    obtained by taking different initial conditions in X_0 is a compact family of functions. The family

    \{\tilde{\Psi}(t, x^0)|x^0\in X_0, t-\text{fixed}\}\subseteq X_0,

    and, therefore, it is bounded. Thus, we have established the boundedness of the set. To show compactness we first see that the remaining conditions of the Frechet-Kolmogorov Theorem [19]. The third condition in the Frechet-Kolmogorov Theorem for compactness in L^1 is trivially satisfied since \tilde{E}_h(\tau, t)=0 for \tau>t and \tilde{I}_h(a, t)=0 for a>t. To see the second condition of that Theorem, we have to bound by two constants the L^1-norms of \partial E_h/\partial\tau and \partial I_h/\partial a. To derive that bound, first notice that

    \begin{array}{ll} \displaystyle \tilde{E}_h(\tau, t)=\left\{\begin{array}{ll} \displaystyle \tilde{B}_E(t-\tau)\pi_1(\tau),\mbox{ } &t>\tau,\\ \displaystyle 0,\mbox{ } &t<\tau,\\ \end{array}\right.\\ \displaystyle \tilde{I}_h(a, t)=\left\{\begin{array}{ll} \displaystyle \tilde{B}_I(t-a)\pi_2(a),\mbox{ } &t>a,\\ \displaystyle 0,\mbox{ } &t<a,\\ \end{array}\right. \end{array} (51)

    where

    \begin{array}{ll} \displaystyle \tilde{B}_E(t)=\displaystyle \beta_hS_h(t)I_v(t),\\ \displaystyle \tilde{B}_I(t)=\displaystyle \int^\infty_0m(\tau)\tilde{E}_h(\tau, t)d\tau=\int^t_0m(\tau)\tilde{B}_E(t-\tau)\pi_1(\tau)d\tau.\\ \end{array} (52)

    First, we notice that for x^0\in X_0, \tilde{B}_E(t) is bounded. We can see that by recalling that S_h and I_v are bounded. Hence, the \tilde{B}_E(t) satisfies

    \begin{array}{ll} \displaystyle \tilde{B}_E(t)\leq k_1. \end{array}

    Then, we obtain

    \displaystyle \tilde{B}_I(t)=\displaystyle\int^t_0m(\tau)\tilde{B}_E(t-\tau)\pi_1(\tau)d\tau \displaystyle\leq k_2\int^t_0\tilde{B}_E(t-\tau)d\tau= k_2\int^t_0\tilde{B}_E(\tau)d\tau \displaystyle \leq k_1k_2t.

    Next, we differentiate (51) with respect to \tau and a:

    \begin{array}{ll} \displaystyle \bigg|\frac{\partial\tilde{E}_h(\tau, t)}{\partial\tau}\bigg|\leq\left\{\begin{array}{ll} |\tilde{B}'_E(t-\tau)|\pi_1(\tau)+\tilde{B}_E(t-\tau)|\pi'_1(\tau)|,\mbox{ } &t>\tau,\\[2ex] \displaystyle 0,\mbox{ } &t<\tau, \end{array}\right.\\[2ex] \displaystyle \bigg|\frac{\partial\tilde{I}_h(a, t)}{\partial a}\bigg|\leq\left\{\begin{array}{ll} |\tilde{B}'_I(t-a)|\pi_2(a)+\tilde{B}_I(t-a)|\pi'_2(a)|,\mbox{ } &t>a,\\[2ex] \displaystyle 0,\mbox{ } &t<a. \end{array}\right. \end{array}

    We have to see that |\tilde{B}'_E(t-\tau)|,\ |\tilde{B}'_I(t-a)| are bounded. Differentiating (52), we obtain

    \begin{array}{ll} \displaystyle \tilde{B}'_E(t)=\displaystyle \beta_h\bigg(S'_h(t)I_v(t)+S_h(t)I'_v(t)\bigg),\\[2ex] \tilde{B}'_I(t)=\displaystyle m(t)\tilde{B}_E(0)\pi_1(t)+\int^t_0m(\tau)\tilde{B}'_E(t-\tau)\pi_1(\tau)d\tau. \end{array} (53)

    Taking an absolute value and bounding all terms, we can rewrite the above equality as the following inequality:

    \begin{array}{ll} \displaystyle |\tilde{B}'_E(t)|\leq k_3, |\tilde{B}'_I(t)|\leq k_4. \end{array}

    Putting all these bounds together, we have

    \begin{array}{ll} \displaystyle \parallel\partial_\tau\tilde{E}_h\parallel &\displaystyle\leq k_3\int^\infty_0\pi_1(\tau)d\tau+k_1(\mu_h+\bar{m})\int^\infty_0\pi_1(\tau)d\tau<\displaystyle \mathfrak{b}_1,\\[2ex] \displaystyle \parallel\partial_a\tilde{I}_h\parallel &\displaystyle\leq k_4\int^\infty_0\pi_2(a)da+k_1k_2(\mu_h+\bar{\alpha}_h+\bar{r}_h)t\int^\infty_0\pi_2(a)da<\displaystyle \mathfrak{b}_2, \end{array}

    where \bar{m}=\sup_\tau\{m(\tau)\},\ \bar{\alpha}_h=\sup_a\{\alpha_h(a)\},\ \bar{r}_h=\sup_a\{r_h(a)\}. To complete the proof, we notice that

    \begin{array}{ll} \displaystyle \int^\infty_0|\tilde{E}_h(\tau+h, t)-\tilde{E}_h(\tau, t)|d\tau\leq\parallel\partial_\tau\tilde{E}_h\parallel |h|\leq \mathfrak{b}_1|h|,\\[2ex] \displaystyle \int^\infty_0|\tilde{I}_h(a+h, t)-\tilde{I}_h(a, t)|d\tau\leq\parallel\partial_a\tilde{I}_h\parallel |h|\leq \mathfrak{b}_2|h|.\\ \end{array}

    Thus, the integral can be made arbitrary small uniformly in the family of functions. That establishes the second requirement of the Frechét-Kolmogorov Theorem. We conclude that the family is compact.

    Now we have all components to establish the uniform strong persistence. The next proposition states the uniform strong persistence of I_v,~E_h and I_h.

    Proposition 5. Assume \mathcal{R}_0>1. Then, for all initial conditions that belong to X_0, The vector-borne disease persists, that is, there exists \gamma>0 such that

    \displaystyle \liminf\limits_t\beta_hI_v(t)\geq\gamma, \liminf\limits_t\int^\infty_0 m(\tau)E_h(\tau, t)d\tau\geq\gamma, \liminf\limits_t\int^\infty_0 \beta_v(a)I_h(a, t)da\geq\gamma.

    Proof. We apply Theorem 2.6 in [29]. We consider the solution semiflow \Psi on X_0. We define three functionals \rho_j: X_0\rightarrowR_+,\ j=1,2,3 as follows:

    \left\{ \begin{array}{ll} \displaystyle \rho_1(\Psi(t, x^0))=\beta_hI_v(t),\\[2ex] \displaystyle \rho_2(\Psi(t, x^0))=\int^\infty_0 m(\tau)\tilde{E}_h(\tau, t)d\tau,\\[2ex] \displaystyle \rho_3(\Psi(t, x^0))=\int^\infty_0 \beta_v(a)\tilde{I}_h(a, t)da. \end{array}\right.

    Proposition 3 implies that the semiflow is uniformly weakly \rho-persistent. Proposition 4 shows that the solution semiflow has a global compact attractor \mathfrak{T}. Total orbits are solutions to the system (2) defined for all times t\in\mathbb{R}. Since the solution semiflow is nonnegative, we have that for any s and any t>s

    \begin{array}{rl} \displaystyle \beta_hI_v(t)&\displaystyle\geq \beta_hI_v(s)e^{-(\mu_v+\alpha_v)(t-s)},\\[2ex] \int^\infty_0 m(\tau)\tilde{E}_h(\tau, t)d\tau &\displaystyle =\tilde{B}_I(t)=\int^t_0m(\tau)\tilde{B}_E(t-\tau)\pi_1(\tau)d\tau \geq k^1\int^t_0\tilde{B}_E(t-\tau)d\tau\\[2ex] &\displaystyle= k^1\int^t_0\tilde{B}_E(\tau)d\tau=k^1\int^t_0\beta_hS_h(\tau)I_v(\tau)d\tau\\[2ex] &\displaystyle\geq k^2\int^t_0I_v(\tau)d\tau=k^2\int^t_0I_v(s)e^{-(\mu_v+\alpha_v)(\tau-s)}d\tau\\[2ex] &\displaystyle=\frac{k^2I_v(s)}{\mu_v+\alpha_v}e^{(\mu_v+\alpha_v)s}(1-e^{-(\mu_v+\alpha_v)t}),\\[2ex] \displaystyle \int^\infty_0 \beta_v(a)\tilde{I}_h(a, t)da&\displaystyle=\int^t_0 \beta_v(a)\tilde{B}_I(t-a)\pi_2(a) da\geq k^3\int^t_0\tilde{B}_I(t-a)da\\[2ex] &\displaystyle= k^3\int^t_0\tilde{B}_I(a)da\\[2ex] &\displaystyle\geq\frac{k^2k^3I_v(s)}{\mu_v+\alpha_v}e^{(\mu_v+\alpha_v)s}\int^t_0(1-e^{-(\mu_v+\alpha_v)a})da. \end{array}

    Therefore, \beta_hI_v(t)>0,\ \int^\infty_0 m(\tau)\tilde{E}_h(\tau, t)d\tau>0,\ \int^\infty_0\beta_v(a)\tilde{I}_h(a, t)da>0 for all t>s, provided I_v(s)>0. Theorem 2.6 in [29] now implies that the semiflow is uniformly strongly \rho-persistent. Hence, there exists \gamma such that

    \displaystyle \liminf\limits_t\beta_hI_v(t)\geq\gamma, \liminf\limits_t\int^\infty_0 m(\tau)E_h(\tau, t)d\tau\geq\gamma, \liminf\limits_t\int^\infty_0 \beta_v(a)I_h(a, t)da\geq\gamma.

    Corollary 1. Assume \mathcal{R}_0>1. There exists constants \vartheta>0 and M>0 such that for each orbit ( S_v(t), I_v, S_h(t), E_h(\tau, t), I_h(a,t)) of \Psi in \mathfrak{T}, we have

    \vartheta\leq S_v(t)\leq M, \vartheta\leq S_h(t)\leq M,\ \ \forall t\in\mathrm{R},

    and

    \displaystyle\vartheta\leq \beta_hI_v(t)\leq M, \displaystyle\vartheta\leq \int^\infty_0 m(\tau)E_h(\tau, t)d\tau\leq M,\displaystyle\vartheta\leq \int^\infty_0 \beta_v(a)I_h(a, t)da\leq M,\forall t\in\mathrm{R}.

    In the next section we show that the endemic equilibrium \mathcal{E}_1 is globally stable.


    5. Global stability of the endemic equilibrium

    Now we are ready to establish the global stability of the equilibrium \mathcal{E}_1. To demonstrate that with the Lyapunov function defined in (35) we have to establish that U'_2(t)\leq0 along the solution curves of system (2). The following Theorem summarizes the result.

    Theorem 5.1. Assume \mathcal{R}_0>1. Then, equilibrium \mathcal{E}_1 is globally asymptotically stable, that is, for any initial condition x^0\in X_0 the solution semiflow converges to \mathcal{E}_1.

    Proof. Since \mathcal{R}_0>1, for any initial condition x^0\in X_0 we can find a complete orbit (S_v(t), I_v( t), S_h(t), E_h(\tau, t), I_h(a, t)) of \Psi in \mathfrak{T} (similarly to the proof of Proposition 4) for which the inequalities in Corollary 1 hold and, consequently, there exist \varepsilon_1>0 and M_1>0 such that

    \displaystyle \varepsilon_1\leq\frac{I_v}{I^*_v}\leq M_1,\varepsilon_1\leq \frac{E_h(\tau, t)}{E_h^*(\tau)}\leq M_1, \varepsilon_1\leq \frac{I_h(a, t)}{I_h^*(a)}\leq M_1.

    This makes the Lyapunov function defined in (35) well defined.

    Because of the complexity of the expressions, we make the derivative of each component of the Lyapunov function separately (see (35)).

    \begin{split} &\displaystyle U'_{2_1}(t)\\[2ex] =&\displaystyle \frac{(1-\frac{S^*_v}{S_v})(\Lambda_v-S_v\int^\infty_0\beta_v(a)I_h(a, t)da-\mu_vS_v)}{S^*_vq(0)\int^\infty_0m(\tau)\pi_1(\tau)d\tau} \\[2ex] =&\displaystyle \frac{(1-\frac{S^*_v}{S_v})[S^*_v\int^\infty_0\beta_v(a)I^*_h(a)da+\mu_vS^*_v-S_v\int^\infty_0\beta_v(a)I_h(a, t)da-\mu_vS_v]}{S^*_vq(0)\int^\infty_0m(\tau)\pi_1(\tau)d\tau} \\[2ex] =&\displaystyle -\frac{\mu_v(S_v-S_v^*)^2}{S^*_vS_vq(0)\displaystyle\int^\infty_0m(\tau)\pi_1(\tau)d\tau}+\frac{\int^\infty_0 \beta_v(a)I^*_h(a)(1-\frac{S^*_v}{S_v}-\frac{S_vI_h(a, t)} {S^*_vI^*_h(a)}+\frac{I_h(a, t)}{I^*_h(a)})da}{q(0)\int^\infty_0m(\tau)\pi_1(\tau)d\tau}.\\ \end{split} (54)

    Next, we need to take the time derivative of U_{2_2}.

    \begin{array}{ll} \displaystyle U'_{2_2}(t) &=\displaystyle \frac{(1-\frac{I^*_v}{I_v})[S_v\int^\infty_0\beta_v(a)I_h(a, t)da-(\mu_v+\alpha_v)I_v]}{S^*_vq(0)\int^\infty_0m(\tau)\pi_1(\tau)d\tau} \\[4ex] &=\displaystyle \frac{(1-\frac{I^*_v}{I_v})(S_v\int^\infty_0\beta_v(a)I_h(a, t)da-\frac{S^*_v\int^\infty_0\beta_v(a)I^*_h(a)da}{I^*_v}I_v)}{S^*_vq(0)\int^\infty_0m(\tau)\pi_1(\tau)d\tau} \\[4ex] &=\displaystyle \frac{(1-\frac{I^*_v}{I_v})S^*_v\int^\infty_0\beta_v(a)I^*_h(a)(\frac{S_vI_h(a, t)} {S^*_vI^*_h(a)}-\frac{I_v}{I^*_v})da}{S^*_vq(0)\int^\infty_0m(\tau)\pi_1(\tau)d\tau} \\[4ex] &=\displaystyle \frac{\int^\infty_0\beta_v(a)I^*_h(a)(\frac{S_vI_h(a, t)} {S^*_vI^*_h(a)}-\frac{I_v}{I^*_v}-\frac{S_vI_h(a, t)I^*_v} {S^*_vI^*_h(a)I_v}+1)da}{q(0)\int^\infty_0m(\tau)\pi_1(\tau)d\tau},\\ \end{array} (55)

    and

    \begin{array}{ll} \displaystyle U'_{2_3}(t) &=\displaystyle \bigg(1-\frac{S^*_h}{S_h}\bigg)\bigg(\Lambda_h-\beta_hS_hI_v-\mu_hS_h\bigg)\\[3ex] &=\displaystyle \bigg(1-\frac{S^*_h}{S_h}\bigg)\bigg(E^*_h(0)+\mu_hS^*_h-E_h(0,t)-\mu_hS_h\bigg)\\[3ex] &=\displaystyle -\frac{\mu_h(S_h-S^*_h)^2}{S_h} \displaystyle+\bigg(E^*_h(0)-E_h(0,t)-\frac{S^*_h}{S_h}E^*_h(0)+\frac{S^*_h}{S_h}E_h(0,t)\bigg).\\ \end{array} (56)

    Differentiating U_{2_4}(t), we have

    \begin{array}{ll} &\displaystyle U'_{2_4}(t)\\[3ex] =&\displaystyle \frac{1}{\mathcal R_0}\int^\infty_0p(\tau)E^*_h(\tau)f'\bigg(\frac{E_h(\tau, t)}{E^*_h(\tau)}\bigg)\frac{1}{E^*_h(\tau)}\frac{\partial E_h(\tau, t)}{\partial t}d\tau\\[3ex] =&\displaystyle -\frac{1}{\mathcal R_0}\int^\infty_0\frac{p(\tau)E^*_h(\tau)}{E^*_h(\tau)}f'\bigg(\frac{E_h(\tau, t)}{E^*_h(\tau)}\bigg)\bigg(\frac{\partial E_h(\tau, t)}{\partial \tau}+(\mu_h+m(\tau))E_h(\tau, t)\bigg)d\tau\\[3ex] =&\displaystyle -\frac{1}{\mathcal R_0}\int^\infty_0p(\tau)E^*_h(\tau)df\bigg(\frac{E_h(\tau, t)}{E^*_h(\tau)}\bigg)\\[3ex] =&\displaystyle -\frac{1}{\mathcal R_0}\bigg[p(\tau)E^*_h(\tau)f\bigg(\frac{E_h(\tau, t)}{E^*_h(\tau)}\bigg)\bigg|^{\infty}_0-\int^\infty_0f\bigg(\frac{E_h(\tau, t)}{E^*_h(\tau)}\bigg)d\bigg(p(\tau)E^*_h(\tau)\bigg)\bigg]\\[3ex] =&\displaystyle \frac{1}{\mathcal R_0}\bigg[p(0)E^*_h(0)f\bigg(\frac{E_h(0, t)}{E^*_h(0)}\bigg)-\Lambda q(0)\int^\infty_0m(\tau)E^*_h(\tau)f\bigg(\frac{E_h(\tau, t)}{E^*_h(\tau)}\bigg)d\tau\bigg]\\[3ex] =&\displaystyle E^*_h(0)f\bigg(\frac{E_h(0, t)}{E^*_h(0)}\bigg)-\frac{\int^\infty_0m(\tau)E^*_h(\tau)f(\frac{E_h(\tau, t)}{E^*_h(\tau)})d\tau}{\int^\infty_0m(\tau)\pi_1(\tau)d\tau}\\[3ex] =&\displaystyle E_h(0, t)-E^*_h(0)-E^*_h(0)\ln\frac{E_h(0, t)}{E^*_h(0)}-\frac{\int^\infty_0m(\tau)E^*_h(\tau)f(\frac{E_h(\tau, t)}{E^*_h(\tau)})d\tau}{\int^\infty_0m(\tau)\pi_1(\tau)d\tau}. \end{array} (57)

    The above equality follows from (35) and the fact

    \begin{align*} &\displaystyle p'(\tau)E^*_h(\tau)+p(\tau)E'^*_h(\tau)\\[3ex] =&\displaystyle \bigg[-\Lambda q(0)m(\tau)+(\mu_h+m(\tau))p(\tau)\bigg]E^*_h(\tau)-p(\tau)(\mu_h+m(\tau))E^*_h(\tau) \\[3ex] =&\displaystyle -\Lambda q(0)m(\tau)E^*_h(\tau). \end{align*}

    We also have

    \begin{array}{ll} &\displaystyle q'(\tau)I^*_h(a)+q(a)I'^*_h(a)\\[3ex] =&\displaystyle \bigg[-\beta_v(a)+(\mu_h+\alpha_h(a)+r_h(a))q(a)\bigg]I^*_h(a)-q(a)(\mu_h+\alpha_h(a)+r_h(a))I^*_h(a) \\[3ex] =&\displaystyle -\beta_v(a)I^*_h(a). \end{array}

    Similar to the differentiation of U_{2_4}(t), we have

    \begin{array}{ll} \displaystyle U'_{2_5}(t) &=\displaystyle\frac{1}{q(0)\int^\infty_0m(\tau)\pi_1(\tau)d\tau}\displaystyle \int^\infty_0q(a)I^*_h(a)f'\bigg(\frac{I_h(a, t)}{I^*_h(a)}\bigg)\frac{1}{I^*_h(a)}\frac{\partial I_h(a, t)}{\partial t}da\\[2ex] &=\displaystyle -\frac{1}{q(0)\int^\infty_0m(\tau)\pi_1(\tau)d\tau}\int^\infty_0q(a)I^*_h(a)df\bigg(\frac{I_h(a, t)}{I^*_h(a)}\bigg)\\[2ex] &=\displaystyle \frac{q(0)I^*_h(0)f(\frac{I_h(0, t)}{I^*_h(0)})-\int^\infty_0\beta_v(a)I^*_h(a)f(\frac{I_h(a, t)}{I^*_h(a)})da}{q(0)\int^\infty_0m(\tau)\pi_1(\tau)d\tau}\\[2ex] &=\displaystyle\frac{\int^\infty_0m(\tau)E^*_h(\tau)(\frac{I_h(0, t)}{I^*_h(0)}-1-\ln\frac{I_h(0, t)}{I^*_h(0)})d\tau}{\int^\infty_0m(\tau)\pi_1(\tau)d\tau}-\frac{\int^\infty_0\beta_v(a)I^*_h(a)f(\frac{I_h(a, t)}{I^*_h(a)})da}{q(0)\int^\infty_0m(\tau)\pi_1(\tau)d\tau}.\\ \end{array} (58)

    Finally, we differentiate U_{2_6}(t),\ U_{2_7}(t) with respect to t, and we have

    \begin{array}{ll} \displaystyle U'_{2_6}(t) =\displaystyle -\frac{S^*_h}{S_h}E_h(0, t),\\ U'_{2_7}(t) =\displaystyle -\frac{S_h}{S^*_h}\frac{(E^*_h(0))^2}{E_h(0, t)}. \end{array} (59)

    Adding all five components of the Lyapunov function, we have

    \begin{array}{ll} &\displaystyle U'_2(t)\\[2ex] =&\displaystyle -\frac{\mu_v(S_v-S_v^*)^2}{S^*_vS_vq(0)\int^\infty_0m(\tau)\pi_1(\tau)d\tau}\\[2ex] &\displaystyle+\frac{1}{q(0)\int^\infty_0m(\tau)\pi_1(\tau)d\tau}\int^\infty_0 \beta_v(a)I^*_h(a)\bigg(1-\frac{S^*_v}{S_v}-\frac{S_vI_h(a, t)} {S^*_vI^*_h(a)}+\frac{I_h(a, t)}{I^*_h(a)}\bigg)da\\[2ex] &\displaystyle+ \frac{\int^\infty_0\beta_v(a)I^*_h(a)\bigg(\frac{S_vI_h(a, t)} {S^*_vI^*_h(a)}-\frac{I_v}{I^*_v}-\frac{S_vI_h(a, t)I^*_v} {S^*_vI^*_h(a)I_v}+1\bigg)da}{q(0)\int^\infty_0m(\tau)\pi_1(\tau)d\tau}\\[2ex] &\displaystyle -\frac{\mu_h(S_h-S^*_h)^2}{S_h} \displaystyle+\bigg(E^*_h(0)-E_h(0,t)-\frac{S^*_h}{S_h}E^*_h(0)+\frac{S^*_h}{S_h}E_h(0,t)\bigg)\\[2ex] &\displaystyle+ E_h(0, t)-E^*_h(0)-E^*_h(0)\ln\frac{E_h(0, t)}{E^*_h(0)}-\frac{\int^\infty_0m(\tau)E^*_h(\tau)f(\frac{E_h(\tau, t)}{E^*_h(\tau)})d\tau}{\int^\infty_0m(\tau)\pi_1(\tau)d\tau}\\[2ex] &\displaystyle+ \frac{\int^\infty_0m(\tau)E^*_h(\tau)(\frac{I_h(0, t)}{I^*_h(0)}-1-\ln\frac{I_h(0, t)}{I^*_h(0)})d\tau}{\int^\infty_0m(\tau)\pi_1(\tau)d\tau}-\frac{ \int^\infty_0\beta_v(a)I^*_h(a)f(\frac{I_h(a, t)}{I^*_h(a)})da}{q(0)\int^\infty_0m(\tau)\pi_1(\tau)d\tau}\\[2ex] & \displaystyle-\frac{S^*_h}{S_h}E_h(0, t)-\frac{S_h}{S^*_h}\frac{(E^*_h(0))^2}{E_h(0, t)}+2E^*_h(0).\\ \end{array} (60)

    Canceling all terms that cancel, we simplify (60):

    \begin{array}{ll} \displaystyle U'_2(t) &=\displaystyle -\frac{\mu_v(S_v-S_v^*)^2}{S^*_vS_vq(0)\int^\infty_0m(\tau)\pi_1(\tau)d\tau} -\frac{\mu_h(S_h-S^*_h)^2}{S_h}\\[2ex] &\displaystyle+ \frac{\int^\infty_0\beta_v(a)I^*_h(a)(3-\frac{S^*_v}{S_v}-\frac{I_v}{I^*_v}-\frac{S_vI_h(a, t)I^*_v} {S^*_vI^*_h(a)I_v}+\ln\frac{I_h(a, t)} {I^*_h(a)})da}{q(0)\int^\infty_0m(\tau)\pi_1(\tau)d\tau}\\[2ex] &\displaystyle-\frac{S^*_h}{S_h}E^*_h(0)-\frac{S_h}{S^*_h}\frac{(E^*_h(0))^2}{E_h(0, t)}-E^*_h(0)\ln\frac{E_h(0, t)}{E^*_h(0)}+2E^*_h(0)\\[2ex] &\displaystyle+ \frac{\int^\infty_0m(\tau)E^*_h(\tau)(\frac{I_h(0,t)}{I^*_h(0)}-\frac{E_h(\tau, t)}{E^*_h(\tau)}+\ln\frac{E_h(\tau, t)}{E^*_h(\tau)}\frac{I^*_h(0)}{I_h(0, t)})d\tau}{\int^\infty_0m(\tau)\pi_1(\tau)d\tau}.\\ \end{array} (61)

    Noting that

    \begin{array}{ll} \displaystyle \int^\infty_0m(\tau)E^*_h(\tau)\bigg(\frac{I_h(0, t)}{I^*_h(0)}-\frac{E_h(\tau, t)}{E^*_h(\tau)}\bigg)d\tau=0,\\[2ex] \displaystyle \int^\infty_0m(\tau)E^*_h(\tau)\bigg(\frac{E_h(\tau, t)}{E^*_h(\tau)}\frac{I^*_h(0)}{I_h(0, t)}-1\bigg)=0. \end{array} (62)

    Indeed,

    \begin{array}{ll} &\displaystyle \int^\infty_0m(\tau)E^*_h(\tau)\bigg(\frac{I_h(0, t)}{I^*_h(0)}-\frac{E_h(\tau, t)}{E^*_h(\tau)}\bigg)d\tau\\[2ex] =&\displaystyle \frac{I_h(0, t)}{I^*_h(0)}\int^\infty_0m(\tau)E^*_h(\tau)d\tau- \int^\infty_0m(\tau)E_h(\tau, t)d\tau,\\[2ex] =&\displaystyle \frac{I_h(0, t)}{I^*_h(0)}I^*_h(0)-I_h(0, t)=0,\\[2ex] &\displaystyle \int^\infty_0m(\tau)E^*_h(\tau)\bigg(\frac{E_h(\tau, t)}{E^*_h(\tau)}\frac{I^*_h(0)}{I_h(0, t)}-1\bigg)\\[2ex] =&\displaystyle \frac{I^*_h(0)}{I_h(0, t)}\int^\infty_0m(\tau)E_h(\tau, t)d\tau-\int^\infty_0m(\tau)E^*_h(\tau)d\tau\\[2ex] =&\displaystyle \frac{I^*_h(0)}{I_h(0, t)}I_h(0, t)-I^*_h(0)=0. \end{array} (63)

    Using (62) to simplify (61) we obtain

    \begin{array}{ll} \displaystyle U'_2(t) &=\displaystyle -\frac{\mu_v(S_v-S_v^*)^2}{S^*_vS_vq(0)\int^\infty_0m(\tau)\pi_1(\tau)d\tau} -\frac{\mu_h(S_h-S^*_h)^2}{S_h}\\[2ex] &\displaystyle -\frac{\int^\infty_0\beta_v(a)I^*_h(a)[f(\frac{S^*_v}{S_v}) +f(\frac{I_v}{I^*_v})+f(\frac{S_vI_h(a, t)I^*_v} {S^*_vI^*_h(a)I_v})]da}{q(0)\int^\infty_0m(\tau)\pi_1(\tau)d\tau} \\[2ex] &\displaystyle -E^*_h(0)\bigg[f\bigg(\frac{S^*_h}{S_h}\bigg)+f\bigg(\frac{S_h}{S^*_h}\frac{E^*_h(0)}{E_h(0, t)}\bigg)\bigg]\\[2ex] &\displaystyle -\frac{1}{\int^\infty_0m(\tau)\pi_1(\tau)d\tau} \int^\infty_0m(\tau)E^*_h(\tau)f\bigg(\frac{E_h(\tau, t)I^*_h(0)}{E^*_h(\tau)I_h(0, t)}\bigg)d\tau.\\ \end{array} (64)

    Hence, U_2'(t)\leq0. Define,

    \displaystyle \Theta_2= \bigg\{(S_v, I_v, S_h, E_h, I_h)\in X_0\bigg|U_2'(t)=0\bigg\}.

    We want to show that the largest invariant set in \Theta_2 is the singleton \mathcal{E}_1. First, we notice that equality in (64) occurs if and only if S_v( t)=S_v^*,\ S_h( t)=S_h^*,\ I_v( t)=I_v^*, and

    \displaystyle \frac{I_h(a, t)} {I^*_h(a)}=1, \frac{E^*_h(0)}{E_h(0, t)}=1, \frac{E_h(\tau, t)I^*_h(0)}{E^*_h(\tau)I_h(0, t)}=1. (65)

    Thus, we obtain

    \displaystyle I_h(a, t)=I_h^*(a), E_h(0, t)=E^*_h(0).

    According to (35),

    \displaystyle E_h(\tau, t)=B_E(t-\tau)\pi_1(\tau)=E_h(0, t-\tau)\pi_1(\tau)=E^*_h(0)\pi_1(\tau)=E^*_h(\tau),\ t>\tau.

    Furthermore, we obtain E_h(\tau, t)=E^*_h(\tau). We conclude that the largest invariant set in \Theta_2 is the singleton \mathcal{E}_1. Reasoning similarly to [7] can show that the compact global attractor \mathfrak{T}=\{\mathcal{E}_1\}.


    6. Discussion

    In this paper, we formulate a partial differential equation (PDE) model describing the transmission dynamics of a vector-borne disease that incorporates both incubation age of the exposed hosts and infection age of the infectious hosts. An explicit formula for the basic reproduction number \mathcal R_0 is obtained for the infection-age structured vector-host epidemic model. We show that if \mathcal R_0 of system (2) is less or equal to one, the disease-free equilibrium is locally and globally asymptotically stable. That means the disease dies out while the endemic equilibrium is not feasible. On the other hand, we show that if \mathcal R_0 is greater than one, system (2) is permanent and the endemic equilibrium is globally asymptotically stable. Therefore the disease becomes endemic. As a result, the global stability of the equilibria of system (2) is completely determined by its basic reproductive number \mathcal R_0. Hence, to control the disease, a strategy should be devised to reduce the reproduction number to below one.

    Examining the reproduction number more closely reveals that the relative impact of the recruitment rate of susceptible vectors \Lambda_v, the transmission rate \beta_h and the specific age-since-infection transmission coefficient \beta_v(a) of the infected host individuals increases \mathcal R_0. It is easy to see that \mathcal R_0 is an decreasing function of the death rate of the vector individuals \mu_v and the recovery rate of the infected vector individuals \alpha_v. It is also evident that \mathcal R_0 decreases with the rates r_h(a) and \alpha_h(a) that give recovery and disease-induced mortality of infected hosts.

    Furthermore, to see the link between \mathcal R_0 and the removal rate of the exposed host individuals with the incubation age \tau, we first need to transform the reproduction number \mathcal R_0. We will use the representation of \mathcal R_0 given in (3) and (5).

    \begin{align*} &\displaystyle \mathcal R_0\\[2ex] =&\displaystyle \frac{\beta_h\Lambda_v\Lambda_h}{\mu_v\mu_h(\mu_v+\alpha_v)}\int^\infty_0m(\tau)\pi_1(\tau)d\tau\int^\infty_0\beta_v(a)\pi_2(a)da\\[2ex] =&\displaystyle \frac{\beta_h\Lambda_v\Lambda_h}{\mu_v\mu_h(\mu_v+\alpha_v)} \int^\infty_0m(\tau)e^{-\mu_h\tau}e^{-\int^\tau_0m(\sigma)d\sigma}d\tau\int^\infty_0\beta_v(a)\pi_2(a)da\\[2ex] =&\displaystyle \frac{\beta_h\Lambda_v\Lambda_h}{\mu_v\mu_h(\mu_v+\alpha_v)}\bigg[ -\int^\infty_0(\mu_h-\mu_h-m(\tau))e^{-\mu_h\tau}e^{-\int^\tau_0m(\sigma)d\sigma}d\tau\bigg]\int^\infty_0\beta_v(a)\pi_2(a)da\\[2ex] =&\displaystyle \frac{\beta_h\Lambda_v\Lambda_h}{\mu_v\mu_h(\mu_v+\alpha_v)}\bigg[1-\mu_h \int^\infty_0e^{-\mu_h\tau}e^{-\int^\tau_0m(\sigma)d\sigma}d\tau\bigg]\int^\infty_0\beta_v(a)\pi_2(a)da.\\ \end{align*}

    Denoting by

    \displaystyle \rho=\int^\infty_0e^{-\mu_h\tau}e^{-\int^\tau_0m(\sigma)d\sigma}d\tau.

    We obtain

    \begin{array}{ll} \displaystyle \mathcal R_0= \frac{\beta_h\Lambda_v\Lambda_h}{\mu_v\mu_h(\mu_v+\alpha_v)}\int^\infty_0\beta_v(a)\pi_2(a)da\bigg(1-\mu_h \rho\bigg).\\ \end{array}

    Taking the \rho derivatives of \mathcal R_0

    \begin{array}{ll} \displaystyle \frac{d\mathcal R_0}{d\rho}&\displaystyle= -\mu_h\frac{\beta_h\Lambda_v\Lambda_h}{\mu_v\mu_h(\mu_v+\alpha_v)}\int^\infty_0\beta_v(a)\pi_2(a)da<0. \end{array}

    We have that \rho decrease with the increase of m(\tau) and \mathcal R_0 decreases with increase of \rho. Thus we have that increasing m(\tau) increases the reproduction number \mathcal R_0.

    In conclusion, our model and its analysis suggest that a better strategy of beginning mosquito control is to remove possible breeding grounds, because the larvae and pupae cycle of the mosquito is aquatic. Mosquitoes lay eggs in stagnant water, that is to say, larvae need standing water to prosper, so we must remove items that retain standing water or construct ways to keep the water moving. Furthermore, we can look for shaded rest areas used by adult mosquitoes and eliminate them. When we are outside during the day and evening hours, we can wear long sleeves and pants to prevent the bites of mosquitoes and the transmission of disease. If the infected host individuals who are in the latent period take an active drug therapy in time, the total number of the infected hosts with the virus may become small. At last it is interesting that the disease prevalence will decrease with the increase of the disease induced death rate \alpha_h(a) at the age of infection a.


    Acknowledgments

    Y. Dang is supported by NSF of Henan Province No.142300 \break 410350, Z. Qiu's research is supported by NSFC grants No. 11671206 and No. 11271190, X. Li is supported by NSF of China grant No.11271314 and Plan For Scientific Innovation Talent of Henan Province No.144200510021, and M. Martcheva is supported partially through grant NSF DMS-1220342. We are very grateful to the anonymous referees for their careful reading, valuable comments and helpful suggestions, which help us to improve the presentation of this work significantly.




    [1] D. F. Anderson, T. G. Kurtz, Continuous time markov chain models for chemical reaction networks, in Design and analysis of biomolecular circuits, Springer, (2011), 3–42. https://doi.org/10.1007/978-1-4419-6766-4_1
    [2] D. Schnoerr, G. Sanguinetti, R. Grima, Approximation and inference methods for stochastic biochemical kineticsa tutorial review, J. Phys. A: Math. Theor., 50 (2017), 093001. https://doi.org/10.1088/1751-8121/aa54d9 doi: 10.1088/1751-8121/aa54d9
    [3] S. Aoki, G. Lillacci, A. Gupta, A. Baumschlager, D. Schweingruber, M. Khammash, A universal biomolecular integral feedback controller for robust perfect adaptation, Nature, 570 (2019), 533–537. https://doi.org/10.1038/s41586-019-1321-1 doi: 10.1038/s41586-019-1321-1
    [4] T. Plesa, G. Stan, T. Ouldridge, W. Bae, Quasi-robust control of biochemical reaction networks via stochastic morphing, J. R. Soc. Interface, 18 (2021), 20200985. https://doi.org/10.1098/rsif.2020.0985 doi: 10.1098/rsif.2020.0985
    [5] J. Kim, G. Enciso, Absolutely robust controllers for chemical reaction networks, J. R. Soc. Interface, 17 (2020), 20200031. https://doi.org/10.1098/rsif.2020.0031 doi: 10.1098/rsif.2020.0031
    [6] H. Kang, T. Kurtz, Separation of time-scales and model reduction for stochastic reaction networks, Ann. Appl. Probab., 23 (2013), 529–583. https://doi.org/10.1214/12-AAP841 doi: 10.1214/12-AAP841
    [7] B. Brook, S. Waters, Mathematical challenges in integrative physiology, J. Math. Biol., 56 (2008), 893–896.
    [8] L. Preziosi, Hybrid and multiscale modelling, J. Math. Biol., 53 (2006), 977–978. https://doi.org/10.1007/s00285-006-0042-x doi: 10.1007/s00285-006-0042-x
    [9] A. Gupta, J. Mikelson, M. Khammash, A finite state projection algorithm for the stationary solution of the chemical master equation, J. Chem. Phys., 147 (2017), 154101. https://doi.org/10.1063/1.5006484 doi: 10.1063/1.5006484
    [10] J. Kuntz, P. Thomas, G. Stan, M. Barahona, Stationary distributions of continuous-time Markov chains: a review of theory and truncation-based approximations, SIAM Rev., 63 (2021), 3–64. https://doi.org/10.1137/19M1289625 doi: 10.1137/19M1289625
    [11] D. Anderson, G. Craciun, T. Kurtz, Product-form stationary distributions for deficiency zero chemical reaction networks, Bull. Math. Biol., 72 (2010), 1947–1970. https://doi.org/10.1007/s11538-010-9517-4 doi: 10.1007/s11538-010-9517-4
    [12] F. Horn, R. Jackson, General mass action kinetics, Arch. Ration. Mech. Anal., 47 (1972), 81–116.
    [13] M. Feinberg, Complex balancing in general kinetic systems, Arch. Ration. Mech. Anal., 49 (1972), 187–194. https://doi.org/10.1007/BF00255665 doi: 10.1007/BF00255665
    [14] G. Craciun, A. Dickenstein, A. Shiu, B. Sturmfels, Toric dynamical systems, J. Symb. Comput., 44 (2009), 1551–1565. https://doi.org/10.1016/j.jsc.2008.08.006 doi: 10.1016/j.jsc.2008.08.006
    [15] M. Gopalkrishnan, E. Miller, A. Shiu, A geometric approach to the global attractor conjecture, SIAM J. Appl. Dyn., 13 (2014), 758–797. https://doi.org/10.1137/130928170 doi: 10.1137/130928170
    [16] B. Joshi, A. Shiu, Atoms of multistationarity in chemical reaction networks, J. Math. Chem., 51 (2013), 153–178. https://doi.org/10.1007/s10910-012-0072-0 doi: 10.1007/s10910-012-0072-0
    [17] B. Joshi, G. Craciun, Reaction network motifs for static and dynamic absolute concentration robustness, preprint, arXiv: 2201.08428.
    [18] M. Pérez Millán, A. Dickenstein, The structure of MESSI biological systems, SIAM J. Appl. Dyn. Syst., 17 (2018), 1650–1682. https://doi.org/10.1137/17M1113722 doi: 10.1137/17M1113722
    [19] D. Cappelletti, A. Gupta, M. Khammash, A hidden integral structure endows absolute concentration robust systems with resilience to dynamical concentration disturbances, J. R. Soc. Interface, 17 (2020), 20200437. https://doi.org/10.1098/rsif.2020.0437 doi: 10.1098/rsif.2020.0437
    [20] G. Craciun, B. Joshi, C. Pantea, I. Tan, Multistationarity in cyclic sequestration-transmutation networks, Bull. Math. Biol., 84 (2022). https://doi.org/10.1007/s11538-022-01021-7 doi: 10.1007/s11538-022-01021-7
    [21] D. Cappelletti, C. Wiuf, Product-form Poisson-like distributions and complex balanced reaction systems, SIAM J. Appl. Math., 76 (2016), 411–432. https://doi.org/10.1137/15M1029916 doi: 10.1137/15M1029916
    [22] B. Joshi, A detailed balanced reaction network is sufficient but not necessary for its Markov chain to be detailed balanced, Discrete Contin. Dyn. Syst. Ser. B, 20 (2015), 1077–1105. https://doi.org/10.3934/dcdsb.2015.20.1077 doi: 10.3934/dcdsb.2015.20.1077
    [23] D. Cappelletti, B. Joshi, Graphically balanced equilibria and stationary measures of reaction networks, SIAM J. Appl. Dyn. Syst., 17 (2018), 2146–2175. https://doi.org/10.1137/17M1153315 doi: 10.1137/17M1153315
    [24] L. Hoessly, C. Mazza, Stationary distributions and condensation in autocatalytic reaction networks, SIAM J. Appl. Math., 79 (2019), 1173–1196. https://doi.org/10.1137/18M1220340 doi: 10.1137/18M1220340
    [25] D. F. Anderson, S. L. Cotter, Product-form stationary distributions for deficiency zero networks with non-mass action kinetics, Bull. Math. Biol., 78 (2016), 2390–2407. https://doi.org/10.1007/s11538-016-0220-y doi: 10.1007/s11538-016-0220-y
    [26] D. Anderson, T. Nguyen, Results on stochastic reaction networks with non-mass action kinetics, Math. Biosci. Eng., 16 (2019), 2118–2140. https://www.aimspress.com/article/doi/10.3934/mbe.2019103
    [27] F. Kelly, Reversibility and Stochastic Networks, Wiley, Chichester, 1979.
    [28] P. Whittle, Systems in Stochastic Equilibrium, John Wiley & Sons, Inc., 1986.
    [29] T. Kurtz, Limit theorems and diffusion approximations for density dependent markov chains, in Stochastic Systems: Modeling, Identification and Optimization, I, Springer, 1976, 67–78.
    [30] T. Kurtz, Strong approximation theorems for density dependent Markov chains, Stoch. Proc. Appl., 6 (1978), 223–240. https://doi.org/10.1016/0304-4149(78)90020-0 doi: 10.1016/0304-4149(78)90020-0
    [31] S. Leite, R. Williams, A constrained Langevin approximation for chemical reaction network, Ann. Appl. Probab., 29 (2019), 1541–1608.
    [32] A. Angius, G. Balbo, M. Beccuti, E. Bibbona, A. Horvath, R. Sirovich, Approximate analysis of biological systems by hybrid switching jump diffusion, Theor. Comput. Sci., 587 (2015), 49–72. https://doi.org/10.1016/j.tcs.2015.03.015 doi: 10.1016/j.tcs.2015.03.015
    [33] J. Norris, Markov Chains, Cambridge university press, 1998.
    [34] D. Anderson, D. Cappelletti, M. Koyama, T. Kurtz, Non-explosivity of stochastically modeled reaction networks that are complex balanced, Bull. Math. Biol., 80 (2018), 2561–2579. https://doi.org/10.1007/s11538-018-0473-8 doi: 10.1007/s11538-018-0473-8
    [35] G. Craciun, Toric differential inclusions and a proof of the global attractor conjecture, preprint, arXiv: 1501.02860.
    [36] L. Paulevé, G. Craciun, H. Koeppl, Dynamical properties of discrete reaction networks, J. Math. Biol., 69 (2014), 55–72. https://doi.org/10.1007/s00285-013-0686-2 doi: 10.1007/s00285-013-0686-2
  • This article has been cited by:

    1. Xia Wang, Yuming Chen, Maia Martcheva, Libin Rong, Asymptotic analysis of a vector-borne disease model with the age of infection, 2020, 14, 1751-3758, 332, 10.1080/17513758.2020.1745912
    2. Suxia Zhang, Hongsen Dong, Xiaxia Xu, Xiaoqin Shen, Analysis of a vector-borne disease model with impulsive perturbation and reinfection, 2019, 5, 2296-9020, 359, 10.1007/s41808-019-00046-9
    3. Sheng-Fu Wang, Lin Hu, Lin-Fei Nie, Global dynamics and optimal control of an age-structure Malaria transmission model with vaccination and relapse, 2021, 150, 09600779, 111216, 10.1016/j.chaos.2021.111216
    4. Sheng-Fu Wang, Lin-Fei Nie, Global Dynamics and Optimal Control of Multi-Age Structured Vector Disease Model with Vaccination, Relapse and General Incidence, 2023, 22, 1575-5460, 10.1007/s12346-022-00724-5
    5. Ting-Ting Zheng, Linfei Nie, Zhidong Teng, Yantao Luo, Shengfu Wang, Analysis of an age-structured dengue model with multiple strains and cross immunity, 2021, 14173875, 1, 10.14232/ejqtde.2021.1.50
    6. Yuan-Shun Tan, Xiao-Xue Li, Jing Yang, R. A. Cheke, Global stability of an age-structured model of smoking and its treatment, 2023, 16, 1793-5245, 10.1142/S1793524522500632
    7. Yanlin Ding, Jianjun Jiao, Global analysis and control for a vector-borne epidemic model with multi-edge infection on complex networks, 2023, 21, 2391-5455, 10.1515/math-2022-0580
    8. SHUANGSHUANG LIANG, SHENGFU WANG, LIN HU, LIN-FEI NIE, GLOBAL DYNAMICS AND OPTIMAL CONTROL FOR A VECTOR-BORNE EPIDEMIC MODEL WITH MULTI-CLASS-AGE STRUCTURE AND HORIZONTAL TRANSMISSION, 2023, 31, 0218-3390, 375, 10.1142/S0218339023500109
  • Reader Comments
  • © 2022 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(2395) PDF downloads(75) Cited by(0)

Article outline

Figures and Tables

Figures(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog