Research article Special Issues

Estimating pinyon and juniper cover across Utah using NAIP imagery

  • Expansion of Pinus L. (pinyon) and Juniperus L. (juniper) (P-J) trees into sagebrush (Artemisia L.) steppe communities can lead to negative effects on hydrology, loss of wildlife habitat, and a decrease in desirable understory vegetation. Tree reduction treatments are often implemented to mitigate these negative effects. In order to prioritize and effectively plan these treatments, rapid, accurate, and inexpensive methods are needed to estimate tree canopy cover at the landscape scale. We used object based image analysis (OBIA) software (Feature AnalystTM for ArcMap 10.1®, ENVI Feature Extraction®, and Trimble eCognition Developer 8.2®) to extract tree canopy cover using NAIP (National Agricultural Imagery Program) imagery. We then compared our extractions with ground measured tree canopy cover (crown diameter and line point intercept) on 309 plots across 44 sites in Utah. Extraction methods did not consistently over- or under-estimate ground measured P-J canopy cover except where tree cover was >45%. Estimates of tree canopy cover using OBIA techniques were strongly correlated with estimates using the crown diameter method (r = 0.93 for ENVI, 0.91 for Feature AnalystTM, and 0.92 for eCognition). Tree cover estimates using OBIA techniques had lower correlations with tree cover measurements using the line-point intercept method (r = 0.85 for ENVI, 0.83 for Feature AnalystTM, and 0.83 for eCognition). All software packages accurately and inexpensively extracted P-J canopy cover from NAIP imagery when the imagery was not blurred, and when P-J cover was not mixed with Amelanchier alnifolia (Utah serviceberry) and Quercus gambelii (Gambel’s oak), which had similar spectral values as P-J.

    Citation: Darrell B. Roundy, April Hulet, Bruce A. Roundy, Ryan R. Jensen, Jordan B. Hinkle, Leann Crook, Steven L. Petersen. Estimating pinyon and juniper cover across Utah using NAIP imagery[J]. AIMS Environmental Science, 2016, 3(4): 765-777. doi: 10.3934/environsci.2016.4.765

    Related Papers:

    [1] Cruz Vargas-De-León, Alberto d'Onofrio . Global stability of infectious disease models with contact rate as a function of prevalence index. Mathematical Biosciences and Engineering, 2017, 14(4): 1019-1033. doi: 10.3934/mbe.2017053
    [2] Yoichi Enatsu, Yukihiko Nakata, Yoshiaki Muroya . Global stability for a class of discrete SIR epidemic models. Mathematical Biosciences and Engineering, 2010, 7(2): 347-361. doi: 10.3934/mbe.2010.7.347
    [3] Jinliang Wang, Hongying Shu . Global analysis on a class of multi-group SEIR model with latency and relapse. Mathematical Biosciences and Engineering, 2016, 13(1): 209-225. doi: 10.3934/mbe.2016.13.209
    [4] Rongjian Lv, Hua Li, Qiubai Sun, Bowen Li . Model of strategy control for delayed panic spread in emergencies. Mathematical Biosciences and Engineering, 2024, 21(1): 75-95. doi: 10.3934/mbe.2024004
    [5] Yu Ji . Global stability of a multiple delayed viral infection model with general incidence rate and an application to HIV infection. Mathematical Biosciences and Engineering, 2015, 12(3): 525-536. doi: 10.3934/mbe.2015.12.525
    [6] 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
    [7] Gonzalo Robledo . Feedback stabilization for a chemostat with delayed output. Mathematical Biosciences and Engineering, 2009, 6(3): 629-647. doi: 10.3934/mbe.2009.6.629
    [8] Pan Yang, Jianwen Feng, Xinchu Fu . Cluster collective behaviors via feedback pinning control induced by epidemic spread in a patchy population with dispersal. Mathematical Biosciences and Engineering, 2020, 17(5): 4718-4746. doi: 10.3934/mbe.2020259
    [9] Dongmei Li, Bing Chai, Weihua Liu, Panpan Wen, Ruixue Zhang . Qualitative analysis of a class of SISM epidemic model influenced by media publicity. Mathematical Biosciences and Engineering, 2020, 17(5): 5727-5751. doi: 10.3934/mbe.2020308
    [10] Jianquan Li, Zhien Ma, Fred Brauer . Global analysis of discrete-time SI and SIS epidemic models. Mathematical Biosciences and Engineering, 2007, 4(4): 699-710. doi: 10.3934/mbe.2007.4.699
  • Expansion of Pinus L. (pinyon) and Juniperus L. (juniper) (P-J) trees into sagebrush (Artemisia L.) steppe communities can lead to negative effects on hydrology, loss of wildlife habitat, and a decrease in desirable understory vegetation. Tree reduction treatments are often implemented to mitigate these negative effects. In order to prioritize and effectively plan these treatments, rapid, accurate, and inexpensive methods are needed to estimate tree canopy cover at the landscape scale. We used object based image analysis (OBIA) software (Feature AnalystTM for ArcMap 10.1®, ENVI Feature Extraction®, and Trimble eCognition Developer 8.2®) to extract tree canopy cover using NAIP (National Agricultural Imagery Program) imagery. We then compared our extractions with ground measured tree canopy cover (crown diameter and line point intercept) on 309 plots across 44 sites in Utah. Extraction methods did not consistently over- or under-estimate ground measured P-J canopy cover except where tree cover was >45%. Estimates of tree canopy cover using OBIA techniques were strongly correlated with estimates using the crown diameter method (r = 0.93 for ENVI, 0.91 for Feature AnalystTM, and 0.92 for eCognition). Tree cover estimates using OBIA techniques had lower correlations with tree cover measurements using the line-point intercept method (r = 0.85 for ENVI, 0.83 for Feature AnalystTM, and 0.83 for eCognition). All software packages accurately and inexpensively extracted P-J canopy cover from NAIP imagery when the imagery was not blurred, and when P-J cover was not mixed with Amelanchier alnifolia (Utah serviceberry) and Quercus gambelii (Gambel’s oak), which had similar spectral values as P-J.


    Infectious diseases (such as influenza, malaria, cholera, tuberculosis, hepatitis, AIDS, etc) have always seriously threatened humans' life and health. With in-depth understanding of infectious diseases, scientists have been continuing to explore effective methods to prevent and control the outbreaks of various infectious diseases. It is well known that mathematical models have played very important roles in analysis of control strategies for disease transmission [1,2,3,4,5,6,7,8,9,10].

    When studying the long-term evolutionary behavior of an ecological system, as pointed in [11], the equilibrium of biological system may not be the desirable one, and smaller value is required. This can be achieved by introducing suitable feedback control variable. The feedback control mechanism might be implemented through harvesting or culling procedures or certain biological control schemes [12]. In addition, in a control system, the time delay factor generally exists in the signal transmission process. Thus, feedback control with coupled time delay may have better biological significance [12] and has been extensively introduced into some important population ecological systems (see, for example, [13,14,15,16] and the references cited therein).

    In recent years, feedback control has also been successfully applied to some infectious disease dynamical systems. For example, in [17], the authors considered the following SI epidemic model with two feedback control variables:

    {˙S(t)=S(t)(raS(t)bI(t)c1u1(t)),˙I(t)=I(t)(bS(t)μfI(t)c2u2(t)),˙u1(t)=e1u1(t)+d1S(t),˙u2(t)=e2u2(t)+d2I(t). (1.1)

    In model (1.1), the state variables S(t) and I(t) represent the numbers of susceptibles and infectives at time t, respectively; u1(t) and u2(t) are feedback control variables. The number of susceptibles grows according to the regulation of a logistic curve with the capacity r/a (r>0, a>0) and a constant recruitment rate r; the constant b>0 is the transmission rate when susceptibles contact with infectives; the constants μ>0 and f>0 are the death rates of the infectives with respect to single and mutiple of infectives, respectively; the constants c1>0, c2>0, d1>0, d2>0, e1>0 and e2>0 are the feedback control parameters. By constructing suitable Lyapunov functions, the authors established threshold dynamics of model (1.1) completely determined by the threshold parameter γ0=(braμ)e1/(c1d1μ). The results in [17] indicate that, by appropriately choosing feedback control parameters, it can make the disease infection endemic or extinct. In [18], the author considered a two-group SI epidemic model with feedback control only in the susceptible individuals, and showed that the disease outbreaks can be controlled by adjusting feedback control parameters. In addition, in [19], the authors further extend model (1.1) to the case of patchy environment.

    Since the authors of [20] have introduced a nonlinear incidence rate g(I)S into classic Kermack-McKendrick SIR model, nonlinear incidence rate has been further introduced into more general SIR/SIRS epidemic models with time delays or infection age etc (see, for example, [21,22,23,24,25,26,27] and the references cited therein). Usually, the function g(I) takes the following two types: (ⅰ) saturated, such as g(I)=bI/(1+kI), or g(I)=bI2/(1+kI2); (ⅱ) unimodal, such as g(I)=bI/(1+kI2), here k>0 is constant. In biology, bI or bI2 measures the infection force of the disease, 1/(1+kI) or 1/(1+kI2) measures the inhibition effect from the behavioral change of the susceptible individuals when their number increases or from the crowding effect of the infective individuals [21,22].

    Recently, in [28], the authors further extended model (1.1) to the following more general case with the saturated incidence rate bSI/(1+kI) and feedback controls:

    {˙S(t)=S(t)(raS(t)bI(t)1+kI(t)c1u1(t)),˙I(t)=I(t)(bS(t)1+kI(t)μfI(t)c2u2(t)),˙u1(t)=e1u1(t)+d1S(t),˙u2(t)=e2u2(t)+d2I(t), (1.2)

    and some sufficient conditions for global asymptotic stability of the disease-free equilibrium and the endemic equilibrium of model (1.2) are established by the method of Lyapunov functions. In addition, the authors also considered permanence and existence of almost periodic solutions for a class of non-autonomous system based on model (1.2).

    Motivated by the above works and model (1.2), we further consider the following SI epidemic model with saturated incidence rate, two feedback control variables and four time delays:

    {˙S(t)=S(t)(raS(t)bI(t)1+kI(t)c1u1(tτ1)),˙I(t)=I(t)(bS(t)1+kI(t)μfI(t)c2u2(tτ2)),˙u1(t)=e1u1(t)+d1S(tτ3),˙u2(t)=e2u2(t)+d2I(tτ4). (1.3)

    The biological significance of all the parameters of model (1.3) are the same as in model (1.2) except time delays τi0 (i=1,2,3,4). In model (1.3), u1(t) and u2(t) are introduced as control variables. Usually, there always exist time delays in the transmission of information. Therefore, τ1 and τ2 can be understood as the result of transmission of information, and while τ3 and τ4 represent usual feedback control delays.

    The purpose in this paper focuses on global dynamics of the equilibria of model (1.3) by constructing appropriate Lyapunov functionals, and our results further extend and improve works in [17,28].

    The rest of the paper is organized as follows. In Section 2, we provide some preliminary results, including the well-posedness and dissipativeness of the solutions of model (1.3), the expression of the basic reproduction number and the classification of the equilibria of model (1.3). In Sections 3 and 4, we establish some sufficient conditions for global asymptotic stability and global attractivity of the disease-free equilibrium and the endemic equilibrium of model (1.3), which are the main results of this paper. In the last section, the conclusions and some numerical simulations are given.

    Let C+=C([τ,0],R4+) be the Banach space of continuous functions mapping the interval [τ,0] into R4+ equipped with the supremum norm, where τ=max{τ1,τ2,τ3,τ4}. The initial condition of model (1.3) is given as follows,

    S(θ)=ϕ1(θ),I(θ)=ϕ2(θ),u1(θ)=ϕ3(θ),u2(θ)=ϕ4(θ),θ[τ,0], (2.1)

    where ϕ=(ϕ1(θ),ϕ2(θ),ϕ3(θ),ϕ4(θ))C+.

    By using the standard theory of delay differential equations (DDEs) (see, for example, [29,30,31]), we can easily establish the following result.

    Theorem 2.1. The solution (S(t),I(t),u1(t),u2(t)) of model (1.3) with the initial condition (2.1) is existent, unique and nonnegative on [0,), and satisfies

    lim supt+S(t)ra,lim supt+I(t)rbaf,lim supt+u1(t)rd1ae1,lim supt+u2(t)rbd2afe2. (2.2)

    Moreover, the following bounded set

    Ω:={ϕC+:ϕ1ra,ϕ2rbaf,ϕ3rd1ae1,ϕ4rbd2afe2}

    is positively invariant with respect to model (1.3).

    Proof. It is not difficult to show that the solution (S(t),I(t),u1(t),u2(t)) of model (1.3) with the initial condition (2.1) is existent, unique and nonnegative on [0,). Let us consider ultimate boundedness of model (1.3). According to the first equation of model (1.3), we have that for t0,

    ˙S(t)S(t)(raS(t)), (2.3)

    which implies lim supt+S(t)ra. For any sufficiently small ε>0, there exists a ˆt>0 such that S(t)<ra+ε for tˆt. Further, according to the second equation of model (1.3), we have for tˆt,

    ˙I(t)I(t)[b(ra+ε)fI(t)],

    which implies

    lim supt+I(t)rbaf+bfε. (2.4)

    Since inequality (2.4) holds for arbitrary ε>0, we obtain lim supt+I(t)rbaf. Similarly, according to the last two equations of model (1.3), we can obtain lim supt+u1(t)rd1ae1, lim supt+u2(t)rbd2afe2.

    Let (S(t),I(t),u1(t),u2(t)) be the solution of model (1.3) with the initial function ϕ=(ϕ1,ϕ2,ϕ3,ϕ4)Ω. { For t0, we have ˙S(t)S(t)(raS(t)), which implies that for t0,

    S(t)raϕ1(0)ϕ1(0)+[raϕ1(0)]ertra,

    where ϕ1(0)ra is used. Further combining the second equation of model (1.3), for t0, we have ˙I(t)I(t)(rbafI(t)), which implies that for t0,

    I(t)rbafϕ2(0)ϕ2(0)+[rbafϕ2(0)]erbatrbaf,

    where ϕ2(0)rbaf is used. Thus, for t0, we have ˙u1(t)rd1ae1u1(t),˙u2(t)rbd2afe2u2(t). This implies that for t0,

    u1(t)rd1ae1+[ϕ3(0)rd1ae1]ee1trd1ae1,u2(t)rbd2afe2+[ϕ4(0)rbd2afe2]ee2trbd2afe2,

    where ϕ3(0)rd1ae1 and ϕ4(0)rbd2afe2 are used. Hence, it has that Ω is positively invariant with respect to model (1.3).

    The proof is completed.

    Obviously, model (1.3) always has a trivial equilibrium ˜E=(0,0,0,0) and a disease-free equilibrium E0=(S0,0,u01,0), where

    S0=re1ae1+d1c1,u01=rd1ae1+d1c1.

    Then, by the methods in [32,33], we can derive the expression of the basic reproduction number of model (1.3) as follows. First we define matrices F and V as

    F=(bS0000),V=(μ0d2e2).

    Then the basic reproduction number R0 is defined as the spectral radius of FV1. Therefore,

    R0:=ρ(FV1)=bS0μ=bre1μ(ae1+d1c1).

    Suppose (S,I,u1,u2) is an endemic equilibrium (positive equilibrium) of model (1.3), where S>0, I>0, u1>0, u2>0 satisfy the following equations

    raSbI1+kIc1u1=0,bS1+kIμfIc2u2=0,e1u1+d1S=0,e2u2+d2I=0. (2.5)

    From Eq (2.5), it is not difficult to obtain the following relationships

    u2=d2e2I,u1=d1e1S,S=e1ae1+d1c1(rbI1+kI). (2.6)

    Through Eq (2.6) and combining the second equation of Eq (2.5), we can obtain that I satisfies the following equation,

    F(I)be1(ae1+d1c1)(1+kI)(rbI1+kI)μ(f+d2c2e2)I=0.

    According to Eq (2.6), in order to ensure that S>0, we need to consider the following two cases:

    (ⅰ) rk<b, 0<I<rbrk˜I;

    (ⅱ) rkb, I>0.

    Clearly, for both case (ⅰ) and case (ⅱ), we have that

    ˙F(I)=be1k(ae1+d1c1)(1+kI)2(rbI1+kI)b2e1(ae1+d1c1)(1+kI)3(f+d2c2e2)<0.

    Hence, F(I) is monotonically decreasing with respect to I and

    limI0+F(I)=bre1ae1+d1c1μ=μ(R01).

    If R01, then limI0+F(I)0 and F(I)=0 has no positive roots. If R0>1, then limI0+F(I)>0.

    For case (ⅰ), we have that

    limI˜IF(I)=μ(f+d2c2e2)˜I<0.

    For case (ⅱ), we have that

    F(rbe1e2(ae1+d1c1)(fe2+d2c2))<rbe1ae1+d1c1(f+d2c2e2)rbe1e2(ae1+d1c1)(fe2+d2c2)=0.

    Therefore, for cases (ⅰ) and (ⅱ), F(I)=0 has a unique positive root I=I if R0>1.

    From the above discussions, we have the following result.

    Theorem 2.2. The following statements are true.

    (i) Model (1.3) always has a trivial equilibrium ˜E=(0,0,0,0).

    (ii) Model (1.3) always has a disease-free equilibrium E0=(S0,0,u01,0).

    (iii) Only for R0>1, model (1.3) has a unique endemic equilibrium E=(S,I,u1,u2), where

    S=e1ae1+d1c1(rbI1+kI),u1=d1e1S,u2=d2e2I,

    and I is the unique positive root of the equation F(I)=0.

    Remak 2.1. In fact, the classifications of the equilibria of models (1.2) and (1.3) are exactly the same for any τi0. Clearly, comparing with the reference [28], our Theorem 2.2 gives more complete classification of the equilibria of model (1.2) and clearer expression of the basic reproduction number R0.

    Let ¯E=(¯S,¯I,¯u1,¯u2) be any equilibrium of model (1.3). In order to investigate local stability of the equilibrium ¯E, we easily have that the characteristic equation of the corresponding linearized system of model (1.3) at ¯E is given by

    |λ(r2a¯Sb¯I1+k¯Ic1¯u1)b¯S(1+k¯I)2c1¯Seλτ10b¯I1+k¯Iλ(b¯S(1+k¯I)2μ2f¯Ic2¯u2)0c2¯Ieλτ2d1eλτ30λ+e100d2eλτ40λ+e2|=0. (3.1)

    At the trivial equilibrium ˜E=(0,0,0,0), the characteristic Eq (3.1) becomes

    (λr)(λ+μ)(λ+e1)(λ+e2)=0,

    which has a positive real root λ=r. Hence, ˜E is unstable for any τi0 (i=1,2,3,4).

    At the disease-free equilibrium E0=(S0,0,u01,0), the characteristic Eq (3.1) becomes

    (λ+μbS0)(λ+e2)[λ2+(aS0+e1)λ+ae1S0+d1c1S0eλ(τ1+τ3)]=0. (3.2)

    It is clear that Eq (3.2) has two real roots λ1=e2<0 and λ2=μ+bS0=μ(R01). Obviously, when R0>1, Eq (3.2) has a positive real root λ2>0, and hence, E0 is unstable for any τi0 (i=1,2,3,4). When R0<1, then λ2<0. When R0=1, then λ2=0 is a simple root of Eq (3.2).

    Let

    F1(λ,τ1,τ3)λ2+(aS0+e1)λ+ae1S0+d1c1S0eλ(τ1+τ3)=0. (3.3)

    The distribution of the roots of Eq (3.3) in the complex plane has been discussed in detail in [30,31,34]. Therefore, we have the following conclusions.

    Lemma 3.1. The following statements are true.

    (i) If d1c1ae1, then all the roots of Eq (3.3) have negative real parts.

    (ii) If d1c1>ae1, then all the roots of Eq (3.3) have negative real parts for τ1+τ3<τ013, and Eq (3.3) has at least one root which has positive real part for τ1+τ3>τ013, where

    τ013=1ωarccos[ω2ae1S0d1c1S0],
    ω=([(aS0)2+e21]+[(aS0)2+e21]24(S0)2(ae1+d1c1)(ae1d1c1)2)12.

    According to the discussions above and Lemma 3.1, it follows from stability theory and Hopf bifurcation theorem for DDEs (see, for example, [29,30,31]) that the following results hold.

    Theorem 3.1. The trivial equilibrium ˜E of model (1.3) is unstable for any τi0(i=1,2,3,4).

    Theorem 3.2. For any τ20 and τ40, the following statements are true.

    (i) If R0>1, then the disease-free equilibrium E0 is unstable for any τ10 and τ30.

    (ii) Assume that d1c1ae1. If R0<1, then the disease-free equilibrium E0 is locally asymptotically stable for any τ10 and τ30; If R0=1, then the disease-free equilibrium E0 is linearly stable for any τ10 and τ30.

    (iii) Assume that d1c1>ae1. If R0<1, then the disease-free equilibrium E0 is locally asymptotically stable for τ1+τ3<τ013, and is unstable for τ1+τ3>τ013. Moreover, model (1.3) undergoes a Hopf bifurcation at the disease-free equilibrium E0 when τ1+τ3=τ013.

    Remak 3.1. Theorem 3.2 indicates that time delays τ2 and τ4 do not affect local asymptotic stability of the disease-free equilibrium E0, and under the condition of d1c1ae1, time delays τ1 and τ3 also do not affect local asymptotic stability of E0. But under the condition d1c1>ae1, for larger time delay τ1 or τ3, stability of E0 will be lost.

    In the following discussions, we establish some sufficient conditions for global asymptotic stability of the disease-free equilibrium E0.

    Theorem 3.3. Assume that d1c1ae1. For any τi0 (i=1,2,3,4), the following statements are true.

    (i) If R0<1, then the disease-free equilibrium E0 is globally asymptotically stable in Ω1:={ϕΩ:ϕ1(0)>0}.

    (ii) If R0=1, then the disease-free equilibrium E0 is globally attractive in Ω1.

    Proof. First, it is easy to show that the set Ω1 is positively invariant for model (1.3). If R0<1, by Theorem 3.2, we only need to show that the disease-free equilibrium E0 is globally attractive.

    Define a Lyapunov functional L1 on Ω1 as follows,

    L1=V1+a20τ3(ϕ1(ξ)S0)2dξ+c1e12d10τ1(ϕ3(ξ)u01)2dξ,

    where

    V1=ϕ1(0)S0S0lnϕ1(0)S0+ϕ2(0)+c12d1(ϕ3(0)u01)2.

    It is clear that L1 is continuous on Ω1 and satisfies the condition (ii) of Lemma 3.1 in [35] on Ω1=¯Ω1Ω1.

    Calculating the derivative of L1 along any solution (S(t),I(t),u1(t),u2(t)) of model (1.3), it follows that, for t0,

    dL1dt=dV1dt+a2(S(t)S0)2a2(S(tτ3)S0)2+c1e12d1(u1(t)u01)2c1e12d1(u1(tτ1)u01)2, (3.4)

    where

    dV1dt=(S(t)S0)[a(S0S(t))bI(t)1+kI(t)+c1(u01u1(tτ1))]+I(t)[bS(t)1+kI(t)μfI(t)c2u2(tτ2)]+c1d1(u1(t)u01)[e1(u1(t)u01)+d1(S(tτ3)S0)]=a(S(t)S0)2[μbS01+kI(t)]I(t)fI2(t)c2I(t)u2(tτ2)c1e1d1(u1(t)u01)2+c1(S(t)S0)(u01u1(tτ1))+c1(u1(t)u01)(S(tτ3)S0), (3.5)

    here r=aS0+c1u01 and e1u01=d1S0 are used. Using the following inequality of arithmetic and geometric means,

    Λ1c1(S(t)S0)(u01u1(tτ1))+c1(u1(t)u01)(S(tτ3)S0)d1c1ae1[a2(S(t)S0)2+c1e12d1(u1(tτ1)u01)2+a2(S(tτ3)S0)2+c1e12d1(u1(t)u01)2],

    we further have that

    dL1dta2(1d1c1ae1)[(S(t)S0)2+(S(tτ3)S0)2][μbS01+kI(t)]I(t)fI2(t)c2I(t)u2(tτ2)c1e12d1(1d1c1ae1)[(u1(t)u01)2+(u1(tτ1)u01)2]. (3.6)

    Note that, if R01, we have that

    [μbS01+kI(t)]I(t)=[μbS0+μkI(t)1+kI(t)]I(t)=[μ(1R0)1+kI(t)+μkI(t)1+kI(t)]I(t)0. (3.7)

    Assume that d1c1<ae1 and R01. By inequalities (3.6) and (3.7), we can obtain that dL1dt0 for t0. Let M be the largest invariant set in the following set G:

    G:={ϕ¯Ω1:L1<anddL1dt=0}.

    Then, it follows from inequalities (3.6) and (3.7) that

    G{ϕΩ1:ϕ1(0)=S0,ϕ2(0)=0,ϕ3(0)=u01}.

    We can easily have from model (1.3) and the invariance of M that M={E0}. Therefore, it follows from Lemma 3.1 in [35] that the disease-free equilibrium E0 is globally attractive.

    Assume that d1c1=ae1 and R01. In this case, it is not easy to conclude that the largest invariant set M is the singleton {E0}. Hence, it is necessary to analyze Eq (3.4).

    Note that

    a2(S(t)S0)2c1e12d1(u1(tτ1)u01)2+c1(S(t)S0)(u01u1(tτ1))=d1c12e1(S(t)S0+e1d1(u1(tτ1)u01))2,a2(S(tτ3)S0)2c1e12d1(u1(t)u01)2+c1(u1(t)u01)(S(tτ3)S0)=d1c12e1(S(tτ3)S0e1d1(u1(t)u01))2=d1c12e1(S(tτ3)e1d1u1(t))2.

    Hence, Eq (3.4) can be rewritten as

    dL1dt=d1c12e1[(S(t)S0+e1d1(u1(tτ1)u01))2+(S(tτ3)e1d1u1(t))2][μ(1R0)1+kI(t)+μkI(t)1+kI(t)]I(t)fI2(t)c2I(t)u2(tτ2). (3.8)

    By inequality (3.7) and Eq (3.8), we can obtain that dL1dt0 for t0. Let M1 be the largest invariant set in the following set G1:

    G1:={ϕ¯Ω1:L1<anddL1dt=0}.

    Then, it follows from Eq (3.8) that

    G1{ϕΩ1:ϕ1(0)+e1d1ϕ3(τ1)=S0+e1d1u01,ϕ1(τ3)e1d1ϕ3(0)=0,ϕ2(0)=0}.

    For any φ=(φ1,φ2,φ3,φ4)M1, let (S(t),I(t),u1(t),u2(t)) be the solution of model (1.3) with the initial function φ. From the invariance of M1, we have that (St,It,u1t,u2t)M1G1 for any tR. Obviously, I(t)0 for any tR, and then from the fourth equation of model (1.3) and the invariance of M1, it is not difficult to obtain u2(t)0 for any tR. In addition, according to the first and third equations of model (1.3), we can obtain, for any tR,

    ˙S(t)=S(t)(raS(t)c1u1(tτ1))=S(t)(rd1c1e1S(t)c1u1(tτ1))=0,˙u1(t)=e1u1(t)+d1S(tτ3)=0.

    Thus, there exist constants δ1 and δ2 such that S(t)δ1 and u1(t)δ2 for any tR. It is not difficult to find that δ1 and δ2 satisfy

    δ1+e1d1δ2=S0+e1d1u01,δ1e1d1δ2=0,

    which imply that δ1=S0 and δ2=u01. Hence, S(t)S0 and u1(t)u01 for any tR. This shows that M1={E0}. Then, it follows from Lemma 3.1 in [35] that the disease-free equilibrium E0 is globally attractive.

    The proof is completed.

    Remak 3.2. Note the conclusion (ii) of Theorem 3.2, where we see that Theorem 3.3 gives complete conclusion of the global dynamics of the disease-free equilibrium E0 in the case of d1c1ae1.

    Now, we continue to discuss global dynamics of the disease-free equilibrium E0 in the absence of condition d1c1ae1. The following lemmas will be used.

    Lemma 3.2. (Barbalat's lemma [37,36]) Let x(t) be a real valued differentiable function defined on some half line [a,+), a(,+). If

    (i) limt+x(t)=α; |α|<.

    (ii) ˙x(t) is uniformly continuous for t>a.

    Then limt+˙x(t)=0.

    Lemma 3.3. Let (S(t),I(t),u1(t),u2(t)) be any solution of model (1.3) with the initial condition (2.1), then the following statements are true.

    (i) If rbaμ (which implies R0<1), then limt+I(t)=0, limt+u2(t)=0.

    (ii) If rb>aμ, then lim supt+I(t)IM, where

    IM={(kμ+f)2+4fk(braμ)(kμ+f)2fk>0,k>0,rbaμaf,k=0.

    Proof. By inequality (2.2), for arbitrary ε>0, there exists a T>0 such that S(t)<ra+ε for t>T. Then, it follows from the second equation of model (1.3) that, for t>T,

    ˙I(t)I(t)(b(ra+ε)1+kI(t)μfI(t))=I(t)1+kI(t)[fkI2(t)+(μk+f)I(t)b(ra+ε)+μ]. (3.9)

    If rbaμ, by inequality (3.9), it follows that, for t>T,

    ˙I(t)I(t)1+kI(t)[(μk+f)I(t)bε],

    which implies

    lim supt+I(t)bεμk+f. (3.10)

    Since I(t)0 and inequality (3.10) holds for arbitrary ε>0, we obtain that limt+I(t)=0. Further, according to the last equation of model (1.3), limt+u2(t)=0 can be easily obtained.

    If rb>aμ, it follows from inequality (3.9) that, for t>T,

    ˙I(t){fkI(t)1+kI(t)(I(t)I1(ε))(I(t)I2(ε)),k>0,f(I(t)I2(ε)),k=0.

    where

    I1(ε)=(μk+f)(μk+f)2+4fk(bra+bεμ)2fk<0,k>0I2(ε)={(μk+f)+(μk+f)2+4fk(bra+bεμ)2fk>0,k>0,b(r+aε)aμaf,k=0.

    Similarly, it follows that

    lim supt+I(t)I2(ε),andlim supt+I(t)I2(0)=IM. (3.11)

    The proof is completed.

    For convenience, let us give the following conditions (H1), (H2) and (H3):

    (H1)c1(e1+2d1)2τ1+c1r2τ3<a.(H2)e12τ1+r2a(a+b+2c1)τ3<e1d1.(H3)rc1b2aτ3<f+μk.

    We can obtain the following result.

    Theorem 3.4. Assume that R01. For any τ20 and τ40, the following statements are true.

    (i) If rbaμ and conditions (H1)(H2) hold, then the disease-free equilibrium E0 is globally attractive in X1:={ϕC+:ϕ1(0)>0}.

    (ii) If rb>aμ and conditions (H1)(H3) hold, then the disease-free equilibrium E0 is globally attractive in X1.

    Proof. It is easy to show that X1 is positively invariant for model (1.3). Let (S(t),I(t),u1(t),u2(t)) be the solution of model (1.3) with any initial function ϕX1. By inequality (2.2), for any sufficiently small ε0>0, there exists a t1>0 such that S(t)<ra+ε0S1(ε0) for t>t1.

    We continue to analyze dV1dt given by Eq (3.5). From Eq (3.5) and inequality (3.7), we have that, for t0,

    dV1dt=a(S(t)S0)2μ(1R0)1+kI(t)I(t)μk1+kI(t)I2(t)fI2(t)c2I(t)u2(tτ2)c1e1d1(u1(t)u01)2+Λ1. (3.12)

    Note that, for t>τ1+τ3, Λ1 can be rewritten as

    Λ1=c1(S(t)S0)(u1(t)u1(tτ1))+c1(u1(t)u01)(S(tτ3)S(t)):=Λ2+Λ3,

    where

    Λ2=c1(S(t)S0)(u1(t)u1(tτ1))=c1(S(t)S0)ttτ1˙u1(ξ)dξ=c1(S(t)S0)ttτ1(e1(u1(ξ)u01)+d1(S(ξτ3)S0))dξ,Λ3=c1(u1(t)u01)(S(tτ3)S(t))=c1(u1(t)u01)ttτ3˙S(ξ)dξ=c1(u1(t)u01)ttτ3S(ξ)(a(S0S(ξ))bI(ξ)1+kI(ξ)+c1(u01u1(ξτ1)))dξ.

    In addition, for t>τ1+τ3, we have

    Λ2c1e12ttτ1((S(t)S0)2+(u1(ξ)u01)2)dξ+c1d12ttτ1((S(t)S0)2+(S(ξτ3)S0)2)dξ=c1(e1+d1)2τ1(S(t)S0)2+c1e12ttτ1(u1(ξ)u01)2dξ+c1d12ttτ1(S(ξτ3)S0)2dξ. (3.13)

    Similarly, for t>t1+τ1+τ3, we have

    Λ3c1S(ε0)|u1(t)u01|ttτ3(a|S0S(ξ)|+bI(ξ)1+kI(ξ)+c1|u01u1(ξτ1)|)dξc1aS(ε0)2ttτ3((u1(t)u01)2+(S0S(ξ))2)dξ+c1bS(ε0)2ttτ3((u1(t)u01)2+I2(ξ)(1+kI(ξ))2)dξ+c21S(ε0)2ttτ3((u1(t)u01)2+(u01u1(ξτ1))2)dξ=c1S(ε0)2(a+b+c1)τ3(u1(t)u01)2+c1aS(ε0)2ttτ3(S0S(ξ))2dξ+c1bS(ε0)2ttτ3I2(ξ)(1+kI(ξ))2dξ+c21S(ε0)2ttτ3(u01u1(ξτ1))2dξ, (3.14)

    here S(t)<S1(ε0) for t>t1 is used. For t>t1+τ1+τ3, let us define the following function L2,

    L2=6i=1Vi,

    where V1 is defined as in the proof of Theorem 3.3,

    V2=c1aS(ε0)2ttτ3tθ(S(ξ)S0)2dξdθ,V3=c1d12ttτ1tθ(S(ξτ3)S0)2dξdθ+c1d1τ12ttτ3(S(ξ)S0)2dξ,V4=c1bS(ε0)2ttτ3tθI2(ξ)(1+kI(ξ))2dξdθ,V5=c1e12ttτ1tθ(u1(ξ)u01)2dξdθ,V6=c21S(ε0)2ttτ3tθ(u1(ξτ1)u01)2dξdθ+c21S(ε0)2τ3ttτ1(u1(ξ)u01)2dξ.

    For t>t1+τ1+τ3, we calculate the derivatives of Vi(i=2,,6) as follows,

    dV2dt=c1aS(ε0)2[τ3(S(t)S0)2ttτ3(S(ξ)S0)2dξ], (3.15)
    dV3dt=c1d12[τ1(S(t)S0)2ttτ1(S(ξτ3)S0)2dξ], (3.16)
    dV4dt=c1bS(ε0)2[τ3I2(t)(1+kI(t))2ttτ3I2(ξ)(1+kI(ξ))2dξ], (3.17)
    dV5dt=c1e12[τ1(u1(t)u01)2ttτ1(u1(ξ)u01)2dξ], (3.18)
    dV6dt=c21S(ε0)2[τ3(u1(t)u01)2ttτ3(u1(ξτ1)u01)2dξ]. (3.19)

    Combining (3.12)–(3.19), for t>t1+τ1+τ3, we finally have

    dL2dt=6i=1dVidt[ac1(e1+2d1)2τ1c1aS(ε0)2τ3](S(t)S0)2μ(1R0)1+kI(t)I(t)+c1bS(ε0)2τ3I2(t)(1+kI(t))2μk1+kI(t)I2(t)fI2(t)c2I(t)u2(tτ2)[c1e1d1c1e12τ1c1S(ε0)2(a+b+2c1)τ3](u1(t)u01)2. (3.20)

    Let us show global attractivity of the disease-free equilibrium E0 in the following two cases.

    Case (ⅰ) rbaμ (which implies R0<1) and (H1)(H2) hold.

    From conditions (H1)(H2), it has that, for sufficiently small ε0>0, the inequalities

    ac1(e1+2d1)2τ1c1aS(ε0)2τ3>0,e1d1e12τ1S(ε0)2(a+b+2c1)τ3>0 (3.21)

    hold. Furthermore, from Lemma 3.3, we have that limtI(t)=0 and limtu2(t)=0. Thus, there exists a t2>0 such that, for t>t2,

    c1bS(ε0)τ3I(t)<μd1c1ae1+d1c1.

    In addition, note that

    1R0=1bre1μ(ae1+d1c1)=μ(ae1+d1c1)bre1μ(ae1+d1c1)d1c1ae1+d1c1,

    from which we have that, for t>t2,

    Λ4:=μ(1R0)1+kI(t)I(t)+c1bS(ε0)2τ3I2(t)(1+kI(t))2I(t)1+kI(t)[μd1c1ae1+d1c1c1bS(ε0)2τ3I(t)]μd1c12(ae1+d1c1)I(t)1+kI(t)0. (3.22)

    Hence, it follows from inequalities (3.20)–(3.22) that dL2dt0 for t>˜T:=max{t1+τ1+τ3,t2}. This indicates that, for t>˜T, the function L2(t) is monotonically decreasing and bounded. Thus, the limitation limt+L2(t) exists.

    In addition, by Theorem 2.1, it is not difficult to show that, for t>˜T, the second derivative L2(t) is also bounded. This implies that L2(t) is uniformly continuous for t>˜T. Therefore, it follows from Lemma 3.1 that limtL2(t)=0. Again from inequalities (3.20)–(3.22), it follows that

    limtS(t)=S0,limtu1(t)=u01.

    This shows that the disease-free equilibrium E0 is globally attractive.

    Case (ⅱ) rb>aμ and (H1)(H3) hold.

    Similarly, from condition (H3), we have that the inequality

    c1bS(ε0)2τ3<f+μk (3.23)

    holds for sufficiently small ε0>0. From Lemma 3.3, there exists a t3>0 such that I(t)<IM+ε0 for t>t3. Then, we have that, for t>t3,

    Λ5:=c1bS(ε0)2τ3I2(t)(1+kI(t))2μk1+kI(t)I2(t)fI2(t)=I2(t)(1+kI(t))2[f(1+kI(t))2+μk(1+kI(t))c1bS(ε0)2τ3]I2(t)(1+k(IM+ε0))2[f+μkc1bS(ε0)2τ3]0. (3.24)

    Hence, by inequalities (3.20), (3.21) and (3.24), we can also obtain dL2dt0 for t>ˆT:=max{t1+τ1+τ3,t3}. By the same arguments as in Case (i), we can obtain

    limtS(t)=S0,limtI(t)=0,limtu1(t)=u01.

    Furthermore, from the last equation of model (1.3), we can easily have limtu2(t)=0. This shows that the disease-free equilibrium E0 is globally attractive.

    The proof is completed.

    From Theorems 3.2 and 3.4, we have the following two corollaries.

    Corollary 3.1. Assume that R0<1, d1c1>ae1 and τ1+τ3<τ013. For any τ20 and τ40, the following statements are true.

    (i) If rbaμ and conditions (H1)(H2) hold, then the disease-free equilibrium E0 is globally asymptotically stable in X1.

    (ii) If rb>aμ and conditions (H1)(H3) hold, then the disease-free equilibrium E0 is globally asymptotically stable in X1.

    Corollary 3.2. Assume that R0<1 and τ1=τ3=0. For any τ20 and τ40, then the disease-free equilibrium E0 is globally asymptotically stable in X1.

    Remak 3.3. If τi=0(i=1,2,3,4), then model (1.3) reduces into model (1.2). Clearly, Theorems 3.2, 3.3 and 3.4 extend and improve Theorem 1 in [28]. Further, if k=0, model (1.2) becomes the model discussed in [17]. Hence, Theorems 3.2, 3.3 and 3.4 also include Theorem 2.1 in [17] as a special case.

    Theoretical analysis of the distribution of the characteristic roots of characteristic equation (3.1) at the endemic equilibrium E=(S,I,u1,u2) usually involves some complicated computations, since there are four time delays τi0(i=1,2,3,4) in model (1.3). However, the numerical simulations in Section 5 show that, each of time delays τi (i=1,2,3,4) can destroy stability of the endemic equilibrium E of model (1.3) by properly choosing parameters. Hence, it is natural to consider the following two problems.

    (ⅰ) Under what conditions, the time delays τi0(i=1,2,3,4) are harmless for global asymptotic stability of the endemic equilibrium E of model (1.3).

    (ⅱ) Under what conditions, the time delays τi0(i=1,2,3,4) may be harmful for global asymptotic stability of the endemic equilibrium E of model (1.3).

    In this section, we study the above problems (i)–(ii) by constructing suitable Lyapunov functionals.

    For convenience, let us denote the following conditions,

    (H4)τ1=τ2=τ3=τ4=0.(H5)d1c1ae1,τ10,τ30,τ2=τ4=0.(H6)d2c2fe2,τ20,τ40,τ1=τ3=0.(H7)d1c1ae1,d2c2fe2,τ10,τ20,τ30,τ40.

    For problem (ⅰ) above, we have the following result.

    Theorem 4.1. Assume that R0>1. If one of conditions (H4)(H7) holds, then the endemic equilibrium E is globally asymptotically stable in X2:={ϕC+:ϕi(0)>0,i=1,2}.

    Proof. It is easy to show that the set X2 is positively invariant for model (1.3). Since E is the equilibrium of model (1.3), the following equalities hold,

    raSbI1+kIc1u1=0,bS1+kIμfIc2u2=0,e1u1+d1S=0,e2u2+d2I=0. (4.1)

    Define a Lyapunov functional W1 on X2 as follows,

    W1=(1+kI)(ϕ1(0)SSlnϕ1(0)S)+ϕ2(0)IIlnϕ2(0)I+c1(1+kI)2d1(ϕ3(0)u1)2+c22d2(ϕ4(0)u2)2.

    It is clear that W1 is continuous on X2 and positive definite with respect to E.

    Calculating the derivative of W1 along any solution (S(t),I(t),u1(t),u2(t)) of model (1.3), it follows that, for t0,

    dW1dt=(1+kI)(S(t)S)[a(SS(t))+bI1+kIbI(t)1+kI(t)+c1(u1u1(tτ1))]+(I(t)I)[bS(t)1+kI(t)bS1+kI(t)+bS1+kI(t)bS1+kI+f(II(t))+c2(u2u2(tτ2))]+c1(1+kI)d1(u1(t)u1)[e1(u1(t)u1)+d1(S(tτ3)S)]+c2d2(u2(t)u2)[e2(u2(t)u2)+d2(I(tτ4)I)],

    here Eq (4.1) is used. Note that

    (1+kI)(S(t)S)(bI1+kIbI(t)1+kI(t))+(I(t)I)(bS(t)1+kI(t)bS1+kI(t))=0,

    from which we have that, for t0,

    dW1dt=a(1+kI)(S(t)S)2[f+bkS(1+kI(t))(1+kI)](I(t)I)2c1e1(1+kI)d1(u1(t)u1)2c2e2d2(u2(t)u2)2+(1+kI)Υ1+Π1, (4.2)

    where

    Υ1=c1(S(t)S)(u1u1(tτ1))+c1(S(tτ3)S)(u1(t)u1),Π1=c2(I(t)I)(u2u2(tτ2))+c2(I(tτ4)I)(u2(t)u2).

    We consider global asymptotic stability of the endemic equilibrium E in the following four cases.

    If condition (H4) holds, it has that Υ1=Π1=0 and dW1dt is negative definite with respect to E. Hence, the endemic equilibrium E is globally asymptotically stable (see, for example, [29,30]).

    If condition (H5) holds, we have that Π1=0. Let us consider another functional as follows,

    W2=W1+a(1+kI)20τ3(ϕ1(ξ)S)2dξ+c1e1(1+kI)2d10τ1(ϕ3(ξ)u1)2dξ.

    W2 is continuous on X2, positive definite with respect to E, and satisfies condition (ii) of Lemma 3.1 in [35] on X2=¯X2X2.

    Calculating the derivative of W2 along any solution (S(t),I(t),u1(t),u2(t)) of model (1.3), it follows that, for t0,

    dW2dt=dW1dt+a(1+kI)2[(S(t)S)2(S(tτ3)S)2]+c1e1(1+kI)2d1[(u1(t)u1)2(u1(tτ1)u1)2]. (4.3)

    By Eqs (4.2) and (4.3), and the following inequality of arithmetic and geometric means:

    Υ1d1c1ae1[a2(S(t)S)2+c1e12d1(u1(tτ1)u1)2+a2(S(tτ3)S)2+c1e12d1(u1(t)u1)2], (4.4)

    we have that, for t0,

    dW2dta(1+kI)2(1d1c1ae1)[(S(t)S)2+(S(tτ3)S)2]c1e1(1+kI)2d1(1d1c1ae1)[(u1(t)u1)2+(u1(tτ1)u1)2][f+bkS(1+kI(t))(1+kI)](I(t)I)2c2e2d2(u2(t)u2)2.

    It follows from condition (H5) that dW2dt0 for t0. Hence, the endemic equilibrium E is stable. Furthermore, dW3dt=0 implies I(t)=I and u2(t)=u2.

    Let M2 be the largest invariant set in the set

    Γ2:={ϕ¯X2:W2<anddW2dt=0}.

    Then, it follows that

    Γ2{ϕX2:ϕ2(0)=I,ϕ4(0)=u2}.

    From model (1.3) and the invariance of M2, we can easily get that M2={E}. Thus, it follows from Lemma 3.1 in [35] that the endemic equilibrium E is globally asymptotically stable.

    If condition (H6) holds, it follows that Υ1=0. Let us consider a functional as follows,

    W3=W1+f20τ4(ϕ2(ξ)I)2dξ+c2e22d20τ2(ϕ4(ξ)u2)2dξ.

    W3 is continuous on X2, positive definite with respect to E, and satisfies condition (ii) of Lemma 3.1 in [35] on X2=¯X2X2.

    Calculating the derivative of W3 along any solution (S(t),I(t),u1(t),u2(t)) of model (1.3), it follows that, for t0,

    dW3dt=dW1dt+f2[(I(t)I)2(I(tτ4)I)2]+c2e22d2[(u2(t)u2)2(u2(tτ2)u2)2]. (4.5)

    By Eqs (4.2) and (4.5), and the following inequality of arithmetic and geometric means:

    Π1d2c2fe2[f2(I(t)I)2+c2e22d2(u2(tτ2)u2)2+f2(I(tτ4)I)2+c2e22d2(u2(t)u2)2], (4.6)

    we have that, for t0,

    dW3dta(1+kI)(S(t)S)2c1e1(1+kI)d1(u1(t)u1)2f2(1d2c2fe2)[(I(t)I)2+(I(tτ4)I)2]bkS(1+kI(t))(1+kI)(I(t)I)2c2e22d2(1d2c2fe2)[(u2(t)u2)2+(u2(tτ2)u2)2].

    It follows from condition (H6) that dW3dt0 for t0. Furthermore, dW3dt=0 implies S(t)=S and u1(t)=u1. By the same arguments as in the situation of condition (H5), we can also show that the endemic equilibrium E is globally asymptotically stable.

    If condition (H7) holds, let us consider the following Lyapunov functional,

    W4=W1+a(1+kI)20τ3(ϕ1(ξ)S)2dξ+f20τ4(ϕ2(ξ)I)2dξ+c1e1(1+kI)2d10τ1(ϕ3(ξ)u1)2dξ+c2e22d20τ2(ϕ4(ξ)u2)2dξ.

    W4 is continuous on X2, positive definite with respect to E, and satisfies condition (ii) of Lemma 3.1 in [35] on X2=¯X2X2.

    Calculating the derivative of W4 along any solution (S(t),I(t),u1(t),u2(t)) of model (1.3), it follows that, for t0,

    dW4dt=dW1dt+a(1+kI)2[(S(t)S)2(S(tτ3)S)2]+f2[(I(t)I)2(I(tτ4)I)2]+c1e1(1+kI)2d1[(u1(t)u1)2(u1(tτ1)u1)2]+c2e22d2[(u2(t)u2)2(u2(tτ2)u2)2].

    Further, by Eq (4.2) and inequalities (4.4) and (4.6), we have that, for t0,

    dW4dta(1+kI)2(1d1c1ae1)[(S(t)S)2+(S(tτ3)S)2]f2(1d2c2fe2)[(I(t)I)2+(I(tτ4)I)2]bkS(1+kI(t))(1+kI)(I(t)I)2c1e1(1+kI)2d1(1d1c1ae1)[(u1(t)u1)2+(u1(tτ1)u1)2]c2e22d2(1d2c2fe2)[(u2(t)u2)2+(u2(tτ2)u2)2].

    It follows from condition (H7) that dW4dt0 for t0. Furthermore, if d1c1<ae1, then dW4dt=0 implies that S(t)=S and u1(t)=u1. By the same arguments as in the situation of condition (H5), we can show that the endemic equilibrium E is globally asymptotically stable. If d1c1=ae1, also by the same arguments as in the situation of d1c1=ae1 in Theorem 3.3, we can show that E is globally asymptotically stable.

    The proof is completed.

    Remak 4.1. In the situation of condition (H7), Theorem 4.1 indicates that the time delays τi0(i=1,2,3,4) are harmless for global asymptotic stability of the endemic equilibrium E.

    Now, let us consider problem (ⅱ). Note that R0=bre1μ(ae1+d1c1)>1 implies rb>aμ. Let us denote the following conditions,

    (H8)c1(e1+2d1)2τ1+c1r2τ3+c2bIM2(1+kI)(1+kIM)τ4<a.(H9)c2(e2+2d2)2τ2+c1br2aτ3+c22[fIM+bkSIM(1+kI)(1+kIM)]τ4<f+bkS(1+kI)(1+kIM).(H10)e12τ1+r2a(a+b1+kI+2c1)τ3<e1d1.(H11)e22τ2+12[bIM1+kIM+fIM+bkSIM(1+kI)(1+kIM)+2c2IM]τ4<e2d2.

    In conditions (H8)(H11) above, the definition of IM is given in Lemma 3.3 (ⅱ). We have the following result.

    Theorem 4.2. Assume that R0>1. If conditions (H8)(H11) hold, then the endemic equilibrium E is globally attractive in X2.

    Proof. Let (S(t),I(t),u1(t),u2(t)) be the solution of model (1.3) with any initial function ϕX2. From Theorem 2.1 and Lemma 3.3, for sufficient small ε1>0, there exists a T1>0 such that, for t>T1,

    S(t)<ra+ε1:=S1(ε1),I(t)<IM+ε1:=IM(ε1).

    Hence, from Eq (4.2), it has that, for t>T1,

    dW1dta(1+kI)(S(t)S)2[f+bkS(1+kIM(ε1))(1+kI)](I(t)I)2c1e1(1+kI)d1(u1(t)u1)2c2e2d2(u2(t)u2)2+(1+kI)|Υ1|+|Π1|. (4.7)

    For simplicity, denote

    A(ε1):=(f+bkS(1+kI)(1+kIM(ε1)))IM(ε1).

    Note that Υ1 and Π1 can be rewritten as

    Υ1=c1(S(t)S)(u1(t)u1(tτ1))+c1(S(tτ3)S(t))(u1(t)u1):=Υ2+Υ3,Π1=c2(I(t)I)(u2(t)u2(tτ2))+c2(I(tτ4)I(t))(u2(t)u2):=Π2+Π3.

    Let us give appropriate estimations on |Υ2|, |Υ3|, |Π2| and |Π3|.

    It follows that, t>τ1+τ2+τ3+τ4+T1,

    |Υ2|=|c1(S(t)S)ttτ1˙u1(ξ)dξ|=|c1(S(t)S)ttτ1(e1(u1(ξ)u1)+d1(S(ξτ3)S))dξ|c1e12ttτ1((S(t)S)2+(u1(ξ)u1)2)dξ+c1d12ttτ1((S(t)S)2+(S(ξτ3)S)2)dξ=c1(e1+d1)2τ1(S(t)S)2+c1e12ttτ1(u1(ξ)u1)2dξ+c1d12ttτ1(S(ξτ3)S)2dξ, (4.8)
    |Υ3|=|c1(u1(t)u1)ttτ3˙S(ξ)dξ|=|c1(u1(t)u1)ttτ3S(ξ)[a(SS(ξ))+bI1+kIbI(ξ)1+kI(ξ)+c1(u1u1(ξτ1))]dξ|=|c1(u1(t)u1)ttτ3S(ξ)[a(SS(ξ))+b(II(ξ))(1+kI)(1+kI(ξ))+c1(u1u1(ξτ1))]dξ|
    c1S1(ε1)|u1(t)u1|ttτ3(a|SS(ξ)|+b1+kI|II(ξ)|+c1|u1u1(ξτ1)|)dξc1aS1(ε1)2ttτ3((u1(t)u1)2+(SS(ξ))2)dξ+c1bS1(ε1)2(1+kI)ttτ3((u1(t)u1)2+(II(ξ))2)dξ+c21S1(ε1)2ttτ3((u1(t)u1)2+(u1u1(ξτ1))2)dξ=c1S1(ε1)2[a+b1+kI+c1]τ3(u1(t)u1)2+c1aS1(ε1)2ttτ3(SS(ξ))2dξ+c1bS1(ε1)2(1+kI)ttτ3(II(ξ))2dξ+c21S1(ε1)2ttτ3(u1u1(ξτ1))2dξ, (4.9)
    |Π2|=|c2(I(t)I)ttτ2˙u2(ξ)dξ|=|c2(I(t)I)ttτ2(e2(u2(ξ)u2)+d2(I(ξτ4)I))dξ|c2e22ttτ2((I(t)I)2+(u2(ξ)u2)2)dξ+c2d22ttτ2((I(t)I)2+(I(ξτ4)I)2)dξc2(e2+d2)2τ2(I(t)I)2+c2e22ttτ2(u2(ξ)u2)2dξ+c2d22ttτ2(I(ξτ4)I)2dξ. (4.10)
    |Π3|=|c2(u2(t)u2)ttτ4˙I(ξ)dξ|=|c2(u2(t)u2)ttτ4I(ξ)[bS(ξ)1+kI(ξ)bS1+kI+f(II(ξ))+c2(u2u2(ξτ2))]dξ|=|c2(u2(t)u2)ttτ4[bI(ξ)1+kI(ξ)(S(ξ)S)+bkSI(ξ)(1+kI)(1+kI(ξ))(II(ξ))+fI(ξ)(II(ξ))+c2I(ξ)(u2u2(ξτ2))]dξ|c2|u2(t)u2|ttτ4[bIM(ε1)1+kIM(ε1)|S(ξ)S|+A(ε1)|II(ξ)|+c2IM(ε1)|u2u2(ξτ2)|]dξc2bIM(ε1)2(1+kIM(ε1))ttτ4((u2(t)u2)2+(S(ξ)S)2)dξ+c22A(ε1)ttτ4((u2(t)u2)2+(I(ξ)I)2)dξ+c22IM(ε1)2ttτ4((u2(t)u2)2+(u2(ξτ2)u2)2)dξ=c22[bIM(ε1)1+kIM(ε1)+A(ε1)+c2IM(ε1)]τ4(u2(t)u2)2+c2bIM(ε1)2(1+kIM(ε1))ttτ4(S(ξ)S)2dξ+c22A(ε1)ttτ4(I(ξ)I)2dξ+c22IM(ε1)2ttτ4(u2(ξτ2)u2)2dξ. (4.11)

    For t>τ1+τ2+τ3+τ4+T1, let us define the following function,

    U=W1+(1+kI)U1+U2,

    where

    U1=c1d12[ttτ1tθ(S(ξτ3)S)2dξdθ+τ1ttτ3(S(ξ)S)2dξ]+c1aS1(ε1)2ttτ3tθ(S(ξ)S)2dξdθ+c1bS1(ε1)2(1+kI)ttτ3tθ(II(ξ))2dξdθ,+c21S1(ε1)2[ttτ3tθ(u1(ξτ1)u1)2dξdθ+τ3ttτ1(u1(ξ)u1)2dξ]+c1e12ttτ1tθ(u1(ξ)u1)2dξdθ,
    U2=c2d22[ttτ2tθ(I(ξτ4)I)2dξdθ+τ2ttτ4(I(ξ)I)2dξ]+c22A(ε1)ttτ4tθ(I(ξ)I)2dξdθ+c2bIM(ε1)2(1+kIM(ε1))ttτ4tθ(S(ξ)S)2dξdθ+c22IM(ε1)2[ttτ4tθ(u2(ξτ2)u2)2dξdθ+τ4ttτ2(u2(ξ)u2)2dξ]+c2e22ttτ2tθ(u2(ξ)u2)2dξdθ.

    Computing the derivatives of U1 and U2, we easily get that, for t>τ1+τ2+τ3+τ4+T1,

    dU1dt=[c1d12τ1+c1aS1(ε1)2τ3](S(t)S)2c1d12ttτ1(S(ξτ3)S)2dξc1aS1(ε1)2ttτ3(S(ξ)S)2dξ+c1bS1(ε1)2(1+kI)τ3(I(t)I)2c1bS1(ε1)2(1+kI)ttτ3(I(ξ)I)2dξ+[c1e12τ1+c21S1(ε1)2τ3](u1(t)u1)2c21S1(ε1)2ttτ3(u1(ξτ1)u1)2dξc1e12ttτ1(u1(ξ)u1)2dξ, (4.12)
    dU2dt=c2bIM(ε1)2(1+kIM(ε1))τ4(S(t)S)2c2bIM(ε1)2(1+kIM(ε1))ttτ4(S(ξ)S)2dξ+[c2d22τ2+c22A(ε1)τ4](I(t)I)2c2d22ttτ2(I(ξτ4)I)2dξc22A(ε1)ttτ4(I(ξ)I)2dξ+[c2e22τ2+c22IM(ε1)2τ4](u2(t)u2)2c22IM(ε1)2ttτ4(u2(ξτ2)u2)2dξc2e22ttτ2(u2(ξ)u2)2dξ. (4.13)

    In summary, by (4.7)–(4.13), we have that, for t>τ1+τ2+τ3+τ4+T1,

    dUdt=dW1dt+(1+kI)dU1dt+dU2dt(1+kI)[a(2d1+e1)c12τ1c1aS1(ε1)2τ3c2bIM(ε1)2(1+kI)(1+kIM(ε1))τ4](S(t)S)2[f+bkS(1+kIM(ε1))(1+kI)(2d2+e2)c22τ2c1bS1(ε1)2τ3c22A(ε1)τ4](I(t)I)2(1+kI)c1[e1d1e12τ1S1(ε1)2(a+b1+kI+2c1)τ3](u1(t)u1)2c2[e2d2e22τ212(bIM(ε1)1+kIM(ε1)+A(ε1)+2c2IM(ε1))τ4](u2(t)u2)2:=(1+kI)Q1(ε1)(S(t)S)2Q2(ε1)(I(t)I)2(1+kI)c1Q3(ε1)(u1(t)u1)2c2Q4(ε1)(u2(t)u2)2. (4.14)

    Furthermore, from conditions (H8)(H11), we see that, for sufficiently small ε1>0, we have that Qi(ε1)>0(i=1,2,3,4). This shows that dUdt0 for t>τ1+τ2+τ3+τ4+T1. By similar arguments as in the proof of Theorem 3.4, we can have

    limt+S(t)=S,limt+I(t)=I,limt+u1(t)=u1,limt+u2(t)=u2.

    This proves that the endemic equilibrium E is globally attractive.

    The proof is completed.

    Remak 4.2. If τi=0(i=1,2,3,4), then model (1.3) reduces into model (1.2). Clearly, Theorem 4.1 extends and improves Theorem 2 in [28]. Further, if k=0, Theorem 4.1 also include Theorem 2.2 in [17] as a special case.

    In this paper, we consider the SI epidemic model (1.3) with two feedback control variables and four time delays. In biology, model (1.3) has more general biological significance. Then, by skillfully constructing appropriate Lyapunov functionals, and combining Lyapunov–LaSalle invariance principle and Barbalat's lemma, some sufficient conditions for global dynamics of the equilibria of model (1.3) are established. In the case of d1c1ae1, Theorem 3.3 gives complete conclusion on global asymptotic stability of the disease-free equilibrium E0. In the case of d1c1>ae1, in Theorem 3.4, global attractivity of the disease-free equilibrium E0 is considered under conditions (H1)(H3). Note that, in the case, it has from Theorem 3.2 that local asymptotic stability of the disease-free equilibrium E0 also depends on the time delays τ1 and τ3. Hence, The set of conditions (H1)(H3) has certain rationality. Furthermore, as a special case, Theorems 3.2, 3.3 and 3.4 improve and generalize Theorem 1 in [28] and Theorem 2.1 in [17].

    We also establish some sufficient conditions for global dynamics of the endemic equilibrium E of model (1.3). Theorem 4.1 shows that, if condition (H5), or (H6), or (H7) holds, the time delays τ1 and τ3, or τ2 and τ4, or τi(i=1,2,3,4) are harmless for global asymptotic stability of the endemic equilibrium E of model (1.3). If condition (H4) holds, i.e., τi=0(i=1,2,3,4), we see that Theorem 4.1 includes Theorem 2 in [28] and Theorem 2.2 in [17] as a special case. Furthermore, in Theorem 4.2, under condition R0>1, a set of sufficient conditions (H8)(H11) is obtained to ensure global attractivity of the endemic equilibrium E of model (1.3). Note that the subsequent numerical simulations imply that any one of time delays τi(i=1,2,3,4) may destroy local asymptotic stability of the endemic equilibrium E of model (1.3) and result in the occurrence of periodic oscillations etc.. Hence, in the set of conditions (H8)(H11), it should be feasible to have some certain limits on the lengths of time delays τi(i=1,2,3,4).

    In the following, let us give some numerical simulations to summarize the applications of Theorem 3.4, Corollary 3.1 and Theorem 4.2.

    Firstly, let us choose the values of a set of parameters as follows,

    r=1.2,a=1,b=0.4,k=1,c1=0.8,μ=0.25,f=1,c2=1,e1=0.6,d1=1.2,e2=1,d2=1. (5.1)

    By computations, we have that 0.96=d1c1>ae1=0.6, 0.48=rb>aμ=0.25, R00.738<1, E0=(0.462,0,0.923,0) and τ0134.542. Furthermore, conditions (H1)(H3) become the following inequalities,

    (ˆH)1.2τ1+0.48τ3<1,0.3τ1+1.8τ3<0.5.

    It has from Theorem 3.2 that the disease-free equilibrium E0 is locally asymptotically stable for τ1+τ3<τ0134.542 and unstable for τ1+τ3>τ0134.542. If time delays τ1 and τ3 satisfy more restrictive condition (ˆH), it has from Theorem 3.4 or Corollary 3.1 that the disease-free equilibrium E0 is also globally asymptotically stable.

    Let us choose τ1=0.75 and τ3=0.15. We see that conditions τ1+τ3<τ0134.542 and (ˆH) are satisfied. Figure 1(a) shows that the disease-free equilibrium E0 is globally asymptotically stable.

    Figure 1.  The phase trajectories and solution curves of the model (1.3) with R00.738<1. (a) (H1)(H3) hold, and E0 is globally asymptotically stable. (b) (H1)(H3) do not hold, but τ1+τ3<τ013. E0 is locally asymptotically stable. (c) (H1)(H3) do not hold, but τ1+τ3>τ013 holds. E0 is unstable and periodic oscillations occur.

    Let us choose τ1=2 and τ3=1.25. We see that condition (ˆH) does not hold, but condition τ1+τ3<τ0134.542 is still satisfied. Figure 1(b) shows that the disease-free equilibrium E0 is locally asymptotically stable.

    Let us further choose τ1=2.75 and τ3=2. We see that condition τ1+τ3>τ0134.542 holds. Figure 1(c) shows that the disease-free equilibrium E0 becomes unstable and periodic oscillations occur. In Figure 1(a)–(c), for simplicity, τ2 and τ4 are fixed as τ2=3 and τ4=6, respectively.

    Secondly, let us give numerical simulations in the situation of the basic reproduction number R0>1.

    Let us choose the values of set of parameters as follows,

    r=2,a=0.8,b=1,k=1.5,c1=1.1,μ=0.25,f=0.8,c2=0.8,e1=0.8,d1=1,e2=0.5,d2=1.8. (5.2)

    By computations, we have that R03.678>1, E(0.870,0.130,1.087,0.467). Furthermore, conditions (H8)(H11) approximately become the following inequalities,

    (˜H){1.54τ1+1.1τ3+0.132τ4<0.8,1.64τ2+1.375τ3+0.481τ4<1.246,0.4τ1+4.796τ3<0.8,0.25τ2+1.570τ4<0.278.

    If time delays τi(i=1,2,3,4) satisfy condition (˜H), it has from Theorem 4.2 that the endemic equilibrium E is globally attractive.

    Let us choose τ1=0.4, τ2=0.35, τ3=0.1 and τ4=0.12. We see that condition (˜H) is satisfied. Figure 2(a) shows that the endemic equilibrium E is globally attractive.

    Figure 2.  The phase trajectories and solution curves of the model (1.3) with R03.678>1. (a) (H8)(H11) hold, and E is globally attractive. (b) (H8)(H11) do not hold, but E may be still attractive. (c) (H8)(H11) do not hold, and E is unstable and periodic oscillations occur.

    Let us choose larger values here, τ1=1.6, τ2=1.4, τ3=0.4 and τ4=0.48. We see that condition (˜H) does not hold. Figure 2(b) show that the endemic equilibrium E may be still attractive.

    Let us further choose τ1=2, τ2=1.75, τ3=0.6 and τ4=0.8. We see that condition (˜H) does not hold. Figure 2(c) shows that the endemic equilibrium E becomes unstable and periodic oscillations occur.

    Moreover, Figure 3(a)–(d) show that any one of time delays τi (i=1,2,3,4) can destroy local asymptotic stability of the endemic equilibrium E of model (1.3). Here, the values of the parameters of model (1.3) are the same as Eq (5.2).

    Figure 3.  The phase trajectories and solution curves of the model (1.3) with R03.678>1.

    At the end of the paper, in view of Figure 1(b) and Figure 2(b) above, we would like to point out that conditions (H1)(H3) in Theorem 3.3 and conditions (H8)(H11) in Theorem 4.2 are actually conservative and worth of further improving. In addition, for global asymptotic stability of the endemic equilibrium E of model (1.3), in the case of d1c1>a1e1 or d2c2>fe2, to give some sufficient conditions which are different from conditions (H8)(H11) may be also interesting, since Figure 4(a)–(d) show that model (1.3) may have richer dynamic behaviors. Here, the values of the parameters of model (1.3) are also the same as Eq (5.2). Further, it would be interesting to extend model (1.3) to the case of non-autonomous model and to consider the uniform persistence and existence of almost periodic solutions etc. [38,39].

    Figure 4.  The phase trajectories and solution curves of the model (1.3) with R03.678>1.

    This paper is supported by Beijing Natural Science Foundation (No.1202019) and National Natural Science Foundation of China (No.11971055). The authors would like to thank the reviewers and the editor for their careful reading, helpful comments and suggestions that greatly improved the paper.

    All authors declare no conflicts of interest in this paper.

    [1] Miller RF, Tausch RJ (2001) The role of fire in pinyon and juniper woodlands: a descriptive analysis. In: K. E. M. Galley and T. P. Wilson [eds]. Proceedings of the Invasive Species Workshop: the Role of Fire in the Control and Spread of Invasive Species. Fire Conference 2000: the First National Congress on Fire Ecology, Prevention, and Management. Miscellaneous Publication No. 11, Tall Timber Research Station, Tallahassee, FL, USA. p. 15-30.
    [2] Miller RF, Bates JD, Svejcar TJ, et al. (2005) Biology, ecology, and management of western juniper. Corvallis, OR, USA: Oregon State University, Technical Bulletin 152. 77 p.
    [3] Miller RF, Svejcar TJ, Rose JA (2000) Impacts of western juniper on plant community composition and structure. J Range Manage 53: 574-585.
    [4] Bates JD, Miller RF, Svejcar TJ (2005) Long-term successional trends following western juniper cutting. Rangeland Ecol Manag 58: 533-541.
    [5] Pierson FB, Williams CJ, Kormos PR, et al. (2010) Hydrologic vulnerability of sagebrush steppe following pinyon and juniper encroachment. Rangeland Ecol Manag 63: 614-629.
    [6] Roundy BA, Miller RF, Tausch RJ, et al. (2014) Understory cover responses to piñon-juniper treatments across tree cover gradients in the Great Basin. Rangeland Ecol Manag 67: 482-494. doi: 10.2111/REM-D-13-00018.1
    [7] Knick ST, Hanser SE, Leu M (2014) Ecological Scale of Bird Community Response to Pinyon-Juniper Removal. Rangeland Ecol Manag 67: 553-562. doi: 10.2111/REM-D-13-00023.1
    [8] Gruell GE (1999) Historical and modern roles of fire in pinyon-juniper. In: Steve B. Monsen and Richard Stevens [EDS]. Proceedings: Ecology and management of pinyon-juniper communities within the Interior West; 15-18 September 1997; Brigham Young University, Provo, Utah, USA. Ogden, UT, USA: U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station, RMRS-P-9. p. 24-28.
    [9] Miller RF, Ratchford J, Roundy BA, et al. (2014) Response of conifer-encroached shrublands in the Great Basin to prescribed fire and mechanical treatments. Rangeland Ecol Manag 67: 468-481. doi: 10.2111/REM-D-13-00003.1
    [10] Young KR, Roundy BA, Eggett DL (2014) Mechanical mastication of Utah juniper encroaching sagebrush steppe increases inorganic soil N. Appl Environ Soil Sci, 2014.
    [11] Young KR, Roundy BA, Bunting SC, et al. (2015) Utah juniper and two- needle piñon reduction alters fuel loads. Int J Wildland Fire 24: 236-248.
    [12] Hulet A, Roundy BA, Petersen SL, et al. (2013) Assessing the relationship between ground measurements and object-based image analysis of land cover classes in pinyon and juniper woodlands. Photogramm Eng Rem S 79: 799-808.
    [13] Hulet A, Roundy BA, Petersen SL, et al. (2013) An object-based image analysis of pinyon and juniper woodlands treated to reduce fuels. Environ Manag 53: 660-671.
    [14] Hulet A, Roundy BA, Petersen SL, et al. (2014) Cover estimations using object-based image analysis rule sets developed across multiple scales in pinyon and juniper woodlands. Rangeland Ecol Manag 67: 318-327.
    [15] Madsen MD, Zvirzdin DL, Davis BD, et al. (2011) Feature extraction techniques for measuring pińon and juniper tree cover and density, and comparison with field-based management surveys. Environ Manag 47: 766-776.
    [16] Hulet A, Roundy BA, Petersen SL, et al. (2014) Utilizing NAIP imagery to estimate tree cover and biomass of pinyon and juniper woodlands. Rangeland Ecol Manag 67: 563-572.
    [17] Weisberg PJ, Lingua E, Pillai RB (2007) Spatial patterns on pinyon-juniper woodland expansion in Central Nevada. Rangeland Ecol Manag 60: 115-124.
    [18] Bybee J, Roundy BA, Young KR, et al. (2016) Vegetation response to pinon and juniper tree shredding. Rangeland Ecol Manag 69: 224-234.
    [19] Booth DT, Cox SE, Fifield C, et al. (2005) Image analysis compared with other methods for measuring ground cover. Arid Land Res Manag 19: 91-100.
    [20] Ko D, Bristow N, Greenwood D, et al. (2009) Canopy cover estimation in semiarid woodlands: comparison of field-based and remote sensing methods. Forest Sci 55: 132-141.
    [21] Laliberte AS, Fredrickson EL, Rango A (2007) Combining decision trees with hierarchical object-oriented image analysis for mapping arid rangelands. Photogramm Eng Rem S 73: 197-207. doi: 10.14358/PERS.73.2.197
    [22] Laliberte AS, Rango A, Herrick JE, et al. (2007) An object-based image analysis approach for determining fractional cover of senescent and green vegetation with digital plot photography. J Arid Environ 69: 1-14. doi: 10.1016/j.jaridenv.2006.08.016
    [23] Laliberte AS, Rango A (2009) Texture and scale in object-based analysis of subdecimeter resolution unmanned aerial vehicle (UAV) imagery. IEEE T Geosci Remote 47: 761-770. doi: 10.1109/TGRS.2008.2009355
    [24] Davies KW, Petersen SL, Johnson DD, et al. (2010) Estimating juniper cover from national agriculture imagery program (NAIP) imagery and evaluating relationships between potential cover and environmental variables. Rangeland Ecol Manag 63: 630-637.
    [25] Visual Information Solutions (2008) ENVI Feature Extraction Module User’s Guide. Available from: http://www.exelisvis.com/portals/0/pdfs/envi/feature_extraction_module.pdf.
    [26] Blundell S, Opitz D, Morris M, et al. (2008) Feature Analyst V5.0. Paper presented at The ASPRS Annual Conference, Portland, Oregon, USA.
    [27] Trimble (2011) eCognition Developer 8.64.1 reference book, version 8.64.1, Trimble Germany GmbH, München, Germany.
    [28] McIver JD, Brunson M (2014) Multidisciplinary, multisite evaluation of alternative sagebrush steppe restoration treatments: the SageSTEP project. Rangeland Ecol Manag 67: 435-439.
    [29] Bybee J, Roundy BA, Young KR, Hulet A, Roundy DB, Crook L, Aanderud Z, Eggett D, Cline NL (2016) Vegetation response to pinon and juniper tree shredding. Rangeland Ecol Manag 69: 224-234.
    [30] Baatz M, Schäpe A (2000) Multiresolution segmentation—an optimization approach for high quality multi-scale image segmentation. In: J. Strobl, et al. [eds]. Angewandte Geographische Informationsverarbeitung XII. Wichmann, Heidelberg, Germany. p. 12-23.
    [31] Congalton RG (2001) Accuracy assessment and validation of remotely sensed and other spatial information. Int J Wildland Fire 10: 321-328.
    [32] Landis J, Koch G (1977) The measurement of observer agreement for categorical data. Biometrics 33: 159-174.
  • This article has been cited by:

    1. Manuel De la Sen, Asier Ibeas, Santiago Alonso-Quesada, On the Supervision of a Saturated SIR Epidemic Model with Four Joint Control Actions for a Drastic Reduction in the Infection and the Susceptibility through Time, 2022, 19, 1660-4601, 1512, 10.3390/ijerph19031512
    2. Ke Guo, Wanbiao Ma, Global dynamics of a delayed air pollution dynamic model with saturated functional response and backward bifurcation, 2022, 21, 1534-0392, 3831, 10.3934/cpaa.2022124
    3. Meng Li, Ke Guo, Wanbiao Ma, Uniform Persistence and Global Attractivity in a Delayed Virus Dynamic Model with Apoptosis and Both Virus-to-Cell and Cell-to-Cell Infections, 2022, 10, 2227-7390, 975, 10.3390/math10060975
    4. Qing Yang, Xinhong Zhang, Daqing Jiang, Asymptotic behavior of a stochastic SIR model with general incidence rate and nonlinear Lévy jumps, 2022, 107, 0924-090X, 2975, 10.1007/s11071-021-07095-7
    5. Karrar Qahtan Al-Jubouri, Raid Kamel Naji, Arvind Kumar Misra, Dynamics Analysis of a Delayed Crimean‐Congo Hemorrhagic Fever Virus Model in Humans, 2024, 2024, 1110-757X, 10.1155/2024/4818840
    6. Luyao Zhao, Mou Li, Wanbiao Ma, Hopf bifurcation and stability analysis of a delay differential equation model for biodegradation of a class of microcystins, 2024, 9, 2473-6988, 18440, 10.3934/math.2024899
  • Reader Comments
  • © 2016 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(5602) PDF downloads(1176) Cited by(2)

Figures and Tables

Figures(4)  /  Tables(1)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog