Processing math: 50%
 

Competitive exclusion in an infection-age structured vector-host epidemic model

  • Received: 01 December 2015 Accepted: 01 December 2016 Published: 01 August 2017
  • MSC : 92D30

  • The competitive exclusion principle means that the strain with the largest reproduction number persists while eliminating all other strains with suboptimal reproduction numbers. In this paper, we extend the competitive exclusion principle to a multi-strain vector-borne epidemic model with age-since-infection. The model includes both incubation age of the exposed hosts and infection age of the infectious hosts, both of which describe the different removal rates in the latent period and the variable infectiousness in the infectious period, respectively. The formulas for the reproduction numbers Rj0 of strain j,j=1,2,···,n , are obtained from the biological meanings of the model. The strain j can not invade the system if Rj0<1 , and the disease free equilibrium is globally asymptotically stable if maxj{Rj0}<1 . If Rj00>1 , then a single-strain equilibrium Ej0 exists, and the single strain equilibrium is locally asymptotically stable when Rj00>1 and Rj00>Rj0,jj0 . Finally, by using a Lyapunov function, sufficient conditions are further established for the global asymptotical stability of the single-strain equilibrium corresponding to strain j0 , which means strain j0 eliminates all other stains as long as Rj0/Rj00<bj/bj0<1,jj0 , where bj denotes the probability of a given susceptible vector being transmitted by an infected host with strain j .

    Citation: Yanxia Dang, Zhipeng Qiu, Xuezhi Li. Competitive exclusion in an infection-age structured vector-host epidemic model[J]. Mathematical Biosciences and Engineering, 2017, 14(4): 901-931. doi: 10.3934/mbe.2017048

    Related Papers:

    [1] Yan-Xia Dang, Zhi-Peng Qiu, Xue-Zhi Li, Maia Martcheva . Global dynamics of a vector-host epidemic model with age of infection. Mathematical Biosciences and Engineering, 2017, 14(5&6): 1159-1186. doi: 10.3934/mbe.2017060
    [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] 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
    [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] Abdennasser Chekroun, Mohammed Nor Frioui, Toshikazu Kuniya, Tarik Mohammed Touaoula . Global stability of an age-structured epidemic model with general Lyapunov functional. Mathematical Biosciences and Engineering, 2019, 16(3): 1525-1553. doi: 10.3934/mbe.2019073
    [6] Xinli Hu, Yansheng Liu, Jianhong Wu . Culling structured hosts to eradicate vector-borne diseases. Mathematical Biosciences and Engineering, 2009, 6(2): 301-319. doi: 10.3934/mbe.2009.6.301
    [7] Xue-Zhi Li, Ji-Xuan Liu, Maia Martcheva . An age-structured two-strain epidemic model with super-infection. Mathematical Biosciences and Engineering, 2010, 7(1): 123-147. doi: 10.3934/mbe.2010.7.123
    [8] 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
    [9] Andrey V. Melnik, Andrei Korobeinikov . Lyapunov functions and global stability for SIR and SEIR models withage-dependent susceptibility. Mathematical Biosciences and Engineering, 2013, 10(2): 369-378. doi: 10.3934/mbe.2013.10.369
    [10] Lu Gao, Yuanshun Tan, Jin Yang, Changcheng Xiang . Dynamic analysis of an age structure model for oncolytic virus therapy. Mathematical Biosciences and Engineering, 2023, 20(2): 3301-3323. doi: 10.3934/mbe.2023155
  • The competitive exclusion principle means that the strain with the largest reproduction number persists while eliminating all other strains with suboptimal reproduction numbers. In this paper, we extend the competitive exclusion principle to a multi-strain vector-borne epidemic model with age-since-infection. The model includes both incubation age of the exposed hosts and infection age of the infectious hosts, both of which describe the different removal rates in the latent period and the variable infectiousness in the infectious period, respectively. The formulas for the reproduction numbers Rj0 of strain j,j=1,2,···,n , are obtained from the biological meanings of the model. The strain j can not invade the system if Rj0<1 , and the disease free equilibrium is globally asymptotically stable if maxj{Rj0}<1 . If Rj00>1 , then a single-strain equilibrium Ej0 exists, and the single strain equilibrium is locally asymptotically stable when Rj00>1 and Rj00>Rj0,jj0 . Finally, by using a Lyapunov function, sufficient conditions are further established for the global asymptotical stability of the single-strain equilibrium corresponding to strain j0 , which means strain j0 eliminates all other stains as long as Rj0/Rj00<bj/bj0<1,jj0 , where bj denotes the probability of a given susceptible vector being transmitted by an infected host with strain j .


    1. Introduction

    In many infectious diseases, such as HIV, schistosomiasis, tuberculosis, the infectiousness of an infected individual can be very different at various stages of infection. Hence, the age of infection may be an important factor to consider in modeling transmission dynamics of infectious diseases. In the epidemic model of Kermack and Mckendrick [9], infectivity is allowed to depend on the age of infection. Because the age-structured epidemic model is described by first order PDEs, it is more difficult to theoretically analyze the dynamical behavior of the PDE models, particularly the global stability. Several recent studies [10,11,18] have focused on age structured models, and the results show that age of infection may play an important role in the transmission dynamics of infectious diseases.

    In our pervious work [4], we formulated an infection-age structured epidemic model to describe the transmission dynamics of vector-borne diseases. The model includes both incubation age of the exposed hosts and infection age of the infectious hosts, both of which describe the different removal rates in the latent period and the variable infectiousness in the infectious period, respectively. The results in [4] show that the basic reproduction number determines transmission dynamics of vector-borne diseases: the disease-free equilibrium is globally asymptotically stable if the basic reproduction number is less than 1, and the endemic equilibrium is globally asymptotically stable if the basic reproduction number is greater than 1. However, the vector-borne epidemic model formulated in [4] only incorporates a single strain. In reality, many diseases are caused by more than one antigenically different strains of the causative agent [15]. For instance, the dengue virus has 4 different serotypes [6], and bacterial pneumonia is caused by more than ninety different serotypes of Streptoccus pneumoniae. Therefore, it is necessary to study infection-age structured epidemic models with multiple strains.

    In this paper, we will extend the model with a single strain to the model with multiple strains, and obtain the following infection-age-structured vector-borne epidemic model with multiple strains:

    {dSvdt=Λvnj=1Sv0βjv(a)Ijh(a,t)daμvSv,dIjvdt=Sv0βjv(a)Ijh(a,t)da(μv+αjv)Ijv,dRvdt=nj=1αjvIjvμvRv,dShdt=Λhnj=1βjhShIjvμhSh,Ejh(τ,t)τ+Ejh(τ,t)t=(μh+mjh(τ))Ejh(τ,t),Ejh(0,t)=βjhShIjv,Ijh(a,t)a+Ijh(a,t)t=(μh+αjh(a)+rjh(a))Ijh(a,t),Ijh(0,t)=0mjh(τ)Ejh(τ,t)dτ,dRhdt=nj=10rjh(a)Ijh(a,t)daμhRh. (1)

    In the model (1), Sh(t), Ejh(τ,t), Ijh(a,t), Rh(t) represent the number/density of the susceptible hosts, infected hosts with strain j but not infectious, infected hosts with strain j and infectious, and recovered hosts at time t, respectively. Sv(t), Ijv(t) and Rv(t) denote the number of the susceptible vectors, infected vectors with strain j and infectious, and recovered vectors at time t, respectively. Λv,Λh are the birth /recruitment rates of the vectors and hosts, respectively; μv,μh are the natural death rates of the vectors and hosts, respectively. The parameter mjh(τ) denotes the removal rate of the infected hosts with strain j of incubation age τ from the latent period; αjh(a) is the additional disease induced death rate due to the strain j at age of infection a; αjv denotes the recovery rate of the infected vectors with strain j; rjh(a) denotes the recovery rate of the infected hosts of infection age a with strain j; βjv(a) is the transmission coefficient of the infected host individuals with strain j at age of infection a, and βjh is the transmission coefficient of strain j from infected vectors to healthy host individuals.

    The dynamics of the epidemic model involving multiple strains has fascinated researchers for a long time (see [3,5,6,7,17] and the references therein), and one of the important results is the competitive exclusion principle. In epidemiology, the competitive exclusion principle states that if multiple strains circulate in the population, only the strain with the largest reproduction number persists and the strains with suboptimal reproduction numbers are eliminated [13]. Using a multiple-strain ODE model Bremermann and Thieme [2] first proved that the principle of competitive exclusion is valid under the assumption that infection with one strain precludes additional infections with other strains. In 2013, Maracheva and Li [13] extended the competitive exclusion principle to a multi-stain age-since-infection structured model of SIR/SI-type. The goal of this paper is to extend this principle to model (1).

    As we all know, the proof of competitive exclusion principle is based on the global stability of the single-strain equilibrium. The stability analysis of nonlinear dynamical systems has always been an important topic theoretically and practically since global stability is one of the most important issues related to their dynamic behaviors. Due to the lack of generically applicable tools proving the global stability is very challenging, especially for the continuous age-structured models which are described by first order PDEs. Although there are various approaches for some general nonlinear systems, the method of Lyapunov functions is the most common tool used to prove the global stability. In this paper, we will apply a class of Lyaponuv functions to study the global dynamics of system (1) and draw on the results to derive the competitive exclusion principle for infinite dimensional systems.

    This paper is organized as follows. In the next section we derive an explicit formula for the basic reproduction number Rj0 of strain j for j=1,,n, and then we will show that strain j will die out if its basic reproduction number is less than one. In section 3, we will define the disease reproduction number R0, and then prove that the disease-free equilibrium (DFE) of the system is globally asymptotically stable if R0<1. In Section 4, we will investigate the existence of single-strain equilibria and their local stabilities. In section 5 we will devote to prove the principle of competitive exclusion. Without loss of generality, we assume that strain one has the maximal reproduction number and R10>1. Under the assumption, we will show that strain one is uniformly strong persistent while the remaining strains become extinct. In Section 6, we use a class of Lyapunov functions to derive the global stability of the strain one equilibrium under the condition that Ri0/R10<bi/b1<1,i1, where bj denotes the probability of a given susceptible vector being transmitted by an infected host with strain j, which implies that complete competitive exclusion holds for the system. Finally, a brief discussion is given in Section 7.


    2. The reproduction numbers and threshold dynamics

    In this section, we mainly derive the reproduction numbers for each strain, and show that the stain will die out if its basic reproduction number is less than one.

    Since the equations for the recovered individuals and the recovered vectors are decoupled from the system, it follows that the dynamical behavior of system (1) is equivalent to the dynamical behavior of the following system:

    {dSvdt=Λvnj=1Sv0βjv(a)Ijh(a,t)daμvSv,dIjvdt=Sv0βjv(a)Ijh(a,t)da(μv+αjv)Ijv,dShdt=Λhnj=1βjhShIjvμhSh,Ejh(τ,t)τ+Ejh(τ,t)t=(μh+mjh(τ))Ejh(τ,t),Ejh(0,t)=βjhShIjv,Ijh(a,t)a+Ijh(a,t)t=(μh+αjh(a)+rjh(a))Ijh(a,t),Ijh(0,t)=0mjh(τ)Ejh(τ,t)dτ. (2)

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

    Sv(0)=Sv0,Ijv(0)=Ijv0,Sh(0)=Sh0,Ejh(τ,0)=φj(τ),Ijh(a,0)=ψj(a).

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

    Assumption 2.1.

    1. The function βjv(a) is bounded and uniformly continuous for every j. When βjv(a) is of compact support, the support has non-zero Lebesgue measure;

    2. The functions mjh(τ), αjh(a), rjh(a) belong to L(0,);

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

    Let us define

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

    It is easily verified that solutions of (2) with nonnegative initial conditions belong to the positive cone for t0. Adding the first and all equations for Ijv yields that

    ddt(Sv(t)+nj=1Ijv(t))Λvμv(Sv(t)+nj=1Ijv(t)).

    Hence,

    lim supt+(Sv(t)+nj=1Ijv(t))Λvμv.

    Similarly, adding the equation for Sh and all equations for Ejh,Ijh, we have

    ddt(Sh(t)+nj=10Ejh(τ,t)dτ+nj=10Ijh(a,t)da)
    Λhμh(Sh(t)+nj=10Ejh(τ,t)dτ+nj=10Ijh(a,t)da),

    and it then follows that

    lim supt+(Sh(t)+nj=10Ejh(τ,t)dτ+nj=10Ijh(a,t)da)Λhμh.

    Therefore, the following set is positively invariant for system (2)

    Ω={(Sv, I1v, , Inv, Sh, E1h, I1h, , Enh, Inh)X+|(Sv(t)+nj=1Ijv(t))Λvμv,(Sh(t)+nj=10Ejh(τ,t)dτ+nj=10Ijh(a,t)da)Λhμh}. (3)

    In what follows, we only consider the solutions of the system (2) with initial conditions which lie in the region Ω. As we all know, the reproduction number is one of most important concepts in epidemiological model. Next, we will express the basic reproduction numbers for each strain. To simplify expression, let us introduce two notations.

    Definition 2.1. The exit rate of exposed host individuals with strain j from the incubation compartment is given by μh+mjh(τ), the probability of still being latent after τ time units, denoted by πj1(τ), is given by

    πj1(τ)=eμhτeτ0mjh(σ))dσ. (4)

    Definition 2.2. The exit rate of infected individuals with strain j from the infective compartment is given by μh+αjh(a)+rjh(a), and it then follows that the probability of still being infectious after a time units, denoted by πj2(a), is given by

    πj2(a)=eμhaea0(αjh(σ)+rjh(σ))dσ. (5)

    Then we can give the expression for the basic reproduction number of strain j which can be expressed as

    Rj0=βjhΛvΛhμvμh(μv+αjv)0mjh(τ)πj1(τ)dτ0βjv(a)πj2(a)da. (6)

    The reproduction number of strain j gives the number of secondary infections produced in an entirely susceptible population by a typical infected individual with strain j during its entire infectious period. Rj0 gives the strength of strain j to invade into the system when rare and alone. The reproduction number of strain j consists of two terms:

    Rjh=Λvμv0βjv(a)πj2(a)da,Rjv=βjhΛhμh(μv+αjv)0mjh(τ)πj1(τ)dτ.

    The first term Rjh represents the reproduction number of human-to-vector transmission of strain j, and the second term Rjv is the reproduction number of vector-to-human transmission of strain j.

    Now we are able to state the results on threshold dynamics of strain j:

    Theorem 2.3. If Rj0<1, strain j will die out.

    Proof. Let

    BjE(t)=Ejh(0,t),BjI(t)=Ijh(0,t).

    Integrating along the characteristic lines of system (2) yields

    Ejh(τ,t)={BjE(tτ)πj1(τ), t>τ,φj(τt)πj1(τ)πj1(τt), t<τ,Ijh(a,t)={BjI(ta)πj2(a), t>a,ψj(at)πj2(a)πj2(at), t<a. (7)

    From the first and the third equations of system (2), we obtain

    lim supt+Sv(t)Λvμv,lim supt+Sh(t)Λhμh. (8)

    Thus, from system (2) and inequalities (8), we have

    {dIjv(t)dtΛvμv0βjv(a)Ijh(a,t)da(μv+αjv)Ijv,Ejh(τ,t)=Ejh(0,tτ)πj1(τ),t>τ,Ijh(a,t)=Ijh(0,ta)πj2(a),t>a. (9)

    From the first inequality of (9), we obtain that

    Ijv(t)Ijv(0)e(μv+αjv)t+Λvμvt0e(μv+αjv)(ts)0βjv(a)Ijh(a,s)dadsIjv(0)e(μv+αjv)t+Λvμvt0e(μv+αjv)(ts)(s0βjv(a)Ijh(0,sa)πj2(a)da+tsβjv(a)ψj(as)πj2(a)πj2(as)da+tβjv(a)Ijh(a,s)da)ds. (10)

    Notice that

    lim supt+t0e(μv+αjv)(ts)s0βjv(a)Ijh(0,sa)πj2(a)dads(lim supt+t0e(μv+αjv)(ts)ds)0βjv(a)πj2(a)da(lim supt+Ijh(0,t))=1μv+αjv0βjv(a)πj2(a)da(lim supt+Ijh(0,t)), (11)
    lim supt+t0e(μv+αjv)(ts)tsβjv(a)ψj(as)πj2(a)πj2(as)dadsˉβlim supt+t0e(μv+αjv)(ts)tsψj(as)eaas(μh+αjh(σ)+rjh(σ))dσdadsˉβlim supt+t0e(μv+αjv)(ts)tsψj(as)eμhsdads
    =ˉβlim supt+(e(μv+αjv)tt0e(μv+αjvμh)sts0ψj(a)dads)=ˉβ0ψj(a)dalim supt+(e(μv+αjv)te(μv+αjvμh)t1μv+αjvμh)=0, (12)

    and

    lim supt+t0e(μv+αjv)(ts)tβjv(a)Ijh(a,s)dads=0. (13)

    It then follows from (11), (12) and (13) that

    lim supt+Ijv(t)Λvμv(μv+αjv)0βjv(a)πj2(a)da(lim supt+Ijh(0,t))Λvμv(μv+αjv)0βjv(a)πj2(a)da(lim supt+0mjh(τ)Ejh(τ,t)dτ)Λvμv(μv+αjv)0βjv(a)πj2(a)dalim supt+(t0mjh(τ)Ejh(τ,t)dτ+tmjh(τ)Ejh(τ,t)dτ)=Λvμv(μv+αjv)0βjv(a)πj2(a)da(lim supt+t0mjh(τ)Ejh(0,tτ)πj1(τ)dτ)Λvμv(μv+αjv)0βjv(a)πj2(a)da0mjh(τ)πj1(τ)dτ(lim supt+Ejh(0,t))βjhΛvΛhμvμh(μv+αjv)0βjv(a)πj2(a)da0mjh(τ)πj1(τ)dτlim supt+Ijv(t)Rj0lim supt+Ijv(t). (14)

    Since Rj0<1 and Ijv(t), j=1,,n, are all bounded, the above expression implies that

    lim supt+Ijv(t)=0,j=1,,n. (15)

    Hence, we have

    lim supt+Ejh(0,t)=0,lim supt+Ejh(τ,t)=lim supt+Ejh(0,tτ)πj1(τ)=0. (16)

    By using the same argument, we have

    lim supt+Ijh(0,t)=0,lim supt+Ijh(a,t)=0. (17)

    Therefore, (Ijv(t),Ejh(τ,t),Ijh(a,t))0 as t. This means that strain j will die out. The proof of Theorem 2.3 is completed.


    3. Global stability of the disease-free equilibrium

    In this section, we mainly define the disease reproduction number and show that the disease free equilibrium is globally asymptotically stable if the disease reproduction number R0 is less than one, where

    R0=max{R10,,Rn0}.

    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,

    and 0=(0,,0) is an n-dimensional zero vector.

    Now let us establish the local stability of the disease-free equilibrium. Let

    Sv(t)=Sv0+xv(t),Ijv(t)=yjv(t),Sh(t)=Sh0+xh(t),Ejh(τ,t)=zjh(τ,t),Ijh(a,t)=yjh(a,t).

    Then the linearized system of system (2) at the disease-free equilibrium E0 can be expressed as

    {dxv(t)dt=nj=1Sv00βjv(a)yjh(a,t)daμvxv(t),dyjv(t)dt=Sv00βjv(a)yjh(a,t)da(μv+αjv)yjv(t),dxh(t)dt=nj=1βjhSh0yjv(t)μhxh(t),zjh(τ,t)τ+zjh(τ,t)t=(μh+mjh(τ))zjh(τ,t),zjh(0,t)=βjhSh0yjv(t),yjh(a,t)a+yjh(a,t)t=(μh+αjh(a)+rjh(a))yjh(a,t),yjh(0,t)=0mjh(τ)zjh(τ,t)dτ. (18)

    Let

    yjv(t)=ˉyjveλt, zjh(τ,t)=ˉzjh(τ)eλt, yjh(a,t)=ˉyjh(a)eλt, (19)

    where ˉyjv,ˉzjh(τ) and ˉyjh(a) are to be determined. Substituting (19) into (18), we obtain

    {λˉyjv=Sv00βjv(a)ˉyjh(a)da(μv+αjv)ˉyjv,dˉzjh(τ)dτ=(λ+μh+mjh(τ))ˉzjh(τ),ˉzjh(0)=βjhSh0ˉyjv,dˉyjh(a)da=(λ+μh+αjh(a)+rjh(a))ˉyjh(a),ˉyjh(0)=0mjh(τ)ˉzjh(τ)dτ. (20)

    Solving the differential equation, we obtain

    ˉzjh(τ)=ˉzjh(0) eλτπj1(τ)=βjhSh0ˉyjv eλτπj1(τ).

    Substituting the expression for ˉzjh(τ) into the equation for ˉyjh(0), expressing ˉyjh(0) in term of ˉzjh(0), and replacing ˉyjh(0) in the equation for ˉyjh(a), we obtain

    ˉyjh(a)=ˉyjh(0) eλaπj2(a)=βjhSh0ˉyjv eλaπj2(a)0mjh(τ) eλτπj1(τ)dτ.

    Substituting the above expression for ˉyjh(a) into the first equation of (20), we can obtain

    λ+μv+αjv=βjhSv0Sh00mjh(τ)eλτπj1(τ)dτ0βjv(a)eλaπj2(a)da. (21)

    Now we are able to state the following result.

    Theorem 3.1. If

    R0=max{R10,,Rn0}<1,

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

    Proof. We first prove the first result. Let us assume R0<1. For ease of notation, set

    LHSdef=λ+μv+αjv,RHSdef=G1(λ)=βjhSv0Sh00mjh(τ)eλτπj1(τ)dτ0βjv(a)eλaπj2(a)da. (22)

    We can easily verify that

    |LHS|μv+αjv,|RHS|G1(λ)G1(0)=βjhSv0Sh00mjh(τ)πj1(τ)dτ0βjv(a)πj2(a)da=βjhΛvΛhμvμh0mjh(τ)πj1(τ)dτ0βjv(a)πj2(a)da=Rj0(μv+αjv)<|LHS|,

    for any λ,λ0. There is a contradiction. The contradiction implies that the equation (21) cannot have any roots with non-negative real parts. Hence, the disease-free equilibrium is locally asymptotically stable.

    Next, let us assume max{R10,,Rn0}=Rj00>1. We rewrite the characteristic equation (21) in the form

    G2(λ)=0, (23)

    where

    G2(λ)=(λ+μv+αj0v)βj0hSv0Sh00mj0h(τ)eλτπj01(τ)dτ0βj0v(a)eλaπj02(a)da.

    It is easily verified that

    G2(0)=(μv+αj0v)βj0hSv0Sh00mj0h(τ)πj01(τ)dτ0βj0v(a)πj02(a)da=(μv+αj0v)(1Rj00)<0,

    and

    limλ+G2(λ)=+.

    Hence, the characteristic equation (23) has a real positive root. Therefore, the disease free equilibrium E0 is unstable. This concludes the proof.

    We have proved that the disease-free equilibrium is locally stable if R0<1. It also follows from Theorem 2.3 that strain j will die out if Rj0<1. Therefore we have the following result.

    Theorem 3.2. If

    R0=max{R10,,Rn0}<1,

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


    4. Existence and stability of boundary equilibria

    In this section, we mainly investigate the existence and stability of the boundary equilibria. For ease of notation, let

    Δj=βjhΛhΛvμhμv(μv+αjv),bj=0mjh(τ)πj1(τ)dτ0βjv(a)πj2(a)da,bj(λ)=0mjh(τ)eλτπj1(τ)dτ0βjv(a)eλaπj2(a)da. (24)

    From Theorem 2.3, it follows that strain j will die out if Rj0<1. Thus in later sections we always assume that Rj0>1 for all j,j=1,2,,n. If Rj0>1, straightforward computation yields that system (2) has a corresponding single-strain equilibrium Ej which is given by

    Ej=(Sjv,0,,0,Ijv,0,,0,Sjh,0,,0,Ejh(τ),Ijh(a),0,,0).

    The non-zero components Ijv,Ejh and Ijh are in positions j+1,n+2j+1 and n+2j+2, respectively, and

    Ijv=μvμh(Rj01)βjh(Λhbj+μv),Sjv=Λv(μv+αjv)Ijvμv=βjhΛv(μv+Λhbj)μvμh(μv+αjv)(Rj01)βjhμv(μv+Λhbj),
    Sjh=ΛhβjhIjv+μh=Λh(μv+Λhbj)μh(μvRj0+Λhbj),Ejh(τ)=Ejh(0)πj1(τ),Ejh(0)=βjhSjhIjv,Ijh(a)=Ijh(0)πj2(a),Ijh(0)=Ejh(0)0mjh(τ)πj1(τ)dτ. (25)

    The results on the local stability of single-strain equilibrium Ej0 are summarized below:

    Theorem 4.1. Assume Rj00>1 for a fixed j0 and

    Rj0<Rj00foralljj0.

    Then single-strain equilibrium Ej0 is locally asymptotically stable. If there exists i0 such that

    Ri00>Rj00,

    then the single-strain equilibrium Ej0 is unstable.

    Proof. Without loss of generality, we assume that R10>1 and Ri0<R10 for i=2,,n. Let

    Sv(t)=S1v+xv(t), Sh(t)=S1h+xh(t),I1v(t)=I1v+y1v(t), E1h(τ,t)=E1h(τ)+z1h(τ,t), I1h(a,t)=I1h(a)+y1h(a,t),Iiv(t)=yiv(t),Eih(τ,t)=zih(τ,t), Iih(a,t)=yih(a,t),

    where i=2,,n. Then the linearization system of system (2) at the equilibrium E1 can be expressed as

    {dxv(t)dt=S1v0β1v(a)y1h(a,t)daxv(t)0β1v(a)I1h(a)dani=2S1v0βiv(a)yih(a,t)daμvxv(t),dy1v(t)dt=S1v0β1v(a)y1h(a,t)da+xv(t)0β1v(a)I1h(a)da(μv+α1v)y1v(t),dyiv(t)dt=S1v0βiv(a)yih(a,t)da(μv+αiv)yiv(t),dxh(t)dt=β1hS1hy1v(t)β1hxh(t)I1vni=2βihS1hyiv(t)μhxh(t),zjh(τ,t)τ+zjh(τ,t)t=(μh+mjh(τ))zjh(τ,t),z1h(0,t)=β1hS1hy1v(t)+β1hxh(t)I1v,zih(0,t)=βihS1hyiv(t),yjh(a,t)a+yjh(a,t)t=(μh+αjh(a)+rjh(a))yjh(a,t),yjh(0,t)=0mjh(τ)zjh(τ,t)dτ. (26)

    An approach similar to [14] (see Appendix B in [14]) can show that the linear stability of the system is determined by the eigenvalues of the linearized system (26). In order to investigate the linear stability of the linearized system (26), we consider exponential solutions (see the case of the disease-free equilibrium) and obtain a linear eigenvalue problem. For the whole system, we only consider the equations for strains i,i=2,,n, and obtain the following eigenvalue problem:

    {dyiv(t)dt=S1v0βiv(a)yih(a,t)da(μv+αiv)yiv(t),zih(τ,t)τ+zih(τ,t)t=(μh+mih(τ))zih(τ,t),zih(0,t)=βihS1hyiv(t),yih(a,t)a+yih(a,t)t=(μh+αih(a)+rih(a))yih(a,t),yih(0,t)=0mih(τ)zih(τ,t)dτ. (27)

    For each i,i1, by using the same argument to equation (21), we obtain the following characteristic equation

    λ+μv+αiv=βihS1vS1h0mih(τ)eλτπi1(τ)dτ0βiv(a)eλaπi2(a)da. (28)

    Notice that Sjv and Sjh satisfy

    βjhSjvSjh0mjh(τ)πj1(τ)dτ0βjv(a)πj2(a)da=μv+αjv, (29)

    for j=1,,n. It then follows from (6) and (24) that we have

    S1vS1h=μv+α1vβ1hb1=ΛvΛhμvμhR10. (30)

    Substituting (30) into the equation (28), we get the following characteristic equation

    λ+μv+αiv=βihΛvΛhμvμhR10bi(λ), (31)

    where bi(λ) is defined in (24).

    First, assume that Ri00>R10 for some i0, and set

    Gi0(λ)def=(λ+μv+αi0v)βi0hΛvΛhμvμhR10bi0(λ).

    Straightforward computation yields that

    Gi0(0)=(μv+αi0v)βi0hΛvΛhμvμhR10bi0=(μv+αi0v)(1Ri00R10)<0.

    Furthermore, for λ real, Gi0(λ) is an increasing function of λ such that limGi0(λ)+ as λ+. Hence Intermediate Value Theorem implies that the equation (31) has a unique real positive solution. We conclude that in that case E1 is unstable.

    Next, assume Ri0<R10 for all i=2,,n, and set

    G3(λ)def=λ+μv+αiv,G4(λ)def=βihΛvΛhμvμhR10bi(λ). (32)

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

    |G3(λ)|μv+αiv,|G4(λ)|G4(λ)G4(0)=1R10βihΛvΛhμvμh0mih(τ)πi1(τ)dτ0βiv(a)πi2(a)da=Ri0R10(μv+αiv)<|G3(λ)|.

    This gives a contradiction. Hence, the equation (31) have no solutions with positive real part and all eigenvalues of these equations have negative real parts. Therefore, the stability of E1 depends on the eigenvalues of the following system

    {λxv=S1v0β1v(a)y1h(a)daxv0β1v(a)I1h(a)daμvxv,λy1v=S1v0β1v(a)y1h(a)da+xv0β1v(a)I1h(a)da(μv+α1v)y1v,λxh=z1h(0)μhxh,dz1h(τ)dτ=(λ+μh+m1h(τ))z1h(τ),z1h(0)=β1hS1hy1v+β1hI1vxh,dy1h(a)da=(λ+μh+α1h(a)+r1h(a))y1h(a),y1h(0)=0m1h(τ)z1h(τ)dτ. (33)

    Solving the differential equation, we have

    z1h(τ)=z1h(0) eλτπ11(τ),y1h(a)=y1h(0) eλaπ12(a)=z1h(0) eλaπ12(a)0m1h(τ) eλτπ11(τ)dτ.

    Substituting the above expression for y1h(a) into the first and the second equations of (33) yileds that

    {(λ+μv+0β1v(a)I1h(a)da)xv+S1vb1(λ)z1h(0)=0,xv0β1v(a)I1h(a)da+(λ+μv+α1v)y1vS1vb1(λ)z1h(0)=0,(λ+μh)xh+z1h(0)=0,β1hI1vxhβ1hS1hy1v+z1h(0)=0. (34)

    Direct calculation yields the following characteristic equation

    (λ+μv+0β1v(a)I1h(a)da)(λ+μv+α1v)(λ+μh+β1hI1v)=β1hS1hS1vb1(λ)(λ+μv)(λ+μh). (35)

    Dividing both sides by (λ+μv)(λ+μh) gives

    G5(λ)=G6(λ), (36)

    where

    G5(λ)=(λ+μv+0β1v(a)I1h(a)da)(λ+μv+α1v)(λ+μh+β1hI1v)(λ+μv)(λ+μh),G6(λ)=β1hS1hS1vb1(λ)=β1hS1hS1v0m1h(τ)eλτπ11(τ)dτ0β1v(a)eλaπ12(a)da. (37)

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

    |G5(λ)|>|λ+μv+α1v|μv+α1v. (38)

    From (29), we have

    |G6(λ)||G6(λ)|G6(0)=β1hS1hS1v0m1h(τ)π11(τ)dτ0β1v(a)π12(a)da=μv+α1v<|G5(λ)|. (39)

    This leads to a contradiction. The contradiction implies that (36) has no roots such that λ0. Thus, the characteristic equation for strain one has only roots with negative real parts. Thus, the single strain equilibrium E1 is locally asymptotically stable if R10>1 and Ri0<R10, i=2,,n. This concludes the proof.


    5. Preliminary results and uniform persistence

    In the previous section, we proved that if the disease reproduction number is less than one, all strains are eliminated and the disease dies out. Our next step is to show that the competitive exclusion principle holds for system (2). In the later sections, we always assume that R0>1. Without loss of generality, we assume that

    R10=max{R10,,Rn0}>1.

    In the following we will show that strain 1 persists while the other strains die out if Ri0/R10<bi/b1<1,i1, where bj denotes the probability of a given susceptible vector being transmitted by an infected host with strain j. Hence, the strain with the maximal reproduction number eliminates all the rest and the competitive exclusion principle will be established for system (2).

    Mathematically speaking, establishing the competitive exclusion principle means establishing the global stability of the single-strain equilibrium E1. From Theorem 4.1 we know that if Ri0/R10<1,i1, the equilibrium E1 is locally asymptotically stable. In the following we only need to show that E1 is a global attractor. The method used here to show this result is similar to the one used in [1,12,13,20].

    Set

    f(x)=x1lnx.

    It is easy to check that f(x)0 for all x>0 and f(x) reaches its global minimum value f(1)=0 when x=1. Next, let us define the following Lyapunov function

    U(t)=U1(t)+U12(t)+ni=2Ui2(t)+U3(t)+U14(t)+ni=2Ui4(t)+U15(t)+ni=2Ui5(t), (40)

    where

    U1(t)=1q1(0)0m1h(τ)π11(τ)dτf(SvS1v),U12(t)=1S1vq1(0)0m1h(τ)π11(τ)dτI1vf(I1vI1v),Ui2(t)=1S1vq1(0)0m1h(τ)π11(τ)dτIiv,U3(t)=S1hf(ShS1h),U14(t)=1R100p1(τ)E1h(τ)f(E1h(τ,t)E1h(τ))dτ,Ui4(t)=1Δiq1(0)0m1h(τ)π11(τ)dτ0pi(τ)Eih(τ,t)dτ,U15(t)=1q1(0)0m1h(τ)π11(τ)dτ0q1(a)I1h(a)f(I1h(a,t)I1h(a))da.Ui5(t)=1q1(0)0m1h(τ)π11(τ)dτ0qi(a)Iih(a,t)da, (41)

    and

    qj(a)=aβjv(s)esa(μh+αjh(σ)+rjh(σ))dσds,pj(τ)=Δjqj(0)τmjh(s)esτ(μh+mjh(σ))dσds. (42)

    Direct computation gives

    pj(0)=Rj0,

    and

    qj(a)=βjv(a)+(μh+αjh(a)+rjh(a))qj(a),pj(τ)=Δjqj(0)mjh(τ)+(μh+mjh(τ))pj(τ). (43)

    The main difficulty with the Lyapunov function U above is that the Lyapunov function U is well defined. Thus in the following we first show that strain one persists both in the hosts and in the vectors as the other strains die out. Let

    ˆX1={φ1L1+(0,)|s0: 0m1h(τ+s)φ1(τ)dτ>0},
    ˆX2={ψ1L1+(0,)|s0: 0β1v(a+s)ψ1(a)da>0},

    and define

    X0=R+×nj=1R+×R+׈X1׈X2×ni=2(L1(0,)×L1(0,)),
    Ω0=ΩX0.

    Note that Ω0 is forward invariant. This is because (3) show that Ω is forward invariant. To see X0 is forward invariant, we firstly demonstrate that ˆX2 is forward invariant. Let us assume that the inequality holds for the initial condition. The inequality says that the support of β1v(a) will intersect the support of the initial condition if it is transferred s units to the right. Since the support of the initial condition only moves to the right, the intersection will take place for any other time if that happens for the initial time. Similarly, \hat{X }_1 is also forward invariant. Therefore, \Omega_0 is forward invariant.

    Now let us recall two important definitions.

    Definition 5.1. Strain one is called uniformly weakly persistence if there exists some \gamma>0 independent of the initial conditions such that

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

    and

    \displaystyle \limsup\limits_{t\rightarrow\infty}I_v^1(t) > \gamma \text{whenever} I_{v_0}^1> 0,

    for all solutions of system (2).

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

    Definition 5.2. Strain one is uniformly strongly persistence if there exists some \gamma>0 independent of the initial conditions such that

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

    and

    \displaystyle \liminf\limits_{t\rightarrow\infty}I_v^1(t) > \gamma \text{whenever} I_{v_0}^1 > 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.

    Now we are able to state the main results in this section.

    Theorem 5.3. Assume \mathcal R^1_0 >1 and \mathcal R^i_0 <\mathcal R^1_0 for i=2, \cdots, n. Furthermore, assume that the other strains except stain 1 will die out, i.e.,

    \displaystyle\limsup\limits_{t\rightarrow +\infty}I_v^i(t)=0,\ \limsup\limits_{t\rightarrow +\infty}\int^\infty_0 E_h^i(\tau, t)d\tau=0 \ \;and\;\ \limsup\limits_{t\rightarrow +\infty}\int^\infty_0I_h^i(a, t)da=0,

    for i=2, \cdots, n. Then strain 1 is uniformly weakly persistent for the initial conditions that belong to \Omega_0, i.e., there exists \gamma>0 such that

    \displaystyle \limsup\limits_{t\rightarrow+\infty}\beta_h^1I_v^1(t)\geq\gamma, \limsup\limits_{t\rightarrow+\infty}\int^\infty_0 m_h^1(\tau)E_h^1(\tau, t)d\tau\geq\gamma, \limsup\limits_{t\rightarrow+\infty}\int^\infty_0 \beta_v^1(a)I_h^1(a, t)da\geq\gamma. }

    Proof. We argue by contradiction. Assume that strain 1 also dies out. For any \varepsilon>0 and every initial condition in \Omega_0 such that

    \displaystyle \limsup\limits_{t\rightarrow+\infty}\beta_h^1I_v^1(t) < \varepsilon, \limsup\limits_{t\rightarrow+\infty}\int^\infty_0 m_h^1(\tau)E_h^1(\tau, t)d\tau < \varepsilon, \limsup\limits_{t\rightarrow+\infty}\int^\infty_0 \beta_v^1(a)I_h^1(a, t)da<\varepsilon.

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

    \displaystyle \beta_h^jI_v^j(t)<\varepsilon,\ \int^\infty_0 m_h^j(\tau)E_h^j(\tau, t)d\tau<\varepsilon,\ \int^\infty_0 \beta_v^j(a)I_h^j(a, t)da<\varepsilon,\ j=1,\cdots,n.

    We may assume that the above inequality holds for all t\geq0 by shifting the dynamical system. From the first equation in (2) we have

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

    Exploiting the comparison principle, we have

    \displaystyle \limsup\limits_{t\rightarrow+\infty}S_v(t)\geq \liminf\limits_{t\rightarrow+\infty}S_v(t)\geq\frac{\Lambda_v}{n\varepsilon+\mu_v},\ \limsup\limits_{t\rightarrow+\infty}S_h(t)\geq \liminf\limits_{t\rightarrow+\infty}S_h(t)\geq\frac{\Lambda_h}{n\varepsilon+\mu_h}.

    Since B_ E^1(t)=E_h^1(0, t),\ B_ I^1(t)=I_h^1(0, t) , it then follows from system (2) that

    \left\{ \begin{array}{ll} \displaystyle B_ E^1(t)=E_h^1(0, t)=\beta_h^1S_hI_v^1(t)\geq\beta_h^1\frac{\Lambda_h}{n\varepsilon+\mu_h}I_v^1(t),\\[2ex] \displaystyle \frac{dI_v^1(t)}{dt}\geq \frac{\Lambda_v}{n\varepsilon+\mu_v}\int^\infty_0 \beta_v^1(a)I_h^1(a, t)da-(\mu_v+\alpha_v^1)I_v^1(t). \end{array}\right. (44)

    By using the equations in (7), we can easily obtain the following inequalities on B_ E^1(t),\ B_ I^1(t) and I_v^1(t):

    \left\{ \begin{array}{ll} \displaystyle B_ E^1(t)\geq\beta_h^1\frac{\Lambda_h}{n\varepsilon+\mu_h}I_v^1(t),\\[2ex] \displaystyle B_ I^1(t)=\int^\infty_0 m_h^1(\tau)E_h^1(\tau, t)d\tau\geq\int^t_0 m_h^1(\tau)B_ E^1( t-\tau)\pi^1_1(\tau)d\tau,\\[2ex] \displaystyle \frac{dI_v^1(t)}{dt}\geq \frac{\Lambda_v}{n\varepsilon+\mu_v}\int^t_0 \beta_v^1(a)B_ I^1( t-a)\pi^1_2(a)da-(\mu_v+\alpha_v^1)I_v^1(t). \end{array}\right. (45)

    Let us take the Laplace transform of both sides of inequalities (45). Since all functions above are bounded, the Laplace transforms of the functions exist for \lambda>0. Denote the Laplace transforms of the functions B_E^1(t), B_I^1(t) and I_v^1(t) by \hat{B}_E^1(\lambda), \hat{B}_I^1(\lambda) and \hat{I}_v^1(\lambda), respectively. Furthermore, set

    \displaystyle \hat{K}_1(\lambda)=\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)e^{-\lambda\tau}d\tau, \displaystyle \hat{K}_2(\lambda)=\int^\infty_0\beta_v^1(a)\pi^1_2(a)e^{-\lambda a}da. (46)

    Using the convolution property of the Laplace transform, we obtain the following inequalities for \hat{B}_E^1(\lambda),\ \hat{B}_I^1(\lambda) and \hat{I}_v^1(\lambda):

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

    Eliminating \hat{B}_I^1(\lambda) and \hat{I}_v^1(\lambda) yields

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

    This is impossible since

    \displaystyle \frac{\beta_h^1\Lambda_v\Lambda_h\hat{K}_1(0)\hat{K}_2(0)} {\mu_v\mu_h(\mu_v+\alpha_v^1)}:=\mathcal R^1_0>1,

    we can choose \varepsilon and \lambda small enough such that

    \displaystyle \frac{\beta_h^1\Lambda_v\Lambda_h\hat{K}_1(\lambda)\hat{K}_2(\lambda)} {(n\varepsilon+\mu_v)(n\varepsilon+\mu_h)(\lambda+\mu_v+\alpha_v^1)}>1.

    The contradiction implies that there exists \gamma>0 such that for any initial condition in \Omega_0, we have

    \displaystyle \limsup\limits_{t\rightarrow+\infty}\beta_h^1I_v^1(t)\geq\gamma, \limsup\limits_{t\rightarrow+\infty}\int^\infty_0 m_h^1(\tau)E_h^1(\tau, t)d\tau\geq\gamma, \limsup\limits_{t\rightarrow+\infty}\int^\infty_0 \beta_v^1(a)I_h^1(a, t)da\geq\gamma.

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

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

    which implies a lower bound for I_v^1. This concludes the proof.

    Next, we claim that system (2) has a global compact attractor \mathfrak{T}. Firstly, define the semiflow \Psi: [0, \infty)\times \Omega_0\rightarrow \Omega_0 generated by the solutions of system (2)

    \begin{array}{ll} \displaystyle \Psi\bigg(t; S_{v_0}, I_{v_0}^1,\cdots,I_{v_0}^n, S_{h_0}, \varphi_1(\cdot), \psi_1(\cdot),\cdots,\varphi_n(\cdot), \psi_n(\cdot)\bigg)\\[2ex] =\displaystyle \bigg( S_v(t), I_v^1(t),\cdots,I_v^n(t),S_h(t), E_h^1(\tau, t),I_h^1(a,t), \cdots, E_h^n(\tau, t),I_h^n(a,t)\bigg). \end{array}

    Definition 5.4. A set \mathfrak{T} in \Omega_0 is called a global compact attractor for \Psi if \mathfrak{T} is a maximal compact invariant set and for all open sets \mathfrak{U} containing \mathfrak{T} and all bounded sets \mathcal{B} of \Omega_0 there exists some T>0 such that \Psi(t, \mathcal{B})\subseteq\mathfrak{U} holds for t>T.

    Theorem 5.5. Under the hypothesis of Theorem 5.3, there exists \mathfrak{T}, a compact subset of \Omega_0, which is a global attractor for the semiflow \Psi on \Omega_0. Moreover, we have

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

    Proof. We split the solution semiflow into two components. For an initial condition x^0\in \Omega_0, }let \Psi(t, x^0)=\hat{\Psi}(t, x^0)+\tilde{\Psi}(t, x^0), where

    \begin{array}{ll} \displaystyle \hat{\Psi}\bigg(t; S_{v_0}, I_{v_0}^1,\cdots,I_{v_0}^n, S_{h_0}, \varphi_1(\cdot), \psi_1(\cdot),\cdots,\varphi_n(\cdot), \psi_n(\cdot)\bigg)\\[2ex] =\displaystyle \bigg( 0, 0,\cdots,0,0, \hat{E}_h^1(\tau, t),\hat{I}_h^1(a,t), \cdots, \hat{E}_h^n(\tau, t),\hat{I}_h^n(a,t)\bigg), \end{array} (49)
    \begin{array}{ll} \displaystyle \tilde{\Psi}\bigg(t; S_{v_0}, I_{v_0}^1,\cdots,I_{v_0}^n, S_{h_0}, \varphi_1(\cdot), \psi_1(\cdot),\cdots,\varphi_n(\cdot), \psi_n(\cdot)\bigg)\\[2ex] =\displaystyle \bigg( S_v(t), I_v^1(t),\cdots,I_v^n(t),S_h(t), \tilde{E}_h^1(\tau, t),\tilde{I}_h^1(a,t), \cdots, \tilde{E}_h^n(\tau, t),\tilde{I}_h^n(a,t)\bigg), \end{array} (50)

    and \displaystyle E_h^j(\tau, t)=\hat{E}_h^j(\tau, t)+\tilde{E}_h^j(\tau, t),\ I_h^j(a, t)=\hat{I}_h^j(a, t)+\tilde{I}_h^j(a, t) for j=1,\cdots,n. \hat{E}_h^j(\tau, t), \hat{I}_h^j(a, t), \tilde{E}_h^j(\tau, t), \tilde{I}_h^j (a, t) are the solutions of the following equations

    \left\{ \begin{array}{ll} \displaystyle \frac{\partial \hat{E}_h^j}{\partial t}+\frac{\partial \hat{E}_h^j}{\partial \tau}=-(\mu_h+m_h^j(\tau))\hat{E}_h^j(\tau, t),\\[2ex] \displaystyle \hat{E}_h^j(0, t)=0,\\[2ex] \displaystyle \hat{E}_h^j(\tau, 0)=\varphi_j(\tau), \end{array}\right. (51)
    \left\{ \begin{array}{ll} \displaystyle \frac{\partial \hat{I}_h^j}{\partial t}+\frac{\partial \hat{I}_h^j}{\partial a}=-(\mu_h+\alpha_h^j(a)+r_h^j(a))\hat{I}_h^j(a, t),\\[2ex] \displaystyle \hat{I}_h^j(0, t)=0,\\[2ex] \displaystyle \hat{I}_h^j(a, 0)=\psi_j(a), \end{array}\right. (52)

    and

    \left\{ \begin{array}{ll} \displaystyle \frac{\partial \tilde{E}_h^j}{\partial t}+\frac{\partial \tilde{E}_h^j}{\partial \tau}=-(\mu_h+m_h^j(\tau))\tilde{E}_h^j(\tau, t),\\[2ex] \displaystyle \tilde{E}_h^j(0, t)=\beta_h^jS_hI_v^j,\\ \displaystyle \tilde{E}_h^j(\tau, 0)=0, \end{array}\right. (53)
    \left\{ \begin{array}{ll} \displaystyle \frac{\partial \tilde{I}_h^j}{\partial t}+\frac{\partial \tilde{I}_h^j}{\partial a}=-(\mu_h+\alpha_h^j(a)+r_h^j(a))\tilde{I}_h^j(a, t),\\[2ex] \displaystyle \tilde{I}_h^j(0, t)=\int^\infty_0m_h^j(\tau)\tilde{E}_h^j(\tau, t)d\tau,\\[2ex] \displaystyle \tilde{I}_h^j(a, 0)=0. \end{array}\right. (54)

    We can easily see that system (51) and (52) are decoupled from the remaining equations. Using the formula (7) to integrate along the characteristic lines, we obtain

    \begin{array}{ll} \displaystyle \hat{E}_h^j(\tau, t)=\left\{\begin{array}{ll} \displaystyle 0,\mbox{ } & t>\tau,\\[2ex] \displaystyle \varphi_j(\tau-t)\frac{\pi^j_1(\tau)}{\pi^j_1(\tau-t)},\mbox{ } & t<\tau, \end{array}\right. \end{array} (55)
    \begin{array}{ll} \displaystyle \hat{I}_h^j(a, t)=\left\{\begin{array}{ll} \displaystyle 0,\mbox{ } & t>a,\\[2ex] \displaystyle \psi_j(a-t)\frac{\pi^j_2(a)}{\pi^j_2(a-t)},\mbox{ } & t<a. \end{array}\right.\\[2ex] \end{array} (56)

    Integrating \hat{E}_h^j with respect to \tau yields

    \displaystyle \int^\infty_t \varphi_j(\tau-t)\frac{\pi^j_1(\tau)}{\pi^j_1(\tau-t)}d\tau=\int^\infty_0 \varphi_j(\tau)\frac{\pi^j_1(t+\tau)}{\pi^j_1(\tau)}d\tau\leq e^{-\mu_h t}\int^\infty_0 \varphi_j(\tau)d\tau\rightarrow 0

    as t\rightarrow\infty. Integrating \hat{I}_h^j with respect to a, we have

    \displaystyle \int^\infty_t \psi_j(a-t)\frac{\pi^j_2(a)}{\pi^j_2(a-t)}da=\int^\infty_0 \psi_j(a)\frac{\pi^j_2(t+a)}{\pi^j_2(a)}da\leq e^{-\mu_h t}\int^\infty_0 \psi_j(a)da \rightarrow 0

    as t\rightarrow\infty. This implies that \hat{\Psi}(t, x^0)\rightarrow0 as t\rightarrow\infty uniformly for every x^0\in\mathcal{B}\subseteq \Omega_0, where \mathcal{B} is a ball of a given radius.

    In the following we need to show \tilde{\Psi}(t, x) is completely continuous. We fix t and let x^0\in \Omega_0. Note that \Omega_0 is bounded. We have to show that the family of functions defined by

    \begin{array}{ll} \displaystyle \tilde{\Psi}(t, x^0)=\displaystyle \bigg( S_v(t), I_v^1(t),\cdots,I_v^n(t),S_h(t), \tilde{E}_h^1(\tau, t),\tilde{I}_h^1(a,t), \cdots, \tilde{E}_h^n(\tau, t),\tilde{I}_h^n(a,t)\bigg) \end{array}

    is a compact family of functions for that fixed t, which are obtained by taking different initial conditions in \Omega_0. The family

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

    and, therefore, it is bounded. Thus, we have established the boundedness of the set. To show that \tilde{\Psi}(t, x) is precompact, we first see the third condition of \lim_{t\rightarrow\infty}\int^\infty_t\tilde{E}_h^j(\tau, t)d\tau=0 and \lim_{t\rightarrow\infty}\int^\infty_t\tilde{I}_h^j(a, t)da=0 in the Frechet-Kolmogorov Theorem of [21]. The third condition in [21] is trivially satisfied since \tilde{E}_h^j(\tau, t)=0 for \tau>t and \tilde{I}_h^j(a, t)=0 for a>t. To use the second condition of the Frechet-Kolmogorov Theorem in [21], we must bound by two constants the L^1-norms of \partial E_h^j/\partial\tau and \partial I_h^j/\partial a. Notice that

    \begin{array}{ll} \displaystyle \tilde{E}_h^j(\tau, t)=\left\{\begin{array}{ll} \displaystyle \tilde{B}_E^j(t-\tau)\pi^j_1(\tau),\mbox{ } & t>\tau,\\[2ex] \displaystyle 0,\mbox{ } & t<\tau,\\[2ex] \end{array}\right.\\[2ex] \displaystyle \tilde{I}_h^j(a, t)=\left\{\begin{array}{ll} \displaystyle \tilde{B}_I^j(t-a)\pi^j_2(a),\mbox{ } & t>a,\\[2ex] \displaystyle 0,\mbox{ } & t<a, \end{array}\right. \end{array} (57)

    where

    \begin{array}{ll} \displaystyle \tilde{B}_E^j(t)=\displaystyle \beta_h^jS_h(t)I_v^j(t),\\[2ex] \displaystyle \tilde{B}_I^j(t)=\displaystyle \int^\infty_0m_h^j(\tau)\tilde{E}_h^j(\tau, t)d\tau=\int^t_0m_h^j(\tau)\tilde{B}_E^j(t-\tau)\pi^j_1(\tau)d\tau.\\[2ex] \end{array} (58)

    \tilde{B}_E^j(t) is bounded because of the boundedness of S_h and I_v^j. Hence, the \tilde{B}_E^j(t) satisfies

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

    Therefore, we obtain

    \begin{array}{ll} \displaystyle \tilde{B}_I^j(t) & \displaystyle=\int^t_0m_h^j(\tau)\tilde{B}_E^j(t-\tau)\pi^j_1(\tau)d\tau\\[2ex] & \displaystyle\leq k_2\int^t_0\tilde{B}_E^j(t-\tau)d\tau= k_2\int^t_0\tilde{B}_E^j(\tau)d\tau\\[2ex] & \displaystyle \leq k_1k_2t. \end{array}

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

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

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

    \begin{array}{ll} \displaystyle (\tilde{B}_E^j(t))'=\displaystyle \beta_h^j\bigg(S'_h(t)I_v^j(t)+S_h(t)(I_v^j(t))'\bigg),\\[2ex] (\tilde{B}_I^j(t))'=\displaystyle m_h^j(t)\tilde{B}_E^j(0)\pi^j_1(t)+\int^t_0m_h^j(\tau)(\tilde{B}_E^j(t-\tau))'\pi^j_1(\tau)d\tau. \end{array} (59)

    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^j(t))'|\leq k_3, |(\tilde{B}_I^j(t))'|\leq k_4. \end{array}

    Putting all these bounds together, we have

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

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

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

    Thus, the integral can be made arbitrary small uniformly in the family of functions. That establishes the second condition of the Frechet-Kolmogorov Theorem. We conclude that the family is asymptotically smooth.

    (3) means that the semigroup \Psi(t) is point dissipative and the forward orbit of boundedness sets is bounded in \Omega_0. Thus, we prove Theorem 5.5 in accordance with Lemma 3.1.3 and Theorem 3.4.6 in [8].

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

    Theorem 5.6. Under the hypothesis of Theorem 5.3 strain one is uniformly strongly persistent for all initial conditions that belong to \Omega_0, that is, there exists \gamma>0 such that

    \displaystyle \liminf\limits_{t\rightarrow+\infty}\beta_h^1I_v^1(t)\geq\gamma, \liminf\limits_{t\rightarrow+\infty}\int^\infty_0 m_h^1(\tau)E_h^1(\tau, t)d\tau\geq\gamma, \liminf\limits_{t\rightarrow+\infty}\int^\infty_0 \beta_v^1(a)I_h^1(a, t)da\geq\gamma. }

    Proof. We apply Theorem 2.6 in [19]. We consider the solution semiflow \Psi on \Omega_0. Let us define three functionals \rho_l: \Omega_0\rightarrowR_+,\ l=1,2,3 as follows:

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

    Theorem 5.3 implies that the semiflow is uniformly weakly \rho-persistent. Theorem 5.5 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

    \displaystyle \beta_h^1I_v^1(t)= \beta_h^1I_v^1(s)e^{-(\mu_v+\alpha_v^1)(t-s)},
    \begin{array}{rl} \displaystyle \int^\infty_0 m_h^1(\tau)\tilde{E}_h^1(\tau, t)d\tau & \displaystyle =\tilde{B}_I^1(t)=\int^t_0m_h^1(\tau)\tilde{B}_E^1(t-\tau)\pi^1_1(\tau)d\tau \\[2ex] & \displaystyle\geq k^1\int^t_0\tilde{B}_E^1(t-\tau)d\tau= k^1\int^t_0\tilde{B}_E^1(\tau)d\tau\\[2ex] & \displaystyle=k^1\int^t_0\beta_h^1S_h(\tau)I_v^1(\tau)d\tau\geq k^2\int^t_0I_v^1(\tau)d\tau \end{array}
    \begin{array}{rl} & \displaystyle=k^2\int^t_0I_v^1(s)e^{-(\mu_v+\alpha_v^1)(\tau-s)}d\tau\\[2ex] & \displaystyle=\frac{k^2I_v^1(s)}{\mu_v+\alpha_v^1}e^{(\mu_v+\alpha_v^1)s}(1-e^{-(\mu_v+\alpha_v^1)t}),\\[2ex] \displaystyle \int^\infty_0 \beta_v^1(a)\tilde{I}_h^1(a, t)da & \displaystyle=\int^t_0 \beta_v^1(a)\tilde{B}_I^1(t-a)\pi^1_2(a) da\geq k^3\int^t_0\tilde{B}_I^1(t-a)da\\[2ex] & \displaystyle= k^3\int^t_0\tilde{B}_I^1(a)da\\[2ex] & \displaystyle\geq\frac{k^2k^3I_v^1(s)}{\mu_v+\alpha_v^1}e^{(\mu_v+\alpha_v^1)s}\int^t_0(1-e^{-(\mu_v+\alpha_v^1)a})da, \end{array}

    for any s and any t>s. Therefore,

    \beta_h^1I_v^1(t)>0, \int^\infty_0 m_h^1(\tau)\tilde{E}_h^1(\tau, t)d\tau>0, \int^\infty_0\beta_v^1(a)\tilde{I}_h^1(a, t)da>0

    for all t>s provided \tilde{I}_v^1(s)>0. Theorem 2.6 in [19] now implies that the semiflow is uniformly strongly \rho-persistent. Hence, there exists \gamma such that

    \displaystyle \liminf\limits_{t \rightarrow +\infty}\beta_h^1I_v^1(t)\geq\gamma, \liminf\limits_{t \rightarrow +\infty}\int^\infty_0 m_h^1(\tau)E_h^1(\tau, t)d\tau\geq\gamma, \liminf\limits_{t \rightarrow +\infty}\int^\infty_0 \beta_v^1(a)I_h^1(a, t)da\geq\gamma.

    According to Theorem 5.6, we obtain that for all initial conditions that belong to \Omega_0, strain 1 persists. Furthermore we had verified that the solutions of (2) with nonnegative initial conditions belong to the positive cone for all t\geq0. All the solutions are in a positively invariant set. Therefore we can obtain the following Theorem 5.7 from Theorem 5.6.

    Theorem 5.7. Under the hypothesis of Theorem 5.3, \forall t\in\mathrm{R}, there exists constants \vartheta>0 and M>0 such that

    \vartheta\leq S_v(t)\leq M, \vartheta\leq S_h(t)\leq M,

    and

    \displaystyle\vartheta\leq \beta_h^1I_v^1(t)\leq M,\ \displaystyle\vartheta\leq \int^\infty_0 m_h^1(\tau)E_h^1(\tau, t)d\tau\leq M,\ \displaystyle\vartheta\leq \int^\infty_0 \beta_v^1(a)I_h^1(a, t)da\leq M,

    for each orbit ( S_v(t), I_v^1(t),\cdots, I_v^n(t), S_h(t), E_h^1(\tau, t), I_h^1(a,t),\cdots,E_h^n(\tau, t), I_h^n(a,t)) of \Psi in \mathfrak{T}.


    6. Principle of competitive exclusion

    In this section we mainly state the main result of the paper.

    Theorem 6.1. Assume \mathcal R^1_0>1,\mathcal R^i_0/\mathcal R^1_0<b_i/b_1<1,i=2,\cdots,n. Then the equilibrium \mathcal{E}_1 is globally asymptotically stable.

    Proof. From Theorem 4.1 we know that the endemic equilibrium \mathcal{E}_1 is locally asymptotically stable. In the following we only need to show that the endemic equilibrium \mathcal{E}_1 is global attractor. From Theorem 5.5 there exists an invariant compact set \mathfrak{T} which is global attractor of system (2). Furthermore, it follows from Theorem 5.7 that there exist \varepsilon_1>0 and M_1>0 such that

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

    for any solution in \Psi. This makes the Lyapunov function defined in (40) well defined.

    After extensive computation, differentiating U(t) along the solution of system (2) yields that

    \begin{split} & \displaystyle \frac{dU_1(t)}{dt}\\ = & \displaystyle \frac{1}{S_v^{1^*}q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\bigg(1-\frac{S_v^{1^*}}{S_v }\bigg) \bigg[S_v^{1^*}\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)da+\mu_vS_v^{1^*}\\ & \displaystyle-S_v\int^\infty_0\beta_v^1(a)I_h^1(a, t)da-\mu_vS_v-\sum^n_{i=2}S_v\int^\infty_0\beta_v^i(a)I_h^i(a, t)da\bigg]\\ = & \displaystyle -\frac{\mu_v(S_v-S_v^{1^*})^2}{S_v^{1^*}S_vq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\\ & \displaystyle+\frac{1}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\int^\infty_0 \beta_v^1(a)I_h^{1^*}(a)\bigg(1-\frac{S_v^{1^*}}{S_v}-\frac{S_vI_h^1(a, t)} {S_v^{1^*}I_h^{1^*}(a)}+\frac{I_h^1(a, t)}{I_h^{1^*}(a)}\bigg)da\\ & \displaystyle-\sum^n_{i=2} \frac{S_v\int^\infty_0\beta_v^i(a)I_h^i(a, t)da-S_v^{1^*}\int^\infty_0\beta_v^i(a)I_h^i(a, t)da}{S_v^{1^*}q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}; \end{split} (60)
    \begin{array}{ll} & \displaystyle \frac{dU^1_2(t)}{dt}\\[3ex] = & \displaystyle \frac{(1-\frac{I_v^{1^*}}{I_v^1})\bigg(S_v\int^\infty_0\beta_v^1(a)I_h^1(a, t)da-\frac{S_v^{1^*}\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)da}{I_v^{1^*}}I_v^1\bigg)}{S_v^{1^*}q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \\[3ex] = & \displaystyle \frac{(1-\frac{I_v^{1^*}}{I_v^1})S_v^{1^*}\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)\bigg(\frac{S_vI_h^1(a, t)} {S_v^{1^*}I_h^{1^*}(a)}-\frac{I_v^1}{I_v^{1^*}}\bigg)da}{S_v^{1^*}q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \\[3ex] = & \displaystyle \frac{\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)\bigg(\frac{S_vI_h^1(a, t)} {S_v^{1^*}I_h^{1^*}(a)}-\frac{I_v^1}{I_v^{1^*}}-\frac{S_vI_h^1(a, t)I_v^{1^*}} {S_v^{1^*}I_h^{1^*}(a)I_v^1}+1\bigg)da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}; \end{array} (61)
    \begin{array}{ll} \displaystyle \frac{dU^i_2(t)}{dt}=\frac{S_v\int^\infty_0\beta_v^i(a)I_h^i(a, t)da-(\mu_v+\alpha_v^i)I_v^i}{S_v^{1^*}q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}; \end{array} (62)
    \begin{array}{ll} \displaystyle \frac{dU_3(t) }{dt} & =\displaystyle \bigg(1-\frac{S_h^{1^*}}{S_h}\bigg)\bigg(E_h^{1^*}(0)+\mu_hS_h^{1^*}-E_h^1(0,t)-\mu_hS_h-\sum^n_{i=2}\beta_h^iS_hI_v^i\bigg)\\[2ex] & =\displaystyle -\frac{\mu_h(S_h-S_h^{1^*})^2}{S_h} \displaystyle+\bigg(E_h^{1^*}(0)-E_h^1(0,t)-\frac{S_h^{1^*}}{S_h}E_h^{1^*}(0)+\frac{S_h^{1^*}}{S_h}E_h^1(0,t)\bigg)\\[2ex] & \displaystyle-\sum^n_{i=2}\bigg(E_h^i(0,t)-\beta_h^iS_h^{1^*}I_v^i\bigg),\\[2ex] \end{array} (63)

    and

    \begin{array}{ll} & \displaystyle \frac{dU^1_4(t)}{dt}\\[2ex] = & \displaystyle \frac{1}{\mathcal R^1_0}\int^\infty_0p_1(\tau)E_h^{1^*}(\tau)f'\bigg(\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)}\bigg)\frac{1}{E_h^{1^*}(\tau)}\frac{\partial E_h^1(\tau, t)}{\partial t}d\tau\\[2ex] = & \displaystyle -\frac{1}{\mathcal R^1_0}\int^\infty_0p_1(\tau)E_h^{1^*}(\tau)f'\bigg(\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)}\bigg)\frac{1}{E_h^{1^*}(\tau)}\bigg(\frac{\partial E_h^1(\tau, t)}{\partial \tau}+(\mu_h+m_h^1(\tau))E_h^1(\tau, t)\bigg)d\tau\\[2ex] = & \displaystyle -\frac{1}{\mathcal R^1_0}\int^\infty_0p_1(\tau)E_h^{1^*}(\tau)df\bigg(\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)}\bigg)\\[2ex] = & \displaystyle -\frac{1}{\mathcal R^1_0}\bigg[p_1(\tau)E_h^{1^*}(\tau)f\bigg(\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)}\bigg)\bigg|^{\infty}_0-\int^\infty_0f\bigg(\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)}\bigg)d\bigg(p_1(\tau)E_h^{1^*}(\tau)\bigg)\bigg]\\[2ex] = & \displaystyle \frac{1}{\mathcal R^1_0}\bigg[p_1(0)E_h^{1^*}(0)f\bigg(\frac{E_h^1(0, t)}{E_h^{1^*}(0)}\bigg)-\Delta_1 q_1(0)\int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)f\bigg(\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)}\bigg)d\tau\bigg]\\[2ex] = & \displaystyle E_h^{1^*}(0)f\bigg(\frac{E_h^1(0, t)}{E_h^{1^*}(0)}\bigg)-\frac{\int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)f(\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)})d\tau}{\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\\[2ex] = & \displaystyle E_h^1(0, t)-E_h^{1^*}(0)-E_h^{1^*}(0)\ln\frac{E_h^1(0, t)}{E_h^{1^*}(0)}-\frac{\int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)f(\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)})d\tau}{\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}. \end{array} (64)

    The above equality follows from (24) and the fact

    \begin{array}{ll} & \displaystyle p_1'(\tau)E_h^{1^*}(\tau)+p_1(\tau)(E_h^{1^*}(\tau))'\\[2ex] = & \displaystyle \bigg[-\Delta_1 q_1(0)m_h^1(\tau)+(\mu_h+m_h^1(\tau))p_1(\tau)\bigg]E_h^{1^*}(\tau)-p_1(\tau)(\mu_h+m_h^1(\tau))E_h^{1^*}(\tau) \\[2ex] = & \displaystyle -\Delta_1 q_1(0)m_h^1(\tau)E_h^{1^*}(\tau). \end{array}

    We also have

    \begin{array}{ll} & \displaystyle q_1'(a)I_h^{1^*}(a)+q_1(a)(I_h^{1^*}(a))'\\[2ex] = & \displaystyle \bigg[-\beta_v^1(a)+(\mu_h+\alpha_h^1(a)+r_h^1(a))q_1(a)\bigg]I_h^{1^*}(a)-q_1(a)(\mu_h+\alpha_h^1(a)+r_h^1(a))I_h^{1^*}(a) \\[2ex] = & \displaystyle -\beta_v^1(a)I_h^{1^*}(a). \end{array}

    Similar to the differentiation of U^1_4(t), we have

    \begin{array}{ll} & \displaystyle \frac{dU^1_5(t)}{dt}\\[2ex] = & \displaystyle\frac{1}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\displaystyle \int^\infty_0q_1(a)I_h^{1^*}(a)f'\bigg(\frac{I_h^1(a, t)}{I_h^{1^*}(a)}\bigg)\frac{1}{I_h^{1^*}(a)}\frac{\partial I_h^1(a, t)}{\partial t}da\\[2ex] = & \displaystyle \frac{1}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\displaystyle \int^\infty_0q_1(a)I_h^{1^*}(a)df\bigg(\frac{I_h^1(a, t)}{I_h^{1^*}(a)}\bigg)\\[2ex] = & \displaystyle \frac{q_1(0)I_h^{1^*}(0)f(\frac{I_h^1(0, t)}{I_h^{1^*}(0)})-\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)f(\frac{I_h^1(a, t)}{I_h^{1^*}(a)})da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\\[2ex] = & \displaystyle \frac{\int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)(\frac{I_h^1(0, t)}{I_h^{1^*}(0)}-1-\ln\frac{I_h^1(0, t)}{I_h^{1^*}(0)})d\tau}{\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}-\frac{\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)f(\frac{I_h^1(a, t)}{I_h^{1^*}(a)})da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}.\\[2ex] \end{array} (65)

    Noting that (43), we differentiate the last two terms with respect to t, and have

    \begin{array}{ll} & \displaystyle \frac{dU^i_4(t)}{dt}\\[2ex] = & \displaystyle \frac{1}{\Delta_iq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\int^\infty_0p_i(\tau)\frac{\partial E_h^i(\tau,t)}{\partial t}d\tau \\[2ex] = & \displaystyle \displaystyle -\frac{1}{\Delta_iq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\int^\infty_0p_i(\tau)\bigg[\frac{\partial E_h^i(\tau,t)}{\partial \tau}+(\mu_h+m_h^i(\tau))E_h^i(\tau,t)\bigg]d\tau \\[2ex] = & \displaystyle \displaystyle -\frac{\int^\infty_0p_i(\tau)d E_h^i(\tau,t) +\int^\infty_0(\mu_h+m_h^i(\tau))p_i(\tau)E_h^i(\tau,t)d\tau}{\Delta_iq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \\[2ex] = & \displaystyle -\frac{p_i(\tau) E_h^i(\tau,t)|^\infty_0 -\int^\infty_0E_h^i(\tau,t)d p_i(\tau) +\int^\infty_0(\mu_h+m_h^i(\tau))p_i(\tau)E_h^i(\tau,t)d\tau}{\Delta_iq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \\[2ex] = & \displaystyle \frac{p_i(0) E_h^i(0,t) -\Delta_i q_i(0)\int^\infty_0m_h^i(\tau)E_h^i(\tau,t)d\tau}{\Delta_iq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\\[2ex] = & \displaystyle \frac{\mathcal R^i_0 E_h^i(0,t)}{\Delta_iq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} -\frac{ q_i(0)I_h^i(0,t)}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\\[2ex] = & \displaystyle\frac{b_i}{b_1}E_h^i(0,t) -\frac{ q_i(0)I_h^i(0,t)}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}. \end{array} (66)

    Similarly, we have

    \begin{array}{ll} & \displaystyle \frac{dU^i_5(t)}{dt}\\[2ex] = & -\displaystyle \frac{\int^\infty_0q_i(a)\bigg[\frac{\partial I_h^i(a,t)}{\partial a}+(\mu_h+\alpha_h^i(a)+r_h^i(a))I_h^i(a,t)\bigg]da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \\[2ex] = & \displaystyle \frac{q_i(0) I_h^i(0,t)}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} - \frac{\int^\infty_0\beta_v^i(a)I_h^i(a,t)da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}.\\[2ex] \end{array} (67)

    Adding all five components of the Lyapunov function, we have

    \displaystyle U'(t)=U^1+U^2,

    where

    \begin{array}{ll} & \displaystyle U^1(t)\\[2ex] = & \displaystyle -\frac{\mu_v(S_v-S_v^{1^*})^2}{S_v^{1^*}S_vq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\\[2ex] & \displaystyle+\frac{\int^\infty_0 \beta_v^1(a)I_h^{1^*}(a)\bigg(1-\frac{S_v^{1^*}}{S_v}-\frac{S_vI_h^1(a, t)} {S_v^{1^*}I_h^{1^*}(a)}+\frac{I_h^1(a, t)}{I_h^{1^*}(a)}\bigg)da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\\[2ex] & +\displaystyle \frac{\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)\bigg(\frac{S_vI_h^1(a, t)} {S_v^{1^*}I_h^{1^*}(a)}-\frac{I_v^1}{I_v^{1^*}}-\frac{S_vI_h^1(a, t)I_v^{1^*}} {S_v^{1^*}I_h^{1^*}(a)I_v^1}+1\bigg)da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\\[2ex] & \displaystyle -\frac{\mu_h(S_h-S_h^{1^*})^2}{S_h} \displaystyle+\bigg(E_h^{1^*}(0)-E_h^1(0,t)-\frac{S_h^{1^*}}{S_h}E_h^{1^*}(0)+\frac{S_h^{1^*}}{S_h}E_h^1(0,t)\bigg) \end{array}
    \begin{array}{ll} & \displaystyle+ E_h^1(0, t)-E_h^{1^*}(0)-E_h^{1^*}(0)\ln\frac{E_h^1(0, t)}{E_h^{1^*}(0)}-\frac{\int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)f(\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)})d\tau}{\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\\[2ex] & \displaystyle +\frac{\int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)(\frac{I_h^1(0, t)}{I_h^{1^*}(0)}-1-\ln\frac{I_h^1(0, t)}{I_h^{1^*}(0)})d\tau}{\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}-\frac{\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)f(\frac{I_h^1(a, t)}{I_h^{1^*}(a)})da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}, \end{array} (68)

    and

    \begin{array}{ll} & \displaystyle U^2(t)\\[2ex] = & \displaystyle-\sum^n_{i=2} \frac{S_v\int^\infty_0\beta_v^i(a)I_h^i(a, t)da-S_v^{1^*} \int^\infty_0\beta_v^i(a)I_h^i(a, t)da}{S_v^{1^*}q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \\[2ex] & +\displaystyle\sum^n_{i=2}\frac{S_v\int^\infty_0\beta_v^i(a)I_h^i(a, t)da-(\mu_v+\alpha_v^i)I_v^i}{S_v^{1^*}q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}-\sum^n_{i=2}\bigg(E_h^i(0,t)-\beta_h^iS_h^{1^*}I_v^i\bigg)\\[2ex] & +\displaystyle\sum^n_{i=2} \bigg(\frac{b_i}{b_1}E_h^i(0,t) -\frac{q_i(0)I_h^i(0,t)}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\bigg)\\[2ex] & + \displaystyle \sum^n_{i=2}\bigg(\frac{q_i(0) I_h^i(0,t)}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} - \frac{\int^\infty_0\beta_v^i(a)I_h^i(a,t)da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\bigg).\\[2ex] \end{array} (69)

    Canceling terms, (68) can be simplified as

    \begin{array}{ll} & \displaystyle U^1(t) \\[2ex] = & \displaystyle -\frac{\mu_v(S_v-S_v^{1^*})^2}{S_v^{1^*}S_vq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} -\frac{\mu_h(S_h-S_h^{1^*})^2}{S_h}\\[2ex] & \displaystyle+ \frac{\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)(3-\frac{S_v^{1^*}}{S_v}-\frac{I_v^1}{I_v^{1^*}}-\frac{S_vI_h^1(a, t)I_v^{1^*}} {S_v^{1^*}I_h^{1^*}(a)I_v^1}+\ln\frac{I_h^1(a, t)} {I_h^{1^*}(a)})da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\\[2ex] & \displaystyle +E_h^{1^*}(0)\bigg(-\frac{S_h^{1^*}}{S_h}+\frac{S_h^{1^*}E_h^1(0, t)}{S_hE_h^{1^*}(0)}-\ln\frac{E_h^1(0, t)}{E_h^{1^*}(0)} \bigg)\\[2ex] & \displaystyle+ \frac{\int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)(\frac{I_h^1(0, t)}{I_h^{1^*}(0)}-\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)}+\ln\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)}\frac{I_h^{1^*}(0)}{I_h^1(0, t)})d\tau}{\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}.\\[2ex] \end{array} (70)

    Direct computation yields that

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

    By using (71), (70) can be simplified as

    \begin{array}{ll} \displaystyle U^1(t) & =\displaystyle -\frac{\mu_v(S_v-S_v^{1^*})^2}{S_v^{1^*}S_vq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} -\frac{\mu_h(S_h-S_h^{1^*})^2}{S_h}\\[2ex] & \displaystyle -\frac{\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)[f(\frac{S_v^{1^*}}{S_v}) +f(\frac{I_v^1}{I_v^{1^*}})+f(\frac{S_vI_h^1(a, t)I_v^{1^*}} {S_v^{1^*}I_h^{1^*}(a)I_v^1})]da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \\[2ex] & \displaystyle -\frac{1}{\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)f\bigg(\frac{E_h^1(\tau, t)I_h^{1^*}(0)}{E_h^{1^*}(\tau)I_h^1(0, t)}\bigg)d\tau\\[2ex] & \displaystyle +E_h^{1^*}(0)\bigg[-f\bigg(\frac{S_h^{1^*}}{S_h}\bigg)+f\bigg(\frac{S_h^{1^*}E_h^1(0, t)}{S_hE_h^{1^*}(0)}\bigg)\bigg].\\[2ex] \end{array} (72)

    Noting that E_h^1(0, t)=\beta_h^1 S_h I_v^1,\ E_h^{1^*}(0)=\beta_h^1 S_h^{1^*} I_v^{1^*}, we get

    \displaystyle \frac{S_h^{1^*}E_h^1(0,t)}{S_hE_h^{1^*}(0)}=\frac{S_h^{1^*}\beta_h^1 S_h I_v^1}{S_h\beta_h^1 S_h^{1^*} I_v^{1^*}}=\frac{I_v^1}{I_v^{1^*}}. (73)

    Furthermore, from (25) and (42) we have

    \begin{array}{rl} \displaystyle\frac{\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)f(\frac{I_v^1}{I_v^{1^*}})da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} & =\displaystyle\frac{\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}f(\frac{I_v^1}{I_v^{1^*}})\\[2ex] & =\displaystyle \frac{I_h^{1^*}(0)}{\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}f(\frac{I_v^1}{I_v^{1^*}})\\[2ex] & =\displaystyle E_h^{1^*}(0)f(\frac{I_v^1}{I_v^{1^*}}). \end{array} (74)

    Finally, simplifying (72) with (73) and (74), we obtain

    \begin{array}{ll} \displaystyle U^1(t) & =\displaystyle -\frac{\mu_v(S_v-S_v^{1^*})^2}{S_v^{1^*}S_vq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} -\frac{\mu_h(S_h-S_h^{1^*})^2}{S_h}\\[2ex] & \displaystyle -\frac{\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)[f(\frac{S_v^{1^*}}{S_v}) +f(\frac{S_vI_h^1(a, t)I_v^{1^*}} {S_v^{1^*}I_h^{1^*}(a)I_v^1})]da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \\[2ex] & \displaystyle -\frac{\int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)f\bigg(\frac{E_h^1(\tau, t)I_h^{1^*}(0)}{E_h^{1^*}(\tau)I_h^1(0, t)}\bigg)d\tau}{\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}-E_h^{1^*}(0)f\bigg(\frac{S_h^{1^*}}{S_h}\bigg) . \end{array} (75)

    Canceling terms, (69) can be simplified as

    \begin{array}{ll} & \displaystyle U^2(t)\\[2ex] = & \displaystyle\sum^n_{i=2}\bigg[ \bigg(\frac{b_i}{b_1}-1\bigg)E_h^i(0,t)+ \bigg(\beta_h^iS_h^{1^*}-\frac{\mu_v+\alpha_v^i}{S_v^{1^*}q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\bigg)I_v^i\bigg]. \end{array} (76)

    Simplifying (76) with (25), we get

    \begin{array}{ll} \displaystyle U^2(t)= & \displaystyle \sum^n_{i=2}\bigg[ \bigg(\frac{b_i}{b_1}-1\bigg)E_h^i(0,t) +\frac{\beta_h^i\Lambda_h(\mu_v+\Lambda_hb_1)}{\mu_h(\mu_v\mathcal R^1_0+\Lambda_hb_1)} \bigg(1-\frac{\mathcal R^1_0b_i}{\mathcal R^i_0b_1}\bigg)I_v^i\bigg]. \end{array} (77)

    Hence, by using (75) and (77) we obtain

    \begin{array}{ll} \displaystyle U'(t) & =\displaystyle -\frac{\mu_v(S_v-S_v^{1^*})^2}{S_v^{1^*}S_vq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} -\frac{\mu_h(S_h-S_h^{1^*})^2}{S_h}\\[2ex] & \displaystyle -\frac{\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)[f(\frac{S_v^{1^*}}{S_v}) +f(\frac{S_vI_h^1(a, t)I_v^{1^*}} {S_v^{1^*}I_h^{1^*}(a)I_v^1})]da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \\[2ex] & \displaystyle -\frac{\int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)f\bigg(\frac{E_h^1(\tau, t)I_h^{1^*}(0)}{E_h^{1^*}(\tau)I_h^1(0, t)}\bigg)d\tau}{\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}-E_h^{1^*}(0)f\bigg(\frac{S_h^{1^*}}{S_h}\bigg) \\[2ex] & \displaystyle+\sum^n_{i=2}\bigg[ \bigg(\frac{b_i}{b_1}-1\bigg)E_h^i(0,t) +\frac{\beta_h^i\Lambda_h(\mu_v+\Lambda_hb_1)}{\mu_h(\mu_v\mathcal R^1_0+\Lambda_hb_1)} \bigg(1-\frac{\mathcal R^1_0b_i}{\mathcal R^i_0b_1}\bigg)I_v^i\bigg]. \end{array} (78)

    Since f(x)\geq0 for x>0, \mathcal R^i_0/\mathcal R^1_0<b_i/b_1<1,i\neq1 we have U'\leq0. Define,

    \displaystyle \Theta_2= \bigg\{(S_v, I_v^1,\cdots,I_v^n, S_h, E_h^1, I_h^1,\cdots,E_h^n, I_h^n)\in \Omega_0\bigg|U'(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 (78) occurs if and only if S_v( t)=S_v^{1^*},\ S_h( t)=S_h^{1^*},\ E_h^i(0,t)=0,\ I_v^i=0, and

    \displaystyle \frac{I_h^1(a, t)I_v^{1^*}} {I_h^{1^*}(a)I_v^1}=1, \frac{E_h^1(\tau, t)I_h^{1^*}(0)}{E_h^{1^*}(\tau)I_h^1(0, t)}=1. (79)

    Thus, we obtain

    \displaystyle \frac{I_h^1(a, t)} {I_h^{1^*}(a)}=\frac{I_v^1(t)} {I_v^{1^*}}. (80)

    It is obvious that the left term \frac{I_h^1(a, t)} {I_h^{1^*}(a)} of (80) is a function with a, t, while the right term \frac{I_v^1(t)} {I_v^{1^*}} is a function with t. So we can assume that I_h^1(a, t)=I_h^{1^*}(a)g(t). Thus we have

    \displaystyle I_v^1=I_v^{1^*}g(t). (81)

    It follows from (2) we can also obtain

    \begin{array}{ll} \displaystyle I_v^{1'}(t) & \displaystyle= S_v\int^\infty_0\beta_v^1(a)I_h^1(a,t)da-(\mu_v+\alpha_v^1)I_v^1,\\[2ex] & \displaystyle= S_v^{1^*}\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)g(t)da-(\mu_v+\alpha_v^1)I_v^1,\\[2ex] & \displaystyle= g(t)S_v^{1^*}\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)da-(\mu_v+\alpha_v^1)I_v^1,\\[2ex] & \displaystyle= g(t)(\mu_v+\alpha_v^1)I_v^{1^*}-(\mu_v+\alpha_v^1)I_v^1,\\[2ex] & \displaystyle=(\mu_v+\alpha_v^1)(I_v^{1^*}g(t)-I_v^1)=0. \end{array} (82)

    Therefore, we can get

    I_v^1=I_v^{1^*}.

    Subsequently, it follows from (80) we have

    I_h^1(a,t)=I_h^{1^*}(a).

    Specially, when a=0, we have I_h^1(0,t)=I_h^{1^*}(0). Thus from (79) we get

    E_h^1(\tau, t)=E_h^{1^*}(\tau).

    Since E_h^i(0,t)=0, then E_h^i(\tau,t)=E_h^i(0,t-\tau)\pi^i_1(\tau)=0 for t>\tau,i=2,\cdots,n. Similarly, we also have I_h^i(0,t)=\int^\infty_0m_h^i(\tau)E_h^i(\tau,t)d\tau=0,\ I_h^i(a,t)=I_h^i(0,t-a)\pi^i_2(a)=0 for t>a. At last we conclude that the largest invariant set in \Theta_2 is the singleton \mathcal{E}_1. Since \Psi(t)\Omega_0^+\subset \Omega^+_0, the global attractor, \mathfrak{T}, is actually contained in \Omega^+_0. Furthermore, the interior global attractor \mathfrak{T} is invariant. By using the above result, we show that the compact global attractor \mathfrak{T}=\{\mathcal{E}_1\}. This completes the proof of Theorem 6.1.


    7. Discussion

    In this paper, we formulate a multi-strain 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, respectively. The formulas for the reproduction number \mathcal R^j_0 of strain j, ~j=1,\cdots,n are obtained from the biological meanings of models. And we define the basic number of the disease as the maximum of the reproduction numbers of each strain. We show that if \mathcal R_0<1, the disease-free equilibrium is locally and globally asymptotically stable. That means the disease dies out and the number of infected with each strain goes to zero. If \mathcal R_0>1, without loss of generality, assuming \mathcal R^1_0=\max\{\mathcal R^1_0,\cdots,\mathcal R^n_0\}>1, we show that the single-strain equilibrium \mathcal E_1 corresponding to strain one exists. The single-strain equilibrium \mathcal E_1 is locally asymptotically stable when \mathcal R^1_0>1 and \mathcal R^i_0<\mathcal R^1_0, ~i=2,\cdots,n.

    The main purpose in this article is to extend the competitive exclusion result established by Bremermann and Thieme in [2], who using a multiple-strain ODE model derives that if multiple strains circulate in the population only the strain with the largest reproduction number persists, the strains with suboptimal reproduction numbers are eliminated. The proof of the competitive exclusion principle is based on the proof of the global stability of the single-strain equilibrium \mathcal E_1. We approach the result by using a Lyapunov function under a stronger condition that

    \displaystyle\frac{\mathcal R^i_0}{\mathcal R^1_0}<\displaystyle\frac{b_i}{b_1}<1, i\neq1. (83)

    Our results do not include the case of

    \max\{\mathcal R^1_0,\cdots,\mathcal R^n_0\}=\mathcal R^1_0=\mathcal R^2_0=\cdots=\mathcal R^m_0>1, m\leq n, m\geq 2.

    According to Proposition 3.3 in [16], where the authors proved and simulated by data that if there is no mutation between two strains and if the basic reproduction numbers corresponding to the two strains are the same, then for the two strain epidemic model there exist many coexistence equilibria, we guess that the coexistence of multi-strains may occur and it is impossible for competitive exclusion in this case.

    From the expression (6) of the basic reproduction number \mathcal R^j_0 corresponding to strain j and the inequality \mathcal R^i_0/\mathcal R^1_0<b_i/b_1,i\neq 1, it follows that

    {r_i}<{r_1},

    where

    r_j=\frac{\beta_h^j}{\mu_v+\alpha^j_v},~\text{for} ~~j=1,2, \cdots, n.

    r_j represents the transmission rate of an infectious vector with strain j during its entire infectious period. The condition (83) implies that the following three inequalities hold at the same time,

    \mathcal R^i_0<\mathcal R^1_0, {r_i}<{r_1}, {b_i}<{b_1}, i\neq1.

    Recall that b_j denotes the probability of a given susceptible vector being transmitted by an infected host with strain j. Then the condition (83) for the occurrence of competition exclusion of strain 1 means that the basic reproduction number corresponding to strain 1, the transmission rate of an infectious vector with strain 1 during its entire infectious period, and the probability of a given susceptible vector being transmitted by an infected host with strain 1 are all biggest comparing to three quantities of other strains.


    Acknowledgments

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


    [1] [ F. Brauer,Z. S. Shuai,P. V. D. Driessche, Dynamics of an age-of-infection cholera model, Math. Biosci. Enger., 10 (2013): 1335-1349.
    [2] [ H. J. Bremermann,H. R. Thieme, A competitive exclusion principle for pathogen virulence, J. Math. Biol, 27 (1989): 179-190.
    [3] [ L. M. Cai,M. Martcheva,X. Z. Li, Competitive exclusion in a vector-host epidemic model with distributed delay, J. Biol. Dynamics, 7 (2013): 47-67.
    [4] [ Y. Dang, Z. Qiu, X. Li and M. Martcheva, Global dynamics of a vector-host epidemic model with age of infection, submitted, 2015.
    [5] [ L. Esteva,A. B. Gumel,C. Vargas-de Leon, Qualitative study of transmission dynamics of drug-resistant malaria, Math. Comput. Modelling, 50 (2009): 611-630.
    [6] [ Z. Feng,J. X. Velasco-Hernandez, Competitive exclusion in a vector host model for the dengue fever, J. Math. Biol, 35 (1997): 523-544.
    [7] [ S. M. Garba,A. B. Gumel, Effect of cross-immunity on the transmission dynamics of two srains of dengue, J. Comput. Math, 87 (2010): 2361-2384.
    [8] [ J. K. Hale, Asymptotic Behavior of Dissipative Systems AMS, Providence, 1988.
    [9] [ W. O. Kermack,A. G. McKendrick, in Contributions to the Mathematical Theory of Epidemics, Contributions to the Mathematical Theory of Epidemics, Royal Society Lond, 115 (1927): 700-721.
    [10] [ A. L. Lloyd, Destabilization of epidemic models with the inclusion of realistic distributions of infectious periods, Royal Society Lond B, 268 (2001): 985-993.
    [11] [ A. L. Lloyd, Realistic distributions of infectious periods in epidemic models: changing patterns of persistence and dynamics, Theor. Popul. Biol., 60 (2001): 59-71.
    [12] [ P. Magal,C. C. McCluskey,G. Webb, Lyapunov functional and global asymptotic stability for an infection-age model, Appl. Anal., 89 (2010): 1109-1140.
    [13] [ M. Martcheva,X. Z. Li, Competitive exclusion in an infection-age structured model with environmental transmission, J. Math. Anal. Appl, 408 (2013): 225-246.
    [14] [ M. Martcheva,H. R. Thieme, Progression age enhanced backward bifurcation in an epidemic model with super-infection, J. Math. Biol, 46 (2003): 385-424.
    [15] [ M. Martcheva,M. Iannelli,X. Z. Li, Subthreshold coexistence of strains: The impact of vaccination and mutation, Math. Biosci. Eng., 4 (2007): 287-317.
    [16] [ D. L. Qian, X. Z. Li and M. Ghosh, Coexistence of the strains induced by mutation Int. J. Biomath. 5 (2012), 1260016, 25 pp.
    [17] [ Z. P. Qiu,Q. K. Kong,X. Z. Li,M. M. Martcheva, The vector host epidemic model with multiple strains in a patchy environment, J. Math. Anal. Appl., 405 (2013): 12-36.
    [18] [ H. R. Thieme,C. Castillo-Chavez, How may infection-age-dependent infectivity affect the dynamics if HIV/AIDs?, SIAM J. Appl. Math., 53 (1993): 1447-1479.
    [19] [ H. R. Thieme, Uniform persistence and permanence for non-autonomous semiflows in population biology, Math. Biosci, 166 (2000): 173-201.
    [20] [ J. X. Yang,Z. P. Qiu,X. Z. Li, Global stability of an age-structured cholera model, Math. Biosci. Enger., 11 (2014): 641-665.
    [21] [ K. Yosida, Functional Analysis Springer-Verlag, Berlin, Heidelberg, New York, 1968.
  • This article has been cited by:

    1. Cuicui Jiang, Wendi Wang, Jiangtao Yang, Threshold conditions for stochastic coexistence of a competition model with Gompertz growth, 2022, 131, 08939659, 108066, 10.1016/j.aml.2022.108066
    2. Xiaoguang Li, Xuan Zou, Liming Cai, Yuming Chen, Global dynamics of a vector-borne disease model with direct transmission and differential susceptibility, 2023, 69, 1598-5865, 381, 10.1007/s12190-022-01745-8
    3. Wendi Wang, Competitive exclusion of two viral strains of COVID-19, 2022, 7, 24680427, 637, 10.1016/j.idm.2022.10.001
    4. Yue Wu, Shenglong Chen, Ge Zhang, Zhiming Li, Dynamic analysis of a stochastic vector-borne model with direct transmission and media coverage, 2024, 9, 2473-6988, 9128, 10.3934/math.2024444
    5. Xiaoguang Li, Liming Cai, Wandi Ding, Modeling the transmission dynamics of a two-strain dengue disease with infection age, 2025, 18, 1793-5245, 10.1142/S1793524524500049
  • Reader Comments
  • © 2017 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(4083) PDF downloads(586) Cited by(5)

Article outline

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog