Processing math: 72%
 

Global stability in a tuberculosis model of imperfect treatment with age-dependent latency and relapse

  • Received: 04 June 2016 Revised: 30 December 2016 Published: 01 October 2017
  • MSC : Primary: 35L60, 92C37; Secondary: 34K20

  • In this paper, an SEIR epidemic model for an imperfect treatment disease with age-dependent latency and relapse is proposed. The model is well-suited to model tuberculosis. The basic reproduction number R0 is calculated. We obtain the global behavior of the model in terms of R0. If R0<1, the disease-free equilibrium is globally asymptotically stable, whereas if R0>1, a Lyapunov functional is used to show that the endemic equilibrium is globally stable amongst solutions for which the disease is present.

    Citation: Shanjing Ren. Global stability in a tuberculosis model of imperfect treatment with age-dependent latency and relapse[J]. Mathematical Biosciences and Engineering, 2017, 14(5&6): 1337-1360. doi: 10.3934/mbe.2017069

    Related Papers:

    [1] P. van den Driessche, Lin Wang, Xingfu Zou . Modeling diseases with latency and relapse. Mathematical Biosciences and Engineering, 2007, 4(2): 205-219. doi: 10.3934/mbe.2007.4.205
    [2] C. Connell McCluskey . Global stability for an SEI epidemiological model with continuous age-structure in the exposed and infectious classes. Mathematical Biosciences and Engineering, 2012, 9(4): 819-841. doi: 10.3934/mbe.2012.9.819
    [3] 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
    [4] 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
    [5] 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
    [6] Xi-Chao Duan, Xue-Zhi Li, Maia Martcheva . Dynamics of an age-structured heroin transmission model with vaccination and treatment. Mathematical Biosciences and Engineering, 2019, 16(1): 397-420. doi: 10.3934/mbe.2019019
    [7] Gang Huang, Edoardo Beretta, Yasuhiro Takeuchi . Global stability for epidemic model with constant latency and infectious periods. Mathematical Biosciences and Engineering, 2012, 9(2): 297-312. doi: 10.3934/mbe.2012.9.297
    [8] 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
    [9] Georgi Kapitanov . A double age-structured model of the co-infection of tuberculosis and HIV. Mathematical Biosciences and Engineering, 2015, 12(1): 23-40. doi: 10.3934/mbe.2015.12.23
    [10] Yuming Chen, Junyuan Yang, Fengqin Zhang . The global stability of an SIRS model with infection age. Mathematical Biosciences and Engineering, 2014, 11(3): 449-469. doi: 10.3934/mbe.2014.11.449
  • In this paper, an SEIR epidemic model for an imperfect treatment disease with age-dependent latency and relapse is proposed. The model is well-suited to model tuberculosis. The basic reproduction number R0 is calculated. We obtain the global behavior of the model in terms of R0. If R0<1, the disease-free equilibrium is globally asymptotically stable, whereas if R0>1, a Lyapunov functional is used to show that the endemic equilibrium is globally stable amongst solutions for which the disease is present.


    1. Introduction

    Mathematical modeling is a very important tool in analyzing the propagation and controlling of infectious diseases. Age structure is an important characteristic in the modeling of some infectious diseases. The first formulation of a partial differential equation(PDE) for the age distribution of a population was due to McKendrick [21]. Since the seminar papers by Kermack and McKendrick [13]-[15], age structure models have been used extensively to study the transmission dynamics of infectious diseases, we refer to the monographs by Hoppensteadt [11], Iannalli [12] and Webb [30] on this topic.

    As an ancient disease, TB peaked and declined by 1940's before it became curable, while the downtrend stopped in the middle 1980's and 1990's. As one of the top 3 deadly infectious diseases, TB would cause a higher death rate if not treated, while the disease would be latent in an individual body for months, years or even decades before it outbreaks. McCluskey [20] pointed out that the risk of activation can be modeled as a function of duration age, and this form can be used to describe more general latent period via introducing the duration age in the latent class as a variable.

    On the other hand, for the infectious tuberculosis, the removed individuals often have a higher relapse rate. Actually, the recurrence as an important feature of some animal and human diseases has been studied extensively, see [4], [23]. For instance, van den Driessche and coauthors in [4], [5] established two models with a constant relapse period and a general relapse distribution respectively, which showed the threshold property of the basic reproduction number. It is interest to investigate the model with age-dependent relapse rate and to determine whether the threshold property can be preserved or not.

    Recently, Wang et al. [24]-[27] considered the global stability of nonlinear age-structured models, Liu et al. [17] introduced age-dependent latency and relapse into an SEIR epidemic model and the local stability and global stability of equilibria are obtained by analyzing the corresponding characteristic equations and constructing the proper Volterra-type Lyapunov functionals, respectively. Wang et al. [28] proposed an SVEIR epidemic model with media impact, age-dependent vaccination and latency, and discussed the global dynamics of the age-structure model.

    However, most of the models assumed that TB would show neither its clinical symptoms nor its infectivity during its latent period, while in fact, TB has many early clinical symptoms such as fever, fatigability, night sweat, chest pain, hemoptysis and so on. Here we formulate and analyze an SEIR epidemic model with continuous age dependent latency and relapse. We assumed, as the development of the disease, TB is infectious during its latent period with less infectivity and incomplete treatment comparing with outbreak period. Although epidemic models with age-dependent have been studied extensively, there have been still inadequate results on the full global stability. In this paper, we employ the method developed by Webb [30] for age-dependent models, namely integrating solutions along the characteristics to obtain an equivalent integral equation. We obtain the basic reproduction number in virtue of the method in [7]. Moreover, we study the asymptotic smoothness of the semi-flow generated by the system and the existence of a global attractor [3], [19]. Finally, we show the global stability of equilibria via constructing the proper Volterra-type Lyapunov functionals. For more details concerning the current Lyapunov functionals approach, we refer the reader to recent work [2].

    This paper will be organized as follows: In Section 2, we formulate our general SEIR tuberculosis model with latent age and relapse age which is described by a coupled system of ODEs and PDEs. In Section 3, we investigate the existence of equilibria and obtain the expression of the basic reproduction number R0. In Section 4, the local asymptotic stability of the equilibria will be derived. In Section 5, we present the results about uniform persistence. In section 6, we deal with the global stability of equilibria. Finally, some numerical simulations and useful discussions are made in the last section.


    2. Model formulation

    The total population is decomposed into four disjoint subclasses, susceptible class S, latent class E, infectious class I, and removed class R. More precisely, let S(t) denote the number of susceptible individuals at time t. Susceptible individuals would become new infected ones after they contact with infectious individuals at a rate β, while they enter a stage when they are infected with the disease but have little infectivity. This stage is often called latent stage, which maybe enter into the stage of removed class R by receiving treatment at a rate μ. The density of individuals in the latent class is denoted by e(t,a) where t is the duration time spenting in this class and a is called the latent-stage progression age, denoting E(t)=+0e(t,a)da the total density of latent individuals. The number of individuals in the class I at time t is I(t). The removal rate from latent class E to infectious class I is given by the function σ(a). Thus, the total rate at which individuals progress into the infectious class alive is +0σ(a)e(t,a)da. Infectious individuals come into the removed class after recovery due to complete treatment. Let r1 be the recovery rate from the infectious class. The density of individuals in removed class is denoted by r(t,c), where c represents the relapse age, denoting R(t)=0r(t,c)dc the total density of removed individuals. In fact, infectious individuals might come into the latent class E due to incomplete treatment at the rate r2. Due to the relapse of the disease, the age-dependent relapse rate in the removed class is given by the function k(c). The total rate at which individuals relapse into the infectious class alive is given by the quantity +0k(c)r(t,c)dc. We also denote Λ, δe,δi as the density of the recruitment into the susceptible class, the additional death rates induced by the latent disease and infectious disease. The parameter b is the natural death rate of all individuals. All recruitment of the population enters the susceptible class and occurs with constant flux Λ. Further, all parameters are assumed to be positive. FIGURE 1 shows the schematic flow diagram of our model which can be described by a system of ordinary and partial differential equations

    Figure 1. Here is the Model of TB.
    {dS(t)dt=ΛbS(t)βS(t)I(t),(t+a)e(t,a)=(b+δe+μ(a)+σ(a))e(t,a),dI(t)dt=(r1+r2+b+δi)I(t)+0σ(a)e(t,a)da+0k(c)r(t,c)dc,(t+c)r(t,c)=(k(c)+b)r(t,c), (1)

    with boundary conditions

    {e(t,0)=βS(t)I(t)+r2I(t),r(t,0)=r1I(t)+0μ(a)e(t,a)da, (2)

    and initial conditions

    S(0)=S0,e(0,a)=e0(a),I(0)=I0,r(0,c)=r0(c), (3)

    where S0,I0R+, and e0(a),r0(c)L1+(0,) which is the nonnegative and Lebesgue integrable space of functions on [0,+).

    In order to simplify the later derivation, we make the following hypotheses about the parameters of the system (1)

    (H1) σ(a),k(c),μ(a)L1+(0,+), with respective essential upper bounds ˉσ,ˉk,ˉμ;

    (H2) σ(a),k(c),μ(a)L1+(0,+) are Lipschitz continuous on R+, with Lipschitzcoefficients Mσ,Mk and Mμ respectively;

    (H3) There exists b0[0,b), such that σ(a),k(c),μ(a)b0,foralla,c>0.

    For a,c>0, we denote

    ε(a)=σ(a)+μ(a)+b+δe,η(c)=k(c)+b,ρ1(a)=ea0ε(s)ds,ρ2(c)=ec0η(s)ds,θ1=+0σ(a)ρ1(a)da,θ2=+0k(c)ρ2(c)dc,θ3=+0μ(a)ρ1(a)da.

    According to Webb [30], by solving the PDE parts of (1) along the characteristic lines ta=const and tc=const respectively, we obtain

    e(t,a)={e(ta,0)ea0ε(s)ds,t>a0,e0(at)eaatε(s)ds,at0, (4)
    r(t,c)={r(tc,0)ec0η(s)ds,t>c0,r0(ct)ecctη(s)ds,ct0. (5)

    Define the space of functions X as

    X:=R+×L1+(0,+)×R+×L1+(0,+)

    equipped with the norm

    (x1,x2,x3,x4)X=|x1|++0|x2(a)|da+|x3|++0|x4(c)|dc.

    The norm has the biological interpretation of giving the total population size. The initial conditions (3) for system (1) can be rewritten as x0=(S0,e(0,),I0,r(0,))X. Using standard methods we can verify the existence and uniqueness of solutions to model (1) with the boundary conditions (2) and initial conditions (3) (see Iannelli [12] and Webb [30]). Meanwhile, we can claim that any solution of system (1) with nonnegative initial conditions remains nonnegative. The nonnegativity of e(t,a) and r(t,c) follows from (4) and (5). Next, we shall show that S(t)>0 for t0 and I(t)>0 for t0. Otherwise, assume that S(t) would lose its positivity for the first time at t1>0, i.e., S(t1)=0. However, from the first equation of (1) we obtain

    S(t1)=ebt1t10βI(τ)dτ{S(0)+t10ebs+s0βI(τ)dτΛds}>0.

    Similarly, assume that I(t) would lose its positivity for the first time at t2>0, i.e., I(t2)=0. However, from the third equation of (1) we obtain

    I(t2)=e(r1+r2+b+δi)t2t20(0σ(a)e(t,a)da+0k(c)r(t,c)dc)e(r1+r2+b+δi)s)ds+e(r1+r2+b+δi)t2I(0)>0.

    Thus S(t)>0 and I(t)>0 are true for t0. This verifies our claim.

    Let us consider a function N(t)=S(t)++0e(t,a)da+I(t)++0r(t,c)dc, which is the total population at time t. We can easily see that the time derivative of N along solutions of model (1) is

    ddtN(t)=ddtS(t)+ddt+0e(t,a)da+ddtI(t)+ddt+0r(t,c)dc.

    due to ρ1(0)=1,dρ1(a)da=ε(a)ρ1(a), we have

    ddt+0e(t,a)da=ddt(t0e(ta,0)ρ1(a)da++te0(at)ρ1(a)ρ1(at)da)=ddtt0(βS(ta)I(ta)+r2I(ta))ρ1(a)da+ddt+te0(at)ρ1(a)ρ1(at)da=ddtt0(βS(τ)I(τ)+r2I(τ))ρ1(tτ)dτ+ddt+0e0(τ)ρ1(t+τ)ρ1(τ)dτ=βS(t)I(t)+r2I(t)+0ε(a)e(t,a)da.

    Similarly, by using ρ2(0)=1,dρ2(c)dc=η(c)ρ2(c), we can get

    ddt+0r(t,c)dc=ddtt0(r1I(tc)++0μ(a)e(tc,a)da)ρ2(c)dc+ddt+tr0(ct)ρ2(c)ρ2(ct)dc=r1I(t)++0μ(a)e(t,a)da+0η(c)r(t,c)dc.

    Hence, we have

    ddt(S(t)++0e(t,a)da+I(t)++0r(t,c)dc)=ΛbS(t)βS(t)I(t)+βS(t)I(t)+r2I(t)+0ε(a)e(t,a)da(r1+r2+b+δi)I(t)+0σ(a)e(t,a)da+0k(c)r(t,c)dc+r1I(t)++0μ(a)e(t,a)da+0η(c)r(t,c)dc=ΛbS(t)+0br(t,c)dcbI(t)+0(b+δe)e(t,a)daδiI(t)ΛbN(t).

    It follows from the variation of constants formula that N(t)Λb, for any t0, which implies that

    Ω={(S(t),e(t,),I(t),r(t,))R+×L1+(0,+)×R+×L1+(0,+):N(t)Λb}

    is positively invariant absorbing set for system (1).

    Proposition 1. If x0X and x0XM for some constant MΛb, then the following statements hold for t0,

    ()0S(t),+0e(t,a)da,I(t),+0r(t,c)dcM,()e(t,0)βM2+r2M,r(t,0)r1M+ˉμM.

    For convenience, we rewrite (1) as follows

    {(t+a)e(t,a)=(b+δe+μ(a)+σ(a))e(t,a),(t+c)r(t,c)=(k(c)+b)r(t,c),dV(t)dt=G(e(t,a),r(t,c),V(t))CV(t),e(t,0)=βS(t)I(t)+r2I(t),e(0,a)=e0(a),r(t,0)=r1I(t)+0μ(a)e(t,a)da,r(0,c)=r0(c),V(0)=V0, (6)

    where

    V(t)=(S(t)I(t)),
    C=(b00r1+r2+b+δi),
    G(e(t,a),r(t,c),V(t))=(ΛβS(t)I(t)0σ(a)e(t,a)da+0k(c)r(t,c)dc,).

    Set Z=Y×R2, where Y=R×L1(R+,R), for any (αϕ)Y, we have (αϕ)=|α|+ϕL1(R+,R). Furthermore, we define

    Z+=Y+×R2+,Z0=Y0×R2,Z0+=Z0Z+,

    where

    Y+=R+×L+(R+,R),Y0={0}×L(R+,R).

    We define A1:Dom(A1)YY, by

    A1(0ϕ1)=(ϕ1(0)ϕ1(b+δe+μ(a)+σ(a))ϕ1)

    with Dom(A1)={0}×W1,1(R+,R). If λ is a complex number with Reλ>(b+δe), then λρ(A1) which is the resolvent set of A1. Moreover, if λρ(A1) and

    (λIA1)1(θ1ψ1)=(0ϕ1),

    then we can get

    ϕ1(a)=e(λ+b+δe)aθ1+a0eas(μ(l)+σ(l))dle(λ+b+δe)(as)ψ1(s)ds.

    Similarly, we define A2:Dom(A2)YY, by

    A2(0ϕ2)=(ϕ2(0)ϕ2(b+k(c))ϕ2)

    with Dom(A2)={0}×W1,1(R+,R), we can obtain

    ϕ1(a)=e(λ+b)cθ2+c0ecs(k(l))dle(λ+b)(cs)ψ2(s)ds.

    Thus (6) can be rewritten as

    {ddt(0e(t,))=A1(0e(t,))+F1((0e(t,)),(0r(t,)),V(t)),ddt(0r(t,))=A2(0r(t,))+F2((0e(t,)),(0r(t,)),V(t)),dV(t)dt=CV(t)+F3((0e(t,)),(0r(t,)),V(t)),e(0,a)=e0(a),r(0,c)=r0(c),V(0)=V0, (7)

    where

    F1((0e(t,)),(0r(t,)),V(t))=(βS(t)I(t)+r2I(t)0),F2((0e(t,)),(0r(t,)),V(t))=(r1I(t)+0μ(a)e(t,a)da0),F3((0e(t,)),(0r(t,)),V(t))=(ΛβS(t)I(t)0σ(a)e(t,a)da+0k(c)r(t,c)dc).

    Let L:D(L)XX be the linear operator defined by

    L(u(t))=(A1(0e(t,)),A2(0r(t,)),CV(t)),

    where D(L)=P×R2 with P={0}×W1,1(R+,R)×{0}×W1,1(R+,R). It follows that X0=¯D(L) and X0+=¯D(L)X+. So ¯D(L)=X0 is not dense in X. We consider the nonlinear map F:¯D(L)X defined by

    F(u(t))=(F1(u(t))F2(u(t))F3(u(t))).

    Therefore (7) can be rewritten as an abstract Cauchy problem

    {du(t)dt=L(u(t))+F(u(t)),u(0)=((0e0()),(0r0()),V0). (8)

    By using the results in Magal [19] and Chen et al. [3], there exists a uniquely determined semiflow {U(t)}t0 on X0+ such that, for each

    u=((0e(t,)),(0r(t,)),V(0)), there exists a unique continuous map UC(R+,X0+), which is an integrated solution of the Cauchy problem (8), that is, for all t0, t0U(s)udsD(L) and U(t)u=u+Lt0U(s)uds+t0F(U(s)u)ds. And Ω is the positively invariant absorbing set under the semi-flow U can be verified, that is, U(t)ΩΩ and for each xX0+, d(U(t),Ω):=infyΩU(t)xy∥→0, as t which means that the semi-flow {U(t)}t0 is bound dissipative on X0+. A semi-flow U(t,x0):R+×XX is called asymptotically smooth if each forward invariant bounded closed set is attracted by a nonempty compact set [8], [22]. In order to obtain global properties of the dynamics of the semi-flow U(t), it is important to prove the asymptotically smooth of semi-flow U(t). First we give the following useful lemma.

    Lemma 2.1. ([1]) Let DR. For j=1,2, suppose fj:DR is a bounded Lipschitz continuous function with bound Kj and Lipschitz coefficient Mj. Then the product function f1f2 is Lipschitz with coefficient K1M2+K2M1.

    Lemma 2.2. ([1]) If the following two conditions hold, then the semi-flow U(t,x0)=ϕ(t,x0)+φ(t,x0):R+×XX is asymptotically smooth in X:

    (ⅰ) There exists a continuous function v:R+×R+R+ such that limt+v(t,h)=0 and ϕ(t,x0)Xv(t,h) if x0Xh;

    (ⅱ) For t0,φ(t,x0) is completely continuous.

    In order to prove Lemma 2.2, we first decompose U:R+×XX into the following two operators ϕ(t,x0),φ(t,x0):R+×XX, ϕ(t,x0)=(0,y2(t,),0,y4(t,)), φ(t,x0)=(S(t),˜y2(t,),I(t),˜y4(t,)), where

    y2(t,a)={0,t>a0,e0(at)ρ1(a)ρ1(at),at0,y4(t,c)={0,t>c0,r0(ct)ρ2(c)ρ2(ct),ct0. (9)
    ˜y2(t,a)={e(ta,0)ρ1(a),t>a0,0,at0,˜y4(t,c)={r(tc,0)ρ2(c),t>c0,0,ct0. (10)

    In order to verify condition (ⅰ) of Lemma 2.2, we need to prove the following proposition.

    Proposition 2. For h>0, let v(t,h)=he(b+2b0+δe)t. Then v(t,h)0 as t+ and ϕ(t,x0)Xv(t,h) if x0Xh.

    Proof. It is obvious that v(t,h)0 as t+, with the help of (9) and (H3), we have

    ϕ(t,x0)X=|0|++0|y2(t,a)|da+|0|++0|y4(t,c)|dc=+t|e0(at)ρ1(a)ρ1(at)|da++t|r0(ct)ρ2(c)ρ2(ct)|dc=+0|e0(τ)eτ+tτε(s)ds|dτ++0|r0(τ)eτ+tτη(s)ds)|dτe(b+2b0+δe)t(|0|++0|e0(τ)|dτ+|0|++0|r0(τ)|dτ)=e(b+2b0+δe)tx0X,

    by the known condition x0Xh,x0Ω, we have ϕ(t,x0)Xhe(b+2b0+δe)tv(t,h).

    Lemma 2.3. ([1]) Let KLp(0,+) be closed and bounded where p1. Then K is compact if the following conditions hold

    () limh0+0|f(z+h)f(z)|pdz=0uniformlyforfK,

    () limh++h|f(z)|pdz=0uniformlyforfK.

    Proposition 3. For t0,ϕ(t,x0) is completely continuous.

    Proof. According to Proposition 1(ⅰ), S(t) and I(t) remain in the compact set [0,Λ/b][0,M], where MΛ/b. Thus, it needs only to show that ˜y2(t,a) and ˜y4(t,c) remain in a precompact subset of L1+(0,+), which is independent of x0Ω. It suffices to verify that (ⅰ) and (ⅱ) in Lemma 2.3 hold. Now, from Proposition 1(ⅱ) and (10) we have

    ˜y2(t,a)(βM2+r2M)e(b+2b0+δe)a. (11)
    +0|˜y2(t,a+h)˜y2(t,a)|da=th0|e(tah,0)ρ1(a+h)e(ta,0)ρ1(a)|da+tth|0e(ta,0)ρ1(a)|dath0e(tah,0)|ρ1(a+h)ρ1(a)|da+th0ρ1(a)|e(tah,0)e(ta,0)|da+tth|e(ta,0)ρ1(a)|da. (12)

    Recall that ρ1(a)=ea0(σ(s)+μ(s)+b+δe)dse(2b0+b+δe)a1, then ρ1(a) is a non-increasing function with respect to a, we have

    th0|ρ1(a+h)ρ1(a)|da=th0ρ1(a)dathρ1(a)da=th0ρ1(a)dathhρ1(a)datthρ1(a)da=h0ρ1(a)datthρ1(a)dah. (13)

    From Proposition 1 and (H1), we find that |dS(t)/dt| is bounded by MS=Λ+bM+βM2 and |dI(t)/dt| is bounded by MI=ˉσM+ˉkM+(b+δi+r1+r2)M. Therefore, S() and I() are Lipschitz on [0,+) with coefficients MS and MI. By Lemma 2.1, S()I() is Lipschitz on [0,+) with coefficient MSI=M(MS+MI). Thus

    th0ρ1(a)|e(tah,0)e(ta,0)|da
    th0ρ1(a)(|βS(tah)I(tah)βS(ta)I(ta)|+|r2I(tah)r2I(ta)|)dath0(βMSI+r2MI)(h)e(b+2b0+δe)ada=(βMSI+r2MI)hb+2b0+δe(1e(b+2b0+δe)(th))(βMSI+r2MI)hb+2b0+δe. (14)

    From (12)-(14), we have

    +0|˜y2(t,a+h)˜y2(t,a)|da(βMSI+r2MIb+2b0+δe+2(βM2+r2M))h

    which converges uniformly to 0 as h0. The condition (i) in Lemma 2.3 is proved for ˜y2(t,a).

    From (11) we have

    limh++h|˜y2(t,a)|dalimh++h(βM2+r2M))e(b+2b0+δe)ada=limh+βM2+r2Mb+2b0+δee(b+2b0+δe)h=0

    which meet the condition (ii) in Lemma 2.3. Similarly, y4(t,c) also satisfies the conditions of Lemma 2.3.

    Theorem 2.4. The semi-flow {U(t)}t0 generated by system (1) is asymptotically smooth, and has a global attractor A contained in X, which attracts the bounded sets of X.


    3. Existence of the equilibria

    Now we consider the existence of equilibria of system (1). The steady state (S,e(),I,r()) of system (1) satisfies the equalities

    {0=ΛbSβSI,ddae(a)=ε(a)e(a),0=(r1+r2+b+δi)I+0σ(a)e(a)da+0k(c)r(c)dc,ddcr(c)=η(c)r(c), (15)

    with boundary conditions

    {e(0)=βSI+r2I,r(0)=r1I+0μ(a)e(a)da. (16)

    From the second equation of (15) and the first equation of (16), we obtain

    e(a)=(βSI+r2I)ea0ε(s)ds. (17)

    Similarly, by using the fourth equation of (15) and the second equation of (16), we get

    r(c)=(r1I++0μ(a)e(a)da)ec0η(s)ds. (18)

    If I=0, then we have e(a)=0 and r(c)=0 respectively from (17) and (18). Furthermore, it is easy to know that S0=Λb from the first equation of (15). Thus, system (15) has a disease-free equilibrium E0, and

    E0=(S0,0,0,0),whereS0=Λb.

    In order to find any endemic equilibrium, we introduce the basic reproduction number R0, which is the average number of new infections generated by a single newly infectious individual during the full infectious period [7]. It is given by the following expression

    R0=βS0(θ1+θ2θ3)r1+r2+b+δir1θ2r2(θ1+θ2θ3).

    Now, if I0, substituting (17) and (18) into the third equation of (15), we have

    0=(r1+r2+b+δi)I++0σ(a)(βSI+r2I)ea0ε(s)dsda++0k(c)(r1I++0μ(a)(βSI+r2I)ea0ε(s)dsda)ec0η(s)ds=(r1+r2+b+δi)I+(βSI+r2I)(θ1+θ2θ3)+r1Iθ2. (19)

    Thus, combining the first equation of (15) and the equation (19), we get

    S=ΛbR0,I=bβ(R01). (20)

    If R0>1, we get e(a)>0 and r(c)>0 from (17) and (18). Therefore, system (15) has a unique positive endemic equilibrium E, and

    E=(S,e(),I,r()),

    where

    S=ΛbR0,e(a)=(βSI+r2I)ρ1(a),I=bβ(R01),r(c)=((βSI+r2I)θ3+r1I)ρ2(c).

    4. Local asymptotic stability of the equilibria

    In this section, sufficient conditions for the local asymptotic stability of the equilibria will be derived.

    Theorem 4.1. The disease-free equilibrium E0 is locally asymptotically stable in the positive variant set Ω if R01 and unstable if R0>1.

    Proof. First, we introduce the change of variables as follows

    x1(t)=S(t)S0,x2(t,a)=e(t,a),x3(t)=I(t),x4(t,c)=r(t,c).

    Linearizing the system (1) about disease-free equilibrium E0, we obtain the following system

    {dx1(t)dt=bx1(t)βΛbx3(t),(t+a)x2(t,a)=(b+δe+μ(a)+σ(a))x2(t,a),dx3(t)dt=(r1+r2+b+δi)x3(t)++0σ(a)x2(t,a)da++0k(c)x4(t,c)dc,(t+c)x4(t,c)=(k(c)+b)x4(t,c),x2(t,0)=(βΛb+r2)x3(t),x4(t,0)=r1x3(t)++0μ(a)x2(t,a)da. (21)

    Set

    x1(t)=x01eλt,x2(t,a)=x02(a)eλt,x3(t)=x03eλt,x4(t,c)=x04(c)eλt, (22)

    where x01,x02(a),x03,x04(c) will be determined. Substituting (22) into (21), we get

    λx01=bx01βΛbx03, (23)
    {λx02(a)+dx02(a)da=(b+δe+μ(a)+σ(a))x02(a),x02(0)=(βΛb+r2)x03, (24)
    λx03=(r1+r2+b+δi)x03++0σ(a)x02(a)da++0k(c)x04(c)dc, (25)
    {λx04(c)+dx04(c)dc=(k(c)+b)x04(c),x04(0)=r1x03++0μ(a)x02(a)da. (26)

    Integrating the first equation of (24) from 0 to a yields

    x02(a)=(βΛb+r2)x03e(λ+b+δe)aa0(σ(s)+μ(s))ds. (27)

    Similarly, we have from (26) that

    x04(c)=(r1x03++0μ(a)x02(a)da)e(λ+b)cc0k(s)ds. (28)

    Substituting (27) and (28) into (25) and solving (25) gives

    λ=(r1+r2+b+δi)++0k(c)r1e(λ+b)cc0k(s)dsdc++0σ(a)(βΛb+r2)e(λ+b+δe)aa0(σ(s)+μ(s)dsda++0k(c)+0μ(a)(βΛb+r2)e(λ+b+δe)aa0(σ(s)+μ(s))dsdae(λ+b)cc0k(s)dsdc (29)

    which is the characteristic equation. Let

    F(λ)=+0k(c)+0μ(a)(βΛb+r2)e(λ+b+δe)aa0(σ(s)+μ(s))dsdae(λ+b)cc0k(s)dsdcλ(r1+r2+b+δi)++0σ(a)(βΛb+r2)e(λ+b+δe)aa0(σ(s)+μ(s)dsda++0k(c)r1e(λ+b)cc0k(s)dsdc.

    Obviously, F(λ) is a continuously differential function and satisfies

    F(λ)=(βΛb+r2)+0aσ(a)e(λ+b+δe)aa0(σ(s)+μ(s))dsdaa+0k(c)e(λ+b)cc0k(s)dsdc+0μ(a)(βΛb+r2)e(λ+b+δe)aa0(σ(s)+μ(s))dsdac+0k(c)e(λ+b)cc0k(s)dsdc+0μ(a)(βΛb+r2)e(λ+b+δe)aa0(σ(s)+μ(s))dsdar1+0ck(c)e(λ+b)cc0k(s)dsdc1<0,

    and

    limλ+F(λ)=,limλF(λ)=+.

    Thus, we know (29) has a unique real root λ. Obviously,

    F(0)=[(r1+r2+b+δi)r1θ2r2(θ1+θ2θ3)](R01),

    we have λ<0,ifR0<1, and λ>0,ifR0>1. Let λ=x+yi be an arbitrary complex root to (29), then

    0=F(λ)=F(x+yi)F(x)

    which means that λ>x. Thus, all the roots of (29) have negative real part if and only if R01 and have at least one eigenvalue with positive real part if R0>1. Therefore we have that the disease-free equilibrium E0 is locally asymptotically stable if R01 and unstable if R0>1.

    Theorem 4.2. The unique endemic equilibrium E is locally asymptotically stable if R0>1.

    Proof. First, we introduce the perturbation variables as follows

    y1(t)=S(t)S,y2(t,a)=e(t,a)e(a),y3(t)=I(t)I,y4(t,c)=r(t,c)r(c).

    Linearizing system (1) at the endemic equilibrium E yields the following system

    {dy1(t)dt=by1(t)βIy1(t)βSy3(t),(t+a)y2(t,a)=(b+δe+μ(a)+σ(a))y2(t,a),dy3(t)dt=(r1+r2+b+δi)y3(t)++0σ(a)y2(t,a)da++0k(c)y4(t,c)dc,(t+c)y4(t,c)=(k(c)+b)y4(t,c),y2(t,0)=βy1(t)I+βSy3(t)+r2y3(t),y4(t,0)=r1y3(t)++0μ(a)y2(t,a)da. (30)

    Set

    y1(t)=y01eλt,y2(t,a)=y02(a)eλt,y3(t)=y03eλt,y4(t,c)=y04(c)eλt, (31)

    Substituting (31) into (30) gives

    λy01=by01βIy01βSy03, (32)
    {dy02(a)da=(λ+b+δe+μ(a)+σ(a))y02(a),y02(0)=βIy01+βSy03+r2y03, (33)
    λy03=(r1+r2+b+δi)y03++0σ(a)y02(a)da++0k(c)y04(c)dc, (34)
    {dy04(c)dc=(λ+k(c)+b)y04(c),y04(0)=r1y03++0μ(a)y02(a)da. (35)

    Integrating the first equation of (33) and (35) from 0 to a and from 0 to c respectively, together with the boundary conditions, yields

    y02(a)=(βIy01+βSy03+r2y03)e(λ+b+δe)aa0(σ(s)+μ(a))ds,y04(c)=(r1y03++0μ(a)y02(a)da)e(λ+b)cc0k(s)ds. (36)

    substituting the above two equations into (34) and solving (34) we get

    λy03=(βIy01+βSy03+r2y03)(K1(λ)+K2(λ)K3(λ))+r1y03K2(λ)(r1+r2+b+δi)y03, (37)

    where

    K1(λ)=+0σ(a)e(λ+b+δe)aa0(σ(s)+μ(s))dsda,K2(λ)=+0k(c)e(λ+b)cc0k(s)dsdc,K3(λ)=+0μ(a)e(λ+b+δe)aa0(σ(s)+μ(a))dsda,

    By combining (37) and (32) we obtain the characteristic equation

    det(λ+b+βIβSβI(K1(λ)+K2(λ)K3(λ))M)=0.

    where M=(βS+r2)(K1(λ)+K2(λ)K3(λ))+r1K2(λ)(λ+r1+r2+b+δi), that is

    M=β2SIλ+b+βI(K1(λ)+K2(λ)K3(λ)). (38)

    It follows from (20) that (38) can also be rewritten as

    (βS0R0+r2)(K1(λ)+K2(λ)K3(λ))+r1K2(λ)=βbS0(R01)(λ+bR0)R0(K1(λ)+K2(λ)K3(λ))+λ+r1+r2+b+δi. (39)

    Note that Ki(λ)<0,i=1,2,3. Thus, Ki(λ),i=1,2,3 is decreasing. Further, Ki(0)=θi,i=1,2,3. Assume that Reλ0, then K1(λ)θ1, K2(λ)θ2 and K3(λ)θ3 hold. Hence, the modulus of the left-hand side of (39) satisfies

    (βS0R0+r2)(K1(λ)+K2(λ)K3(λ))+r1K2(λ)(βS0R0+r2)(θ1+θ2θ3)+r1θ2=r1+r2+b+δi

    which, together with (39), leads to

    βbS0(R01)(λ+bR0)R0(K1(λ)+K2(λ)K3(λ))+λ+r1+r2+b+δi∣≤r1+r2+b+δi. (40)

    Since R0>1, hence

    βbS0(R01)(λ+bR0)R0(K1(λ)+K2(λ)K3(λ))+λ0. (41)

    that is Reλ0. There is a contradiction. This means that all roots of (39) have negative real parts. Consequently, the endemic equilibrium E of (1) is locally asymptotically stable if R0>1.


    5. Uniform persistence

    In this section, we study the uniform persistence of system (1). Define M0={(S,I,0,0,e,r)TX0+:I+0e(a)da+0r(c)dc>0}, and M0=X0+M0.

    Theorem 5.1. M0 and M0 are both positively invariant under the semiform {U(t)}t0 generated by system (1) on X0+. Moreover, the infection-free equilibrium E0=(S0,0,0,0,0L1,0L1) is globally asymptotically stable for the semiflow {U(t)}t0 restricted to M0.

    Proof. Let (S0,I0,0,0,e0,r0)M0,T(t)=I(t)+0e(t,a)da+0r(t,c)dc. It follows that

    T(t)max{(r1+r2+b+δi),(b+δe+μmax)}T(t),

    where μmax=esssupa(0,)μ(a). Then,

    T(t)emax{(r1+r2+b+δi),(b+δe+μmax)}tT(0).

    This completes the fact that U(t)M0M0. Now let (S0,I0,0,0,e0,r0)M0, using (4) and (5), we easily find that I(t)=0, for t0, and

    0e(t,a)da=t0e(ta,0)ea0ε(s)dsda+te0(at)eaatε(s)dsda=t0[βS(ta)I(ta)+r2I(ta)]ea0ε(s)dsda+te0(at)eaatε(s)dsdaeεminte0L10,t.

    Similarly,

    0r(t,c)dc=t0r(tc,0)ec0η(s)dsdc+tr0(ct)ecctη(s)dsdceηmintr0L10,t.

    Thus U(t)M0M0. Let (S0,I0,0,0,e0,r0)M0, we obtain

    {dIdt=(r1+r2+b+δi)I(t)+0σ(a)e(t,a)da+0k(c)r(t,c)dc,(t+a)e(t,a)=(b+δe+μ(a)+σ(a))e(t,a),e(t,0)=βSI+r2I,(t+c)r(t,c)=(k(c)+b)r(t,c),r(t,0)=r1I(t)+0μe(t,a)da,I(0)=0,e(0,a)=e0(a),r(0,c)=r0(c).

    Since S(t)S0 as t is large enough, we get I(t)˜I(t),e(t,a)˜e(t,a) and r(t,c)˜r(t,c), where

    {d˜Idt=(r1+r2+b+δi)˜I+0σ(a)˜e(t,a)da+0k(c)˜r(t,c)dc,(t+a)˜e(t,a)=(b+δe+μ(a)+σ(a))˜e(t,a),˜e(t,0)=βS0˜I+r2˜I,(t+c)˜r(t,c)=(k(c)+b)˜r(t,c),˜r(t,0)=r1˜I+0μ(a)˜e(t,a)da,˜I(0)=0,˜e(0,a)=e0(a),˜r(0,c)=r0(c). (42)

    By the formulations (4), (5), we have

    ˜e(t,a)={˜e(ta,0)ea0ε(s)ds,t>a0,e0(at)eaatε(s)ds,at0. (43)
    ˜r(t,c)={˜r(tc,0)ec0η(s)ds,t>c0,r0(ct)ecctη(s)ds,ct0. (44)

    Substituting (43) and (44) into the first equation of (42), with the help of the third and the fifth equations of (42), we obtain

    {d˜I(t)dt=(H1+H2+H3+H4)˜I(t)+Fe(t)+Fr(t)+Fer(t),˜I(0)=0, (45)

    where

    H1=(r1+r2+b+δi),H2=t0σ(a)(βS0+r2)ea0ε(s)dsda,H3=t0k(c)r1ec0η(s)dsdc,H4=t0k(c)t0μ(a)(βS0+r2)ea0ε(s)dsdaec0η(s)dsdc,Fe(t)=tσ(a)e0(at)eaatε(s)dsda,Fr(t)=tk(c)r0(ct)ecctη(s)dsdc,Fer(t)=t0k(c)tμ(a)e0(at)eaatε(s)dsdaec0η(s)dsdc.

    It's simple to know that, for each t,

    Fe(t)eεminttσ(a)e0(at)da=0,Fr(t)eηminttk(c)r0(ct)dc=0,Fer(t)t0k(c)eεminttμ(a)e0(at)daec0η(s)dsdc=0.

    Then, we know that equation (45) has a unique solution ˜I(t)=0 and we obtain ˜e(t,0)=0, ˜r(t,0)=0 from the third and fifth equations of (42). If 0a<t, according to (43), we have ˜e(t,a)=0. Similarly, If 0\leq c < t , according to (44), we have \tilde{r}(t, c)=0. If a\geq t, \parallel\tilde{e}(t, a)\parallel_{L^{1}}\leq e^{-\varepsilon_{min}t}\parallel e_{0}\parallel_{L^{1}}, if c\geq t, \parallel\tilde{r}(t, c)\parallel_{L^{1}}\leq e^{-\eta_{min}t}\parallel r_{0}\parallel_{L^{1}}, which yields that \tilde{e}(t, a)\rightarrow 0 as t\rightarrow\infty, and \tilde{r}(t, c)\rightarrow 0 as t\rightarrow\infty. By using I(t)\leq \tilde{I}(t), e(t, a)\leq \tilde{e}(t, a) and r(t, c)\leq \tilde{r}(t, c), we have I(t)\rightarrow 0, e(t, a)\rightarrow 0 and r(t, c)\rightarrow 0.

    Theorem 5.2. Assume R_{0}>1, the semiflow \{U(t)\}_{t\geq 0} generated by system (1) is uniformly persistent with respect to the pair (\partial M_{0}, M_{0}), that is there exists \varepsilon>0 , such that for each y\in M_{0},

    \displaystyle \lim\limits_{t\rightarrow +\infty}inf d(U(t)y, \partial M _{0})\geq \varepsilon.

    Furthermore, there exists a compact subset A_{0}\subset M _{0} which is a global attractor for \{U(t)\}_{t\geq 0} in M_{0}.

    Proof. Since the infection-free equilibrium E_{0}=(S^{0}, 0, 0, 0, 0_{L^{1}}, 0_{L^{1}}) is globally asymptotically stable in \partial M_{0}, applying Theorem 4.2 in Hale and Waltman [10], we only need to investigate the behavior of the solutions starting in M_{0} in some neighborhood of E_{0}. Then, we will show that W^{s}(\{E_{0}\})\bigcap M_{0}=\emptyset, where W^{s}(\{E_{0}\})=\{y\in X_{0+}:\lim_{t\rightarrow +\infty}U(t)y=E_{0}\}. Assume there exists y\in W^{s}(\{E_{0}\})\bigcap M_{0}, it follows that there exists t_{0}>0 such that I(t_{0})+\int_{0}^{\infty}e(t_{0}, a)da+\int_{0}^{\infty}r(t_{0}, c)dc > 0. Using the same argument in the proof of Lemma 3.6(i) in Demasse and Ducrot [6], we have that I(t)>0 for t\geq 0, and e(t, a)>0 for (t, a)\in [0, \infty)\times [0, \infty), r(t, c)>0 for(t, c)\in [0, \infty)\times [0, \infty). By means of the method of Brauer et al. [2], we define the following functions

    \begin{split} \displaystyle A(a)=&\int_{a}^{\infty}(\sigma(\theta)+B(0)\mu(\theta))e^{-\int_{a}^{\theta}\varepsilon(s)ds}d\theta, \\ \displaystyle B(c)=&\int_{c}^{\infty}k(\theta)e^{-\int_{c}^{\theta}\eta(s)ds}d\theta. \end{split} (46)

    For, \forall a, c>0, A(a), B(c)\geq 0, and A(0)=\theta_{1}+\theta_{2}\theta_{3}, B(0)=\theta_{2}. Furthermore, for \forall a\geq 0, c\geq 0 ,

    \begin{split} \displaystyle A'(a)&=-\sigma(a)+\varepsilon(a)A(a)-\theta_{2}\mu(a), \\ \displaystyle B'(c)&=-k(c)+\eta(c)B(c) . \end{split} (47)

    Consider the function

    \displaystyle\Phi(t)=I(t)+\int_{0}^{\infty}A(a)e(t, a)da+\int_{0}^{\infty}B(c)r(t, c)dc,

    which satisfies

    \displaystyle\frac{d\Phi(t)}{dt}=\beta I(\theta_{1}+\theta_{2}\theta_{3})(S-\frac{ S^{0}}{R_{0}}).

    Since y\in W^{s}(\{E_{0}\}), we have S(t)\rightarrow S^{0}, I(t)\rightarrow 0, as t\rightarrow \infty. When R_{0}>1, we know that the function \Phi(t) is not decreasing for t large enough. Thus there exists t_{0}>0 such that \Phi(t)\geq\Phi(t_{0}) for all t\geq t_{0}. Since \Phi(t_{0})>0, this prevents that the function (I(t), e(t, a), r(t, c)) converges to (0, 0_{L^{1}}, 0_{L^{1}}) as t\rightarrow \infty. A contradiction with S(t)\rightarrow S^{0}.


    6. Global asymptotic stability of the endemic equilibrium

    Let

    \displaystyle g(x)=x-1- \ln x,

    denote g'(x)=1-\frac{1}{x}. Thus, g: R^{+}\rightarrow R^{+} is concave up. Also, the function g has only one extremum which is a global minimum at 1, satisfying g(1)=0 and \forall x, y\in R, g(xy)\geq g(x)+g(y).

    Theorem 6.1. The unique endemic equilibrium E^{*} is globally asymptotically stable if R_0>1.

    Proof. Constructing the Lyapunov functional as follows

    \displaystyle V_{*}=W_{s}+W_{e}+W_{i}+W_{r},

    where

    \begin{split} \displaystyle W_{s}=&(\theta_{1}+\theta_{2}\theta_{3})S^{*}g(\frac{S}{S^{*}}), W_{i}=I^{*}g(\frac{I}{I^{*}}), \\ \displaystyle W_{e}=&\int_{0}^{+\infty}A(a)e^{*}(a)g(\frac{e(t, a)}{e^{*}(a)})da, W_{r}=\int_{0}^{+\infty}B(c)r^{*}(c)g(\frac{r(t, c)}{r^{*}(c)})dc. \end{split}

    Since \Lambda=bS^{*}+\beta S^{*}I^{*}, then the derivative of W_{s} along with the solutions of (1) is

    \begin{split} \frac{dW_{s}}{dt}=&(\theta_{1}+\theta_{2}\theta_{3})bS^{*}(2-\frac{S}{S^{*}}-\frac{S^{*}}{S})\\ &+(\theta_{1}+\theta_{2}\theta_{3})\beta S^{*}I^{*}(1-\frac{SI}{S^{*}I^{*}}-\frac{S^{*}}{S}+\frac{I}{I^{*}}). \end{split}

    Calculating the derivative of W_{e} along with the solutions of system (1) yields

    \begin{split} \displaystyle \frac{dW_{e}}{dt}=&\int_{0}^{+\infty}A(a)e^{*}(a)\frac{\partial}{\partial t}g(\frac{e(t, a)}{e^{*}(a)})da\\ \displaystyle =&\int_{0}^{+\infty}A(a)e^{*}(a)\frac{\partial}{\partial t}(\frac{e(t, a)}{e^{*}(a)}-1-\ln\frac{e(t, a)}{e^{*}(a)})da\\ \displaystyle =&\int_{0}^{+\infty}A(a)e^{*}(a)(\frac{1}{e^{*}(a)}-\frac{1}{e(t, a)})\frac{\partial}{\partial t}e(t, a)da\\ \displaystyle =&\int_{0}^{+\infty}A(a)e^{*}(a)(\frac{1}{e^{*}(a)}-\frac{1}{e(t, a)})(-\frac{\partial}{\partial a}e(t, a)-\varepsilon(a)e(t, a))da\\ \displaystyle =&-\int_{0}^{+\infty}A(a)e^{*}(a)(\frac{e(t, a)}{e^{*}(a)}-1)(\frac{e_{a}(t, a)}{e(t, a)}+\varepsilon(a))da. \end{split}

    Note that

    \begin{split} \displaystyle \frac{\partial}{\partial a}g(\frac{e(t, a)}{e^{*}(a)})=&\frac{e_{a}(t, a)+e(t, a)\varepsilon(a)}{e^{*}(a)}-\frac{e_{a}(t, a)}{e(t, a)}+\frac{e^{*}(a)(-\varepsilon(a))}{e^{*}(a)}\\ \displaystyle =&(\frac{e(t, a)}{e^{*}(a)}-1)(\frac{e_{a}(t, a)}{e(t, a)}+\varepsilon(a)). \end{split}

    And

    \begin{split} \displaystyle \frac{dA(a)}{da}=&A(a)\varepsilon(a)-\sigma(a)-\mu(a)B(0), \\ \displaystyle \frac{de^{*}(a)}{da}=&-\varepsilon(a)e^{*}(a). \end{split}

    Hence, using integration by parts, we have

    \begin{split} \displaystyle \frac{dW_{e}}{dt}=&-\int_{0}^{+\infty}A(a)e^{*}(a)\frac{\partial}{\partial a}g(\frac{e(t, a)}{e^{*}(a)})da\\ \displaystyle =&-A(a)e^{*}(a)g(\frac{e(t, a)}{e^{*}(a)})\mid_{0}^{+\infty}+\int_{0}^{+\infty}(\frac{d}{d a}A(a))e^{*}(a)g(\frac{e(t, a)}{e^{*}(a)})da \end{split}
    \begin{split}\displaystyle &+\int_{0}^{+\infty}A(a)(\frac{d}{da}e^{*}(a))g(\frac{e(t, a)}{e^{*}(a)})da\\ \displaystyle=&-A(a)e^{*}(a)g(\frac{e(t, a)}{e^{*}(a)})\mid_{+\infty}+A(0)e^{*}(0)g(\frac{e(t, 0)}{e^{*}(0)})\\ \displaystyle &-\int_{0}^{+\infty}e^{*}(a)(\sigma(a)+\mu(a) B(0))g(\frac{e(t, a)}{e^{*}(a)})da. \end{split}

    Note A(0)=\theta_{1}+\theta_{2}\theta_{3}, B(0)=\theta_{2}, e^{*}(0)=\beta S^{*}I^{*}+r_{2}I^{*}, e(t, 0)=\beta S(t)I(t)+r_{2}I(t), thus

    \begin{split} \displaystyle\frac{dW_{e}}{dt}=&-A(a)e^{*}(a)g(\frac{e(t, a)}{e^{*}(a)})\mid_{+\infty}+(\theta_{1}+\theta_{2}\theta_{3})(\beta S^{*}I^{*}+r_{2}I^{*})g(\frac{e(t, 0)}{e^{*}(0)})\\ \displaystyle&-\int_{0}^{+\infty}e^{*}(a)(\sigma(a)+\mu(a) \theta_{2})g(\frac{e(t, a)}{e^{*}(a)})da. \end{split}

    Further, it follows from (r_{1}+r_{2}+b+\delta_{i})I^{*}=\int_{0}^{+\infty}\sigma(a)e^{*}(a)da+\int_{0}^{+\infty}k(c)r^{*}(c)dc , that the derivative of W_{i} along with the solutions of system (1) gives

    \begin{split} \displaystyle\frac{dW_{i}}{dt}=&I^{*}(\frac{I_{t}}{I^{*}}-\frac{I_{t}}{I})\\ \displaystyle=&I^{*}(\frac{1}{I^{*}}-\frac{1}{I})[-(r_{1}+r_{2}+b+\delta_{i})I\\ \displaystyle&+\int_{0}^{+\infty}\sigma(a)e(t, a)da+\int_{0}^{+\infty}k(c)r(t, c)dc]\\ \displaystyle=&\int_{0}^{+\infty}\sigma(a)e^{*}(a)(1-\frac{I}{I^{*}}-\frac{I^{*}e(t, a)}{Ie^{*}(a)}+\frac{e(t, a)}{e^{*}(a)})da\\ \displaystyle&+\int_{0}^{+\infty}k(c)r^{*}(c)(1-\frac{I}{I^{*}}-\frac{I^{*}r(t, c)}{Ir^{*}(c)}+\frac{r(t, c)}{r^{*}(c)})dc. \end{split}

    Similar to W_{e}, by using B(0)=\theta_{2}, r^{*}(0)=r_{1}I^{*}+\int_{0}^{+\infty}\mu(a) e^{*}(a)da, and r(t, 0)=r_{1}I(t)+\int_{0}^{+\infty}\mu(a) e(t, a)da, the derivative of W_{r} along with the solutions of system (1) reads

    \begin{split} \displaystyle\frac{dW_{r}}{dt}&=\int_{0}^{+\infty}B(c)r^{*}(c)\frac{\partial}{\partial t}g(\frac{r(t, c)}{r^{*}(c)})dc\\ \displaystyle&=\int_{0}^{+\infty}B(c)r^{*}(c)\frac{\partial}{\partial t}[\frac{r(t, c)}{r^{*}(c)}-1-\ln\frac{r(t, c)}{r^{*}(c)}]dc\\ \displaystyle&=\int_{0}^{+\infty}B(c)r^{*}(c)[(\frac{1}{r^{*}(c)}-\frac{1}{r(t, c)})\frac{\partial}{\partial t}r(t, c)]dc\\ \displaystyle&=\int_{0}^{+\infty}B(c)r^{*}(c)[(\frac{1}{r^{*}(c)}-\frac{1}{r(t, c)})(-\frac{\partial}{\partial c}r(t, c)-\eta(c)r(t.c))]dc\\ \displaystyle&=-\int_{0}^{+\infty}B(c)r^{*}(c)(\frac{r(t, c)}{r^{*}(c)}-1)(\frac{r_{c}(t, c)}{r(t, c)}+\eta(c))dc. \end{split}

    Note

    \begin{split} \displaystyle \frac{\partial}{\partial c}g(\frac{r(t, c)}{r^{*}(c)}) =(\frac{r(t, c)}{r^{*}(c)}-1)(\frac{r_{c}(t, c)}{r(t, c)}+\eta(c)), \end{split}

    and

    \displaystyle \frac{dB(c)}{dc}=B(c)\eta(c)-k(c), \frac{dr^{*}(c)}{dc}=-\eta(c)r^{*}(c).

    Hence, using integration by parts, we have

    \begin{split} \displaystyle\frac{dW_{r}}{dt}=&-\int_{0}^{+\infty}B(c)r^{*}(c)\frac{\partial}{\partial c}g(\frac{r(t, c)}{r^{*}(c)})dc\\ \displaystyle =&-B(c)r^{*}(c)g(\frac{r(t, c)}{r^{*}(c)})\mid_{+\infty}+B(0)r^{*}(0)g(\frac{r(t, 0)}{r^{*}(0)})\\ \displaystyle &-\int_{0}^{+\infty}r^{*}(c)k(c)g(\frac{r(t, c)}{r^{*}(c)})dc\\ \displaystyle =&-B(c)r^{*}(c)g(\frac{r(t, c)}{r^{*}(c)})\mid_{+\infty}-\int_{0}^{+\infty}r^{*}(c)k(c)g(\frac{r(t, c)}{r^{*}(c)})dc\\ \displaystyle &+\theta_{2}(r_{1}I^{*}+\int_{0}^{+\infty}\mu(a) e^{*}(a)da)g(\frac{r(t, 0)}{r^{*}(0)}). \end{split}

    Note

    \begin{split} \displaystyle&\int_{0}^{+\infty}(\sigma(a)+\mu(a)\theta_{2})e^{*}(a)da=(\theta_{1}+\theta_{2}\theta_{3})(\beta S^{*}I^{*}+r_{2}I^{*}), \\ \displaystyle&\int_{0}^{+\infty}k(c)r^{*}(c)dc=\theta_{2}r_{1}I^{*}+\theta_{2}\theta_{3}(\beta S^{*}I^{*}+r_{2}I^{*}). \end{split}

    We derive

    \begin{split} \displaystyle\frac{dV_{*}}{dt}=&(\theta_{1}+\theta_{2}\theta_{3})bS^{*}(2-\frac{S}{S^{*}}-\frac{S^{*}}{S})-A(a)e^{*}(a)g(\frac{e(t, a)}{e^{*}(a)})\mid_{+\infty}\\ \displaystyle&-B(c)r^{*}(c)g(\frac{r(t, c)}{r^{*}(c)})\mid_{+\infty}+\int_{0}^{+\infty}\sigma(a)e^{*}(a)dag(\frac{e(t, 0)}{e^{*}(0)})\\ &-\int_{0}^{+\infty}\mu(a)\theta_{2}e^{*}(a)(g(\frac{e(t, a)}{e^{*}(a)})-g(\frac{e(t, 0)}{e^{*}(0)}))da\\ \displaystyle& +\int_{0}^{+\infty}k(c)r^{*}(c)dcg(\frac{r(t, 0)}{r^{*}(0)})+H_{1}+H_{2}+H_{3} \end{split}

    where

    \begin{split} \displaystyle H_{1}=&(\theta_{1}+\theta_{2}\theta_{3})\beta S^{*}I^{*}[-g(\frac{SI}{S^{*}I^{*}})-g(\frac{S^{*}}{S})+g(\frac{I}{I^{*}})]\\ \displaystyle\leq&(\theta_{1}+\theta_{2}\theta_{3})(\beta S^{*}I^{*}+r_{2}I^{*})[-g(\frac{SI}{S^{*}I^{*}})-g(\frac{S^{*}}{S})+g(\frac{I}{I^{*}})]\\ \displaystyle=&\int_{0}^{+\infty}(\sigma(a)+\mu(a)\theta_{2})e^{*}(a)[-g(\frac{SI}{S^{*}I^{*}})-g(\frac{S^{*}}{S})+g(\frac{I}{I^{*}})]da, \\ \displaystyle H_{2}=&\int_{0}^{+\infty}\sigma(a)e^{*}(a)(1-\frac{I}{I^{*}}-\frac{I^{*}e(t, a)}{Ie^{*}(a)}+\frac{e(t, a)}{e^{*}(a)}-\frac{e(t, a)}{e^{*}(a)}+1+\ln\frac{e(t, a)}{e^{*}(a)})da\\ \displaystyle=&-\int_{0}^{+\infty}\sigma(a)e^{*}(a)(g(\frac{I}{I^{*}})+g(\frac{I^{*}e(t, a)}{Ie^{*}(a)})da\\ \displaystyle\leq& -\int_{0}^{+\infty}\sigma(a)e^{*}(a)(g(\frac{I}{I^{*}})+g(\frac{e(t, a)}{e^{*}(a)})da, \end{split}
    \begin{split} \displaystyle H_{3}=&\int_{0}^{+\infty}k(c)r^{*}(c)(1-\frac{I}{I^{*}}-\frac{I^{*}r(t, c)}{Ir^{*}(c)}+\frac{r(t, c)}{r^{*}(c)}-\frac{r(t, c)}{r^{*}(c)}+1+\ln\frac{r(t, c)}{r^{*}(c)})dc\\ \displaystyle=&-\int_{0}^{+\infty}k(c)r^{*}(c)(g(\frac{I}{I^{*}})+g(\frac{I^{*}r(t, c)}{Ir^{*}(c)}))dc \end{split}
    \begin{split}\displaystyle\leq&-\int_{0}^{+\infty}k(c)r^{*}(c)(g(\frac{I}{I^{*}})+g(\frac{r(t, c)}{r^{*}(c)}))dc. \end{split}

    Hence, dV_{*}/dt \leq 0 holds true. Furthermore, the strict equality holds only if S=S^{*}, I=I^{*}, e(t, a)=e^{*}(a), r(t, c)=r^{*}(c). Consequently, the endemic equilibrium E^{*} of (1) is globally asymptotically stable if R_{0} > 1.


    7. Numerical simulations

    In the following, we provide some numerical simulations to illustrate the global stability of the disease-free equilibrium and the endemic equilibrium for system (1). We choose parameters \Lambda=3; b=0.065; n=0.02; \alpha_{1}=0.01; \alpha_{2}=0.03; r_{1}=0.1; r_{2}=0.2; and

    \sigma(a)=\left\{ \begin{array}{ccc} 0.3&a\geq \tau\\ 0&\tau\geq a\geq 0 \end{array}\right., k(c)=\left\{ \begin{array}{ccc} 0.1 &c\geq \tau\\ 0&\tau\geq c\geq 0 \end{array}\right., \mu(a)=\left\{ \begin{array}{ccc} 0.25 &a\geq \tau\\ 0&\tau\geq a\geq 0 \end{array}\right..

    Under the initial values

    \displaystyle S(0)=30, e(0, a)=6e^{-0.3a}, I(0)=10, r(0, c)=6e^{-0.3c}.

    In Figure 2, we choose \tau=12, then R_{0} < 1, while in Figure 3, we choose \tau=1, then R_{0}>1. The figures show the series of S(t) and I(t) which converge to their equilibrium values, and the age distribution and time series of e(t, a) and r(t, c), respectively.

    Figure 2. The time series of S(t) and I(t), and the age distributions of e(t, a) and r(t, c) when \tau=12.

    8. Discussion

    In this section, we briefly summarize our results. First, a PDE tuberculosis model (1) is proposed here to incorporate the latent-stage progression age of latent individuals and the relapse age of removed individuals. In addition, we assumed that infectious individuals might come into the latent class E due to incomplete treatment, and the relapse in the removed class. Under our assumptions, the expression of the basic reproduction number R_{0} is given, and we proved that if R_0 < 1 the disease-free equilibrium E_{0} is globally asymptotically stable, while if R_0>1 the unique endemic equilibrium E^{*} is globally asymptotically stable. Figure 2 and Figure 3 further verify our results.

    Figure 3. he time series of S(t) and I(t), and the age distributions of e(t, a) and r(t, c) when \tau=1.

    Acknowledgments

    The author is very grateful to Professors Shigui Ruan and Xingan Zhang for their supervision and assistance. The author also thanks the editor and the anonymous reviewers for their constructive comments that help to improve an early version of this paper.


    [1] [ R. A. Adams and J. J. Fournier, Sobolev Spaces, Academic Press, New York, 2003.
    [2] [ F. Brauer,Z. Shuai,P. van den Driessche, Dynamics of an age-of-infection cholera model, Mathematical Biosciences and Engineering, 10 (2013): 1335-1349.
    [3] [ Y. Chen,J. Yang,F. Zhang, The global stability of an SIRS model with infection age, Mathematical Biosciences and Engineering, 11 (2014): 449-469.
    [4] [ P. van den Driessche,L. Wang,X. Zou, Modeling diseases with latency and relapse, Mathematical Biosciences and Engineering, 4 (2007): 205-219.
    [5] [ P. van den Driessche,X. Zou, Modeling relapse in infectious diseases, Mathematical Biosciences, 207 (2007): 89-103.
    [6] [ R. D. Demasse,A. Ducrot, An age-structured within-host model for multistrain malaria infections, Siam Journal on Applied Mathematics, 73 (2013): 572-593.
    [7] [ O. Diekmann,J. A. P. Heesterbeek,J. A. J. Metz, On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations, Mathematical Biology, 28 (1990): 365-382.
    [8] [ J. K. Hale, Asymptotic Behavior of Dissipative Systems, American Mathematical Society, Providence, 1988.
    [9] [ J. K. Hale, Functional Differential Equations, Springer, New York, 1971.
    [10] [ J. K. Hale,P. Waltman, Persistence in infinite dimensional systems, Siam Journal on Mathematical Analysis, 20 (1989): 388-395.
    [11] [ F. Hoppensteadt, Mathematical Theories of Populations: Deomgraphics, Genetics, and Epidemics, SIAM, Philadelphia, 1975.
    [12] [ M. Iannelli, Mathematical theory of age-structured population dynamics, Applied Mathematics Monagraphs CNR, 7 (1994).
    [13] [ W. O. Kermack,A. G. McKendrick, A contribution to the mathematical theory of epidemics, Proc. R. Soc. London Ser. A, 115 (1927): 700-721.
    [14] [ W. O. Kermack,A. G. McKendrick, Contributions to the mathematical theory of epidemics Ⅱ. The problem of endemicity, Proc.R. Soc.London Ser. A, 138 (1932): 55-83.
    [15] [ W. O. Kermack,A. G. McKendrick, Contributions to the mathematical theory of epidemics Ⅲ. Further studies of the problem of endemicity, Proc.R. Soc.London Ser. A, 141 (1933): 94-122.
    [16] [ M. Y. Li,J. S. Muldowney, Global stability for the SEIR model in epidemilogy, Mathematical Biosciences, 125 (1995): 155-164.
    [17] [ L. Liu,J. Wang,X. Liu, Global stability of an SEIR epidemic model with age-dependent latency and relapse, Nonlinear Analysis, 24 (2015): 18-35.
    [18] [ P. Magal,C. C. McCluskey,G. F. Webb, Lyapunov functional and global asymptotic stability for an infection-age model, Applicable Analysis, 89 (2010): 1109-1140.
    [19] [ P. Magal, Compact attractors for time-periodic age-structured population models, Electronic Journal of Differential Equations, 65 (2001): 229-262.
    [20] [ C. C. McCluskey, Global stability for an SEI epidemiological model with continuous age-structure in the exposed and infectious classes, Mathematical Biosciences and Engineering, 9 (2012): 819-841.
    [21] [ A. G. McKendrick, Application of mathematics to medical problems, Proceedings of the Edinburgh Mathematical Society, 44 (1925): 98-130.
    [22] [ H. L. Smith and H. R. Thieme, Dynamical Systems and Population Persistence, American Mathematical Society, Providence, 2011.
    [23] [ C. Vargas-De-León, On the Global Stability of Infectious Diseases Models with Relapse, Abstr. Appl, 9 (2013): 50-61.
    [24] [ J. Wang,R. Zhang,T. Kuniya, A note on dynamics of an age-of-infection cholera model, Mathematical Biosciences and Engineering, 13 (2016): 227-247.
    [25] [ J. Wang,R. Zhang,T. Kuniya, Global dynamics for a class of age-infection HIV models with nonlinear infection rate, Mathematical Analysis and Applications, 432 (2015): 289-313.
    [26] [ J. Wang,R. Zhang,T. Kuniya, The dynamics of an SVIR epidemiological model with infection age, Applied Mathematics, 81 (2016): 321-343.
    [27] [ J. Wang,R. Zhang,T. Kuniya, The stability analysis of an SVEIR model with continuous age-structure in the exposed and infectious classes, Biological Dynamics, 9 (2015): 73-101.
    [28] [ L. Wang,Z. Liu,X. Zhang, Global dynamics for an age-structured epidemic model with media impact and incomplete vaccination, Nonlinear Analysis: Real World Applications, 32 (2016): 136-158.
    [29] [ J. Wang,S. S. Gao,X. Z. Li, A TB model with infectivity in latent period and imperfect treatment, Discrete Dynamics in Nature and Society, 2012 (2012): 267-278.
    [30] [ G. F. Webb, Theory of Nonlinear Age-Dependent Population Dynamics, Marcel Dekker, New York, 1985.
    [31] [ Y. Yang,S. Ruan,D. Xiao, Global stability of age-structured virus dynamics model with Beddinaton-Deangelis infection function, Mathematical Biosciences and Engineering, 12 (2015): 859-877.
    [32] [ J. H. Zhang, Qualitative Analysis and Data Simulation of Tuberculosis Transmission Models, Ph. D thesis, Central China Normal University in Wuhan, 2014.
    [33] [ Centers for Disease Control and Prevention, Vaccines and Immunizations-Measles Epidemiology and Prevention of Vaccine-Preventable Diseases Available from: http://www.cdc.gov/vaccines/pubs/pinkbook/meas.html.
  • This article has been cited by:

    1. Lili Liu, Xiaomei Feng, A multigroup SEIR epidemic model with age-dependent latency and relapse, 2018, 41, 01704214, 6814, 10.1002/mma.5193
    2. Y. Ma, C. R. Horsburgh, L. F. White, H. E. Jenkins, Quantifying TB transmission: a systematic review of reproduction number and serial interval estimates for tuberculosis, 2018, 146, 0950-2688, 1478, 10.1017/S0950268818001760
    3. Zhong-Kai Guo, Hong Xiang, Hai-Feng Huo, Analysis of an age-structured tuberculosis model with treatment and relapse, 2021, 82, 0303-6812, 10.1007/s00285-021-01595-1
    4. Yong Li, Xianning Liu, Yiyi Yuan, Jiang Li, Lianwen Wang, Global analysis of tuberculosis dynamical model and optimal control strategies based on case data in the United States, 2022, 422, 00963003, 126983, 10.1016/j.amc.2022.126983
    5. Zhong-Kai Guo, Hai-Feng Huo, Hong Xiang, Optimal control of TB transmission based on an age structured HIV-TB co-infection model, 2022, 359, 00160032, 4116, 10.1016/j.jfranklin.2022.04.005
    6. Shanjing Ren, Lingling Li, Global stability mathematical analysis for virus transmission model with latent age structure, 2022, 19, 1551-0018, 3337, 10.3934/mbe.2022154
    7. Dandan Sun, Yingke Li, Zhidong Teng, Tailei Zhang, Jingjing Lu, Dynamical properties in an SVEIR epidemic model with age‐dependent vaccination, latency, infection, and relapse, 2021, 44, 0170-4214, 12810, 10.1002/mma.7583
    8. Dipo Aldila, Joseph Páez Chávez, Karunia Putra Wijaya, Naleen Chaminda Ganegoda, Gracia Monalisa Simorangkir, Hengki Tasman, Edy Soewono, A tuberculosis epidemic model as a proxy for the assessment of the novel M72/AS01E vaccine, 2023, 120, 10075704, 107162, 10.1016/j.cnsns.2023.107162
    9. Xiaogang Liu, Yuming Chen, Xiaomin Li, Jianquan Li, Global stability of latency-age/stage-structured epidemic models with differential infectivity, 2023, 86, 0303-6812, 10.1007/s00285-023-01918-4
    10. Zhongkai Guo, Liang Zhang, Global Dynamics of an Age-Structured Tuberculosis Model with Vaccine Failure and Nonlinear Infection Force, 2023, 12, 2075-1680, 805, 10.3390/axioms12090805
    11. Zhong-Kai Guo, Hai-Feng Huo, Hong Xiang, Qiu-Yan Ren, Global dynamics of a tuberculosis model with age-dependent latency and time delays in treatment, 2023, 87, 0303-6812, 10.1007/s00285-023-01999-1
    12. Riya Das, Dhiraj Kumar Das, Tapan Kumar Kar, Qualitative analysis of TB transmission dynamics considering both the age since latency and relapse, 2023, 03784754, 10.1016/j.matcom.2023.09.021
    13. Zhong-Kai Guo, Hai-Feng Huo, Hong Xiang, TRANSMISSION DYNAMICS AND OPTIMAL CONTROL OF AN AGE-STRUCTURED TUBERCULOSIS MODEL, 2024, 14, 2156-907X, 1434, 10.11948/20230248
    14. Jun Zhang, Yasuhiro Takeuchi, Yueping Dong, Zhihang Peng, Modelling the preventive treatment under media impact on tuberculosis: A comparison in four regions of China, 2024, 24680427, 10.1016/j.idm.2024.02.006
    15. Tao-Li Kang, Hai-Feng Huo, Hong Xiang, Dynamics and optimal control of tuberculosis model with the combined effects of vaccination, treatment and contaminated environments, 2024, 21, 1551-0018, 5308, 10.3934/mbe.2024234
    16. Jinliang Wang, Guoyang Lyu, Analysis of an age‐space structured tuberculosis model with treatment and relapse, 2024, 0022-2526, 10.1111/sapm.12700
    17. Saduri Das, Prashant K. Srivastava, Pankaj Biswas, Tuberculosis transmission with multiple saturated exogenous reinfections, 2024, 17, 1793-5245, 10.1142/S179352452350064X
    18. Bashir Al‐Hdaibat, Mahmoud H. DarAssi, Irfan Ahmad, Muhammad Altaf Khan, Reem Algethamie, Investigating Tuberculosis Dynamics Under Various Control Strategies: A Comprehensive Analysis Using Real Statistical Data, 2025, 0170-4214, 10.1002/mma.10779
    19. Muhammad Altaf Khan, Mahmoud H. DarAssi, Irfan Ahmad, Noha Mohammad Seyam, Ebraheem Alzahrani, Modeling the Dynamics of Tuberculosis with Vaccination, Treatment, and Environmental Impact: Fractional Order Modeling, 2024, 141, 1526-1506, 1365, 10.32604/cmes.2024.053681
  • 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(3632) PDF downloads(698) Cited by(19)

Article outline

Figures and Tables

Figures(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog