
Citation: Robert Russell Monteith Paterson. Depletion of Indonesian oil palm plantations implied from modeling oil palm mortality and Ganoderma boninense rot under future climate[J]. AIMS Environmental Science, 2020, 7(5): 366-379. doi: 10.3934/environsci.2020024
[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 |
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)(r−aS(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=(br−aμ)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)(r−aS(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)(r−aS(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 τi≥0 (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+:‖ϕ1‖≤ra,‖ϕ2‖≤rbaf,‖ϕ3‖≤rd1ae1,‖ϕ4‖≤rbd2afe2} |
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 t≥0,
˙S(t)≤S(t)(r−aS(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 t≥0, we have ˙S(t)≤S(t)(r−aS(t)), which implies that for t≥0,
S(t)≤raϕ1(0)ϕ1(0)+[ra−ϕ1(0)]e−rt≤ra, |
where ϕ1(0)≤ra is used. Further combining the second equation of model (1.3), for t≥0, we have ˙I(t)≤I(t)(rba−fI(t)), which implies that for t≥0,
I(t)≤rbafϕ2(0)ϕ2(0)+[rbaf−ϕ2(0)]e−rbat≤rbaf, |
where ϕ2(0)≤rbaf is used. Thus, for t≥0, we have ˙u1(t)≤rd1a−e1u1(t),˙u2(t)≤rbd2af−e2u2(t). This implies that for t≥0,
u1(t)≤rd1ae1+[ϕ3(0)−rd1ae1]e−e1t≤rd1ae1,u2(t)≤rbd2afe2+[ϕ4(0)−rbd2afe2]e−e2t≤rbd2afe2, |
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=(μ0−d2e2). |
Then the basic reproduction number R0 is defined as the spectral radius of FV−1. Therefore,
R0:=ρ(FV−1)=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
r−aS−bI1+kI−c1u1=0,bS1+kI−μ−fI−c2u2=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(r−bI1+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)(r−bI1+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<rb−rk≡˜I;
(ⅱ) rk≥b, I>0.
Clearly, for both case (ⅰ) and case (ⅱ), we have that
˙F(I)=−be1k(ae1+d1c1)(1+kI)2(r−bI1+kI)−b2e1(ae1+d1c1)(1+kI)3−(f+d2c2e2)<0. |
Hence, F(I) is monotonically decreasing with respect to I and
limI→0+F(I)=bre1ae1+d1c1−μ=μ(R0−1). |
If R0≤1, then limI→0+F(I)≤0 and F(I)=0 has no positive roots. If R0>1, then limI→0+F(I)>0.
For case (ⅰ), we have that
limI→˜I−F(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∗,u∗1,u∗2), where
S∗=e1ae1+d1c1(r−bI∗1+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 τi≥0. 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
|λ−(r−2a¯S−b¯I1+k¯I−c1¯u1)b¯S(1+k¯I)2c1¯Se−λτ10−b¯I1+k¯Iλ−(b¯S(1+k¯I)2−μ−2f¯I−c2¯u2)0c2¯Ie−λτ2−d1e−λτ30λ+e100−d2e−λτ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 τi≥0 (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=μ(R0−1). Obviously, when R0>1, Eq (3.2) has a positive real root λ2>0, and hence, E0 is unstable for any τi≥0 (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 d1c1≤ae1, 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[ω2−ae1S0d1c1S0], |
ω=(−[(aS0)2+e21]+√[(aS0)2+e21]2−4(S0)2(ae1+d1c1)(ae1−d1c1)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 τi≥0(i=1,2,3,4).
Theorem 3.2. For any τ2≥0 and τ4≥0, the following statements are true.
(i) If R0>1, then the disease-free equilibrium E0 is unstable for any τ1≥0 and τ3≥0.
(ii) Assume that d1c1≤ae1. If R0<1, then the disease-free equilibrium E0 is locally asymptotically stable for any τ1≥0 and τ3≥0; If R0=1, then the disease-free equilibrium E0 is linearly stable for any τ1≥0 and τ3≥0.
(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 d1c1≤ae1, 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 d1c1≤ae1. For any τi≥0 (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+a2∫0−τ3(ϕ1(ξ)−S0)2dξ+c1e12d1∫0−τ1(ϕ3(ξ)−u01)2dξ, |
where
V1=ϕ1(0)−S0−S0lnϕ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 t≥0,
dL1dt=dV1dt+a2(S(t)−S0)2−a2(S(t−τ3)−S0)2+c1e12d1(u1(t)−u01)2−c1e12d1(u1(t−τ1)−u01)2, | (3.4) |
where
dV1dt=(S(t)−S0)[a(S0−S(t))−bI(t)1+kI(t)+c1(u01−u1(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)(u01−u1(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,
Λ1≡c1(S(t)−S0)(u01−u1(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
dL1dt≤−a2(1−√d1c1ae1)[(S(t)−S0)2+(S(t−τ3)−S0)2]−[μ−bS01+kI(t)]I(t)−fI2(t)−c2I(t)u2(t−τ2)−c1e12d1(1−√d1c1ae1)[(u1(t)−u01)2+(u1(t−τ1)−u01)2]. | (3.6) |
Note that, if R0≤1, we have that
−[μ−bS01+kI(t)]I(t)=−[μ−bS0+μkI(t)1+kI(t)]I(t)=−[μ(1−R0)1+kI(t)+μkI(t)1+kI(t)]I(t)≤0. | (3.7) |
Assume that d1c1<ae1 and R0≤1. By inequalities (3.6) and (3.7), we can obtain that dL1dt≤0 for t≥0. 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 R0≤1. 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)2−c1e12d1(u1(t−τ1)−u01)2+c1(S(t)−S0)(u01−u1(t−τ1))=−d1c12e1(S(t)−S0+e1d1(u1(t−τ1)−u01))2,−a2(S(t−τ3)−S0)2−c1e12d1(u1(t)−u01)2+c1(u1(t)−u01)(S(t−τ3)−S0)=−d1c12e1(S(t−τ3)−S0−e1d1(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]−[μ(1−R0)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 dL1dt≤0 for t≥0. 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)∈M1⊂G1 for any t∈R. Obviously, I(t)≡0 for any t∈R, 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 t∈R. In addition, according to the first and third equations of model (1.3), we can obtain, for any t∈R,
˙S(t)=S(t)(r−aS(t)−c1u1(t−τ1))=S(t)(r−d1c1e1S(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 t∈R. It is not difficult to find that δ1 and δ2 satisfy
δ1+e1d1δ2=S0+e1d1u01,δ1−e1d1δ2=0, |
which imply that δ1=S0 and δ2=u01. Hence, S(t)≡S0 and u1(t)≡u01 for any t∈R. 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 d1c1≤ae1.
Now, we continue to discuss global dynamics of the disease-free equilibrium E0 in the absence of condition d1c1≤ae1. 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 rb≤aμ (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,rb−aμ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 rb≤aμ, 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 R0≤1. For any τ2≥0 and τ4≥0, the following statements are true.
(i) If rb≤aμ 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+ε0≡S1(ε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 t≥0,
dV1dt=−a(S(t)−S0)2−μ(1−R0)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(S0−S(ξ))−bI(ξ)1+kI(ξ)+c1(u01−u1(ξ−τ1)))dξ. |
In addition, for t>τ1+τ3, we have
Λ2≤c1e12∫tt−τ1((S(t)−S0)2+(u1(ξ)−u01)2)dξ+c1d12∫tt−τ1((S(t)−S0)2+(S(ξ−τ3)−S0)2)dξ=c1(e1+d1)2τ1(S(t)−S0)2+c1e12∫tt−τ1(u1(ξ)−u01)2dξ+c1d12∫tt−τ1(S(ξ−τ3)−S0)2dξ. | (3.13) |
Similarly, for t>t1+τ1+τ3, we have
Λ3≤c1S(ε0)|u1(t)−u01|∫tt−τ3(a|S0−S(ξ)|+bI(ξ)1+kI(ξ)+c1|u01−u1(ξ−τ1)|)dξ≤c1aS(ε0)2∫tt−τ3((u1(t)−u01)2+(S0−S(ξ))2)dξ+c1bS(ε0)2∫tt−τ3((u1(t)−u01)2+I2(ξ)(1+kI(ξ))2)dξ+c21S(ε0)2∫tt−τ3((u1(t)−u01)2+(u01−u1(ξ−τ1))2)dξ=c1S(ε0)2(a+b+c1)τ3(u1(t)−u01)2+c1aS(ε0)2∫tt−τ3(S0−S(ξ))2dξ+c1bS(ε0)2∫tt−τ3I2(ξ)(1+kI(ξ))2dξ+c21S(ε0)2∫tt−τ3(u01−u1(ξ−τ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=6∑i=1Vi, |
where V1 is defined as in the proof of Theorem 3.3,
V2=c1aS(ε0)2∫tt−τ3∫tθ(S(ξ)−S0)2dξdθ,V3=c1d12∫tt−τ1∫tθ(S(ξ−τ3)−S0)2dξdθ+c1d1τ12∫tt−τ3(S(ξ)−S0)2dξ,V4=c1bS(ε0)2∫tt−τ3∫tθI2(ξ)(1+kI(ξ))2dξdθ,V5=c1e12∫tt−τ1∫tθ(u1(ξ)−u01)2dξdθ,V6=c21S(ε0)2∫tt−τ3∫tθ(u1(ξ−τ1)−u01)2dξdθ+c21S(ε0)2τ3∫tt−τ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)2−∫tt−τ3(S(ξ)−S0)2dξ], | (3.15) |
dV3dt=c1d12[τ1(S(t)−S0)2−∫tt−τ1(S(ξ−τ3)−S0)2dξ], | (3.16) |
dV4dt=c1bS(ε0)2[τ3I2(t)(1+kI(t))2−∫tt−τ3I2(ξ)(1+kI(ξ))2dξ], | (3.17) |
dV5dt=c1e12[τ1(u1(t)−u01)2−∫tt−τ1(u1(ξ)−u01)2dξ], | (3.18) |
dV6dt=c21S(ε0)2[τ3(u1(t)−u01)2−∫tt−τ3(u1(ξ−τ1)−u01)2dξ]. | (3.19) |
Combining (3.12)–(3.19), for t>t1+τ1+τ3, we finally have
dL2dt=6∑i=1dVidt≤−[a−c1(e1+2d1)2τ1−c1aS(ε0)2τ3](S(t)−S0)2−μ(1−R0)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)−[c1e1d1−c1e12τ1−c1S(ε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 (ⅰ) rb≤aμ (which implies R0<1) and (H1)–(H2) hold.
From conditions (H1)–(H2), it has that, for sufficiently small ε0>0, the inequalities
a−c1(e1+2d1)2τ1−c1aS(ε0)2τ3>0,e1d1−e12τ1−S(ε0)2(a+b+2c1)τ3>0 | (3.21) |
hold. Furthermore, from Lemma 3.3, we have that limt→∞I(t)=0 and limt→∞u2(t)=0. Thus, there exists a t2>0 such that, for t>t2,
c1bS(ε0)τ3I(t)<μd1c1ae1+d1c1. |
In addition, note that
1−R0=1−bre1μ(ae1+d1c1)=μ(ae1+d1c1)−bre1μ(ae1+d1c1)≥d1c1ae1+d1c1, |
from which we have that, for t>t2,
Λ4:=−μ(1−R0)1+kI(t)I(t)+c1bS(ε0)2τ3I2(t)(1+kI(t))2≤−I(t)1+kI(t)[μd1c1ae1+d1c1−c1bS(ε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 dL2dt≤0 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 L″2(t) is also bounded. This implies that L′2(t) is uniformly continuous for t>˜T. Therefore, it follows from Lemma 3.1 that limt→∞L′2(t)=0. Again from inequalities (3.20)–(3.22), it follows that
limt→∞S(t)=S0,limt→∞u1(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+μk−c1bS(ε0)2τ3]≤0. | (3.24) |
Hence, by inequalities (3.20), (3.21) and (3.24), we can also obtain dL2dt≤0 for t>ˆT:=max{t1+τ1+τ3,t3}. By the same arguments as in Case (i), we can obtain
limt→∞S(t)=S0,limt→∞I(t)=0,limt→∞u1(t)=u01. |
Furthermore, from the last equation of model (1.3), we can easily have limt→∞u2(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 τ2≥0 and τ4≥0, the following statements are true.
(i) If rb≤aμ 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 τ2≥0 and τ4≥0, 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∗,u∗1,u∗2) usually involves some complicated computations, since there are four time delays τi≥0(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 τi≥0(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 τi≥0(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)d1c1≤ae1,τ1≥0,τ3≥0,τ2=τ4=0.(H6)d2c2≤fe2,τ2≥0,τ4≥0,τ1=τ3=0.(H7)d1c1≤ae1,d2c2≤fe2,τ1≥0,τ2≥0,τ3≥0,τ4≥0. |
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,
r−aS∗−bI∗1+kI∗−c1u∗1=0,bS∗1+kI∗−μ−fI∗−c2u∗2=0,−e1u∗1+d1S∗=0,−e2u∗2+d2I∗=0. | (4.1) |
Define a Lyapunov functional W1 on X2 as follows,
W1=(1+kI∗)(ϕ1(0)−S∗−S∗lnϕ1(0)S∗)+ϕ2(0)−I∗−I∗lnϕ2(0)I∗+c1(1+kI∗)2d1(ϕ3(0)−u∗1)2+c22d2(ϕ4(0)−u∗2)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 t≥0,
dW1dt=(1+kI∗)(S(t)−S∗)[a(S∗−S(t))+bI∗1+kI∗−bI(t)1+kI(t)+c1(u∗1−u1(t−τ1))]+(I(t)−I∗)[bS(t)1+kI(t)−bS∗1+kI(t)+bS∗1+kI(t)−bS∗1+kI∗+f(I∗−I(t))+c2(u∗2−u2(t−τ2))]+c1(1+kI∗)d1(u1(t)−u∗1)[−e1(u1(t)−u∗1)+d1(S(t−τ3)−S∗)]+c2d2(u2(t)−u∗2)[−e2(u2(t)−u∗2)+d2(I(t−τ4)−I∗)], |
here Eq (4.1) is used. Note that
(1+kI∗)(S(t)−S∗)(bI∗1+kI∗−bI(t)1+kI(t))+(I(t)−I∗)(bS(t)1+kI(t)−bS∗1+kI(t))=0, |
from which we have that, for t≥0,
dW1dt=−a(1+kI∗)(S(t)−S∗)2−[f+bkS∗(1+kI(t))(1+kI∗)](I(t)−I∗)2−c1e1(1+kI∗)d1(u1(t)−u∗1)2−c2e2d2(u2(t)−u∗2)2+(1+kI∗)Υ1+Π1, | (4.2) |
where
Υ1=c1(S(t)−S∗)(u∗1−u1(t−τ1))+c1(S(t−τ3)−S∗)(u1(t)−u∗1),Π1=c2(I(t)−I∗)(u∗2−u2(t−τ2))+c2(I(t−τ4)−I∗)(u2(t)−u∗2). |
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∗)2∫0−τ3(ϕ1(ξ)−S∗)2dξ+c1e1(1+kI∗)2d1∫0−τ1(ϕ3(ξ)−u∗1)2dξ. |
W2 is continuous on X2, positive definite with respect to E∗, and satisfies condition (ii) of Lemma 3.1 in [35] on ∂X2=¯X2∖X2.
Calculating the derivative of W2 along any solution (S(t),I(t),u1(t),u2(t)) of model (1.3), it follows that, for t≥0,
dW2dt=dW1dt+a(1+kI∗)2[(S(t)−S∗)2−(S(t−τ3)−S∗)2]+c1e1(1+kI∗)2d1[(u1(t)−u∗1)2−(u1(t−τ1)−u∗1)2]. | (4.3) |
By Eqs (4.2) and (4.3), and the following inequality of arithmetic and geometric means:
Υ1≤√d1c1ae1[a2(S(t)−S∗)2+c1e12d1(u1(t−τ1)−u∗1)2+a2(S(t−τ3)−S∗)2+c1e12d1(u1(t)−u∗1)2], | (4.4) |
we have that, for t≥0,
dW2dt≤−a(1+kI∗)2(1−√d1c1ae1)[(S(t)−S∗)2+(S(t−τ3)−S∗)2]−c1e1(1+kI∗)2d1(1−√d1c1ae1)[(u1(t)−u∗1)2+(u1(t−τ1)−u∗1)2]−[f+bkS∗(1+kI(t))(1+kI∗)](I(t)−I∗)2−c2e2d2(u2(t)−u∗2)2. |
It follows from condition (H5) that dW2dt≤0 for t≥0. Hence, the endemic equilibrium E∗ is stable. Furthermore, dW3dt=0 implies I(t)=I∗ and u2(t)=u∗2.
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)=u∗2}. |
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+f2∫0−τ4(ϕ2(ξ)−I∗)2dξ+c2e22d2∫0−τ2(ϕ4(ξ)−u∗2)2dξ. |
W3 is continuous on X2, positive definite with respect to E∗, and satisfies condition (ii) of Lemma 3.1 in [35] on ∂X2=¯X2∖X2.
Calculating the derivative of W3 along any solution (S(t),I(t),u1(t),u2(t)) of model (1.3), it follows that, for t≥0,
dW3dt=dW1dt+f2[(I(t)−I∗)2−(I(t−τ4)−I∗)2]+c2e22d2[(u2(t)−u∗2)2−(u2(t−τ2)−u∗2)2]. | (4.5) |
By Eqs (4.2) and (4.5), and the following inequality of arithmetic and geometric means:
Π1≤√d2c2fe2[f2(I(t)−I∗)2+c2e22d2(u2(t−τ2)−u∗2)2+f2(I(t−τ4)−I∗)2+c2e22d2(u2(t)−u∗2)2], | (4.6) |
we have that, for t≥0,
dW3dt≤−a(1+kI∗)(S(t)−S∗)2−c1e1(1+kI∗)d1(u1(t)−u∗1)2−f2(1−√d2c2fe2)[(I(t)−I∗)2+(I(t−τ4)−I∗)2]−bkS∗(1+kI(t))(1+kI∗)(I(t)−I∗)2−c2e22d2(1−√d2c2fe2)[(u2(t)−u∗2)2+(u2(t−τ2)−u∗2)2]. |
It follows from condition (H6) that dW3dt≤0 for t≥0. Furthermore, dW3dt=0 implies S(t)=S∗ and u1(t)=u∗1. 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∗)2∫0−τ3(ϕ1(ξ)−S∗)2dξ+f2∫0−τ4(ϕ2(ξ)−I∗)2dξ+c1e1(1+kI∗)2d1∫0−τ1(ϕ3(ξ)−u∗1)2dξ+c2e22d2∫0−τ2(ϕ4(ξ)−u∗2)2dξ. |
W4 is continuous on X2, positive definite with respect to E∗, and satisfies condition (ii) of Lemma 3.1 in [35] on ∂X2=¯X2∖X2.
Calculating the derivative of W4 along any solution (S(t),I(t),u1(t),u2(t)) of model (1.3), it follows that, for t≥0,
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)−u∗1)2−(u1(t−τ1)−u∗1)2]+c2e22d2[(u2(t)−u∗2)2−(u2(t−τ2)−u∗2)2]. |
Further, by Eq (4.2) and inequalities (4.4) and (4.6), we have that, for t≥0,
dW4dt≤−a(1+kI∗)2(1−√d1c1ae1)[(S(t)−S∗)2+(S(t−τ3)−S∗)2]−f2(1−√d2c2fe2)[(I(t)−I∗)2+(I(t−τ4)−I∗)2]−bkS∗(1+kI(t))(1+kI∗)(I(t)−I∗)2−c1e1(1+kI∗)2d1(1−√d1c1ae1)[(u1(t)−u∗1)2+(u1(t−τ1)−u∗1)2]−c2e22d2(1−√d2c2fe2)[(u2(t)−u∗2)2+(u2(t−τ2)−u∗2)2]. |
It follows from condition (H7) that dW4dt≤0 for t≥0. Furthermore, if d1c1<ae1, then dW4dt=0 implies that S(t)=S∗ and u1(t)=u∗1. 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 τi≥0(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+bkS∗IM(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+bkS∗IM(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,
dW1dt≤−a(1+kI∗)(S(t)−S∗)2−[f+bkS∗(1+kIM(ε1))(1+kI∗)](I(t)−I∗)2−c1e1(1+kI∗)d1(u1(t)−u∗1)2−c2e2d2(u2(t)−u∗2)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)−u∗1):=Υ2+Υ3,Π1=c2(I(t)−I∗)(u2(t)−u2(t−τ2))+c2(I(t−τ4)−I(t))(u2(t)−u∗2):=Π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(ξ)−u∗1)+d1(S(ξ−τ3)−S∗))dξ|≤c1e12∫tt−τ1((S(t)−S∗)2+(u1(ξ)−u∗1)2)dξ+c1d12∫tt−τ1((S(t)−S∗)2+(S(ξ−τ3)−S∗)2)dξ=c1(e1+d1)2τ1(S(t)−S∗)2+c1e12∫tt−τ1(u1(ξ)−u∗1)2dξ+c1d12∫tt−τ1(S(ξ−τ3)−S∗)2dξ, | (4.8) |
|Υ3|=|−c1(u1(t)−u∗1)∫tt−τ3˙S(ξ)dξ|=|c1(u1(t)−u∗1)∫tt−τ3S(ξ)[a(S∗−S(ξ))+bI∗1+kI∗−bI(ξ)1+kI(ξ)+c1(u∗1−u1(ξ−τ1))]dξ|=|c1(u1(t)−u∗1)∫tt−τ3S(ξ)[a(S∗−S(ξ))+b(I∗−I(ξ))(1+kI∗)(1+kI(ξ))+c1(u∗1−u1(ξ−τ1))]dξ| |
≤c1S1(ε1)|u1(t)−u∗1|∫tt−τ3(a|S∗−S(ξ)|+b1+kI∗|I∗−I(ξ)|+c1|u∗1−u1(ξ−τ1)|)dξ≤c1aS1(ε1)2∫tt−τ3((u1(t)−u∗1)2+(S∗−S(ξ))2)dξ+c1bS1(ε1)2(1+kI∗)∫tt−τ3((u1(t)−u∗1)2+(I∗−I(ξ))2)dξ+c21S1(ε1)2∫tt−τ3((u1(t)−u∗1)2+(u∗1−u1(ξ−τ1))2)dξ=c1S1(ε1)2[a+b1+kI∗+c1]τ3(u1(t)−u∗1)2+c1aS1(ε1)2∫tt−τ3(S∗−S(ξ))2dξ+c1bS1(ε1)2(1+kI∗)∫tt−τ3(I∗−I(ξ))2dξ+c21S1(ε1)2∫tt−τ3(u∗1−u1(ξ−τ1))2dξ, | (4.9) |
|Π2|=|c2(I(t)−I∗)∫tt−τ2˙u2(ξ)dξ|=|c2(I(t)−I∗)∫tt−τ2(−e2(u2(ξ)−u∗2)+d2(I(ξ−τ4)−I∗))dξ|≤c2e22∫tt−τ2((I(t)−I∗)2+(u2(ξ)−u∗2)2)dξ+c2d22∫tt−τ2((I(t)−I∗)2+(I(ξ−τ4)−I∗)2)dξ≤c2(e2+d2)2τ2(I(t)−I∗)2+c2e22∫tt−τ2(u2(ξ)−u∗2)2dξ+c2d22∫tt−τ2(I(ξ−τ4)−I∗)2dξ. | (4.10) |
|Π3|=|−c2(u2(t)−u∗2)∫tt−τ4˙I(ξ)dξ|=|c2(u2(t)−u∗2)∫tt−τ4I(ξ)[bS(ξ)1+kI(ξ)−bS∗1+kI∗+f(I∗−I(ξ))+c2(u∗2−u2(ξ−τ2))]dξ|=|c2(u2(t)−u∗2)∫tt−τ4[bI(ξ)1+kI(ξ)(S(ξ)−S∗)+bkS∗I(ξ)(1+kI∗)(1+kI(ξ))(I∗−I(ξ))+fI(ξ)(I∗−I(ξ))+c2I(ξ)(u∗2−u2(ξ−τ2))]dξ|≤c2|u2(t)−u∗2|∫tt−τ4[bIM(ε1)1+kIM(ε1)|S(ξ)−S∗|+A(ε1)|I∗−I(ξ)|+c2IM(ε1)|u∗2−u2(ξ−τ2)|]dξ≤c2bIM(ε1)2(1+kIM(ε1))∫tt−τ4((u2(t)−u∗2)2+(S(ξ)−S∗)2)dξ+c22A(ε1)∫tt−τ4((u2(t)−u∗2)2+(I(ξ)−I∗)2)dξ+c22IM(ε1)2∫tt−τ4((u2(t)−u∗2)2+(u2(ξ−τ2)−u∗2)2)dξ=c22[bIM(ε1)1+kIM(ε1)+A(ε1)+c2IM(ε1)]τ4(u2(t)−u∗2)2+c2bIM(ε1)2(1+kIM(ε1))∫tt−τ4(S(ξ)−S∗)2dξ+c22A(ε1)∫tt−τ4(I(ξ)−I∗)2dξ+c22IM(ε1)2∫tt−τ4(u2(ξ−τ2)−u∗2)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−τ1∫tθ(S(ξ−τ3)−S∗)2dξdθ+τ1∫tt−τ3(S(ξ)−S∗)2dξ]+c1aS1(ε1)2∫tt−τ3∫tθ(S(ξ)−S∗)2dξdθ+c1bS1(ε1)2(1+kI∗)∫tt−τ3∫tθ(I∗−I(ξ))2dξdθ,+c21S1(ε1)2[∫tt−τ3∫tθ(u1(ξ−τ1)−u∗1)2dξdθ+τ3∫tt−τ1(u1(ξ)−u∗1)2dξ]+c1e12∫tt−τ1∫tθ(u1(ξ)−u∗1)2dξdθ, |
U2=c2d22[∫tt−τ2∫tθ(I(ξ−τ4)−I∗)2dξdθ+τ2∫tt−τ4(I(ξ)−I∗)2dξ]+c22A(ε1)∫tt−τ4∫tθ(I(ξ)−I∗)2dξdθ+c2bIM(ε1)2(1+kIM(ε1))∫tt−τ4∫tθ(S(ξ)−S∗)2dξdθ+c22IM(ε1)2[∫tt−τ4∫tθ(u2(ξ−τ2)−u∗2)2dξdθ+τ4∫tt−τ2(u2(ξ)−u∗2)2dξ]+c2e22∫tt−τ2∫tθ(u2(ξ)−u∗2)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∗)2−c1d12∫tt−τ1(S(ξ−τ3)−S∗)2dξ−c1aS1(ε1)2∫tt−τ3(S(ξ)−S∗)2dξ+c1bS1(ε1)2(1+kI∗)τ3(I(t)−I∗)2−c1bS1(ε1)2(1+kI∗)∫tt−τ3(I(ξ)−I∗)2dξ+[c1e12τ1+c21S1(ε1)2τ3](u1(t)−u∗1)2−c21S1(ε1)2∫tt−τ3(u1(ξ−τ1)−u∗1)2dξ−c1e12∫tt−τ1(u1(ξ)−u∗1)2dξ, | (4.12) |
dU2dt=c2bIM(ε1)2(1+kIM(ε1))τ4(S(t)−S∗)2−c2bIM(ε1)2(1+kIM(ε1))∫tt−τ4(S(ξ)−S∗)2dξ+[c2d22τ2+c22A(ε1)τ4](I(t)−I∗)2−c2d22∫tt−τ2(I(ξ−τ4)−I∗)2dξ−c22A(ε1)∫tt−τ4(I(ξ)−I∗)2dξ+[c2e22τ2+c22IM(ε1)2τ4](u2(t)−u∗2)2−c22IM(ε1)2∫tt−τ4(u2(ξ−τ2)−u∗2)2dξ−c2e22∫tt−τ2(u2(ξ)−u∗2)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τ1−c1aS1(ε1)2τ3−c2bIM(ε1)2(1+kI∗)(1+kIM(ε1))τ4](S(t)−S∗)2−[f+bkS∗(1+kIM(ε1))(1+kI∗)−(2d2+e2)c22τ2−c1bS1(ε1)2τ3−c22A(ε1)τ4](I(t)−I∗)2−(1+kI∗)c1[e1d1−e12τ1−S1(ε1)2(a+b1+kI∗+2c1)τ3](u1(t)−u∗1)2−c2[e2d2−e22τ2−12(bIM(ε1)1+kIM(ε1)+A(ε1)+2c2IM(ε1))τ4](u2(t)−u∗2)2:=−(1+kI∗)Q1(ε1)(S(t)−S∗)2−Q2(ε1)(I(t)−I∗)2−(1+kI∗)c1Q3(ε1)(u1(t)−u∗1)2−c2Q4(ε1)(u2(t)−u∗2)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 dUdt≤0 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)=u∗1,limt→+∞u2(t)=u∗2. |
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 d1c1≤ae1, 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, R0≈0.738<1, E0=(0.462,0,0.923,0) and τ013≈4.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<τ013≈4.542 and unstable for τ1+τ3>τ013≈4.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<τ013≈4.542 and (ˆH) are satisfied. Figure 1(a) shows that the disease-free equilibrium E0 is globally asymptotically stable.
Let us choose τ1=2 and τ3=1.25. We see that condition (ˆH) does not hold, but condition τ1+τ3<τ013≈4.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>τ013≈4.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 R0≈3.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.
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).
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].
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] | Paterson RRM, Lima N (2018) Climate change affecting oil palm agronomy, and oil palm cultivation increasing climate change, require amelioration. Ecol Evol 8: 452-461. |
[2] | Paterson RRM (2019) Ganoderma boninense disease of oil palm is expected to significantly reduce production after 2050 in Sumatra if projected climate change occurs. Microorganisms 7: 24. |
[3] | Paterson RRM (2019) Ganoderma boninense disease deduced from simulation modelling with large data sets of future Malaysian oil palm climate. Phytoparasitica 47: 255-262. |
[4] | Paterson RRM (2020). Oil palm survival under climate change in Kalimantan and alternative SE Asian palm oil countries with future basal stem rot assessments. Forest Pathol 50: e12604. |
[5] | Paterson RRM (2020) Future scenarios for oil palm mortality and infection by Phytophthora palmivora in Colombia, Ecuador and Brazil, extrapolated to Malaysia and Indonesia. Phytoparasitica 48: 513-523. |
[6] | Sarkar SK, Begum RA, Pereira JJ (2020) Impacts of climate change on oil palm production in Malaysia. Environ Sci Pollut R 27: 9760-9770. |
[7] | Corley RHV, Tinker PB (2015) The Oil Palm, Wiley Blackwell. |
[8] | Lam YW, Kulak M, Sim S, et al. (2019) Greenhouse gas footprints of palm oil production in Indonesia over space and time. Sci Total Environ 688: 827-837. |
[9] | Rianto B (2010) Overview of palm oil industry in Indonesia. Pricewaterhouse Coopers Indonesia. Available from: https://www.pwc.com/id/en/publications/assets/palm-oil-plantation.pdf |
[10] | Paterson RRM, Kumar L, Taylor S, et al. (2015) Future climate effects on suitability for growth of oil palms in Malaysia and Indonesia. Sci Report 5. |
[11] | Paterson RRM, Kumar L, Shabani F, et al. (2017) World climate suitability projections to 2050 and 2100 for growing oil palm. J Agric Sci 155: 689-702. |
[12] | Dislich C, Keyel AC, Salecker J, et al. (2017) Wiegand, K. A review of the ecosystem functions in oil palm plantations, using forests as a reference system. Biol Rev 49: 1539-1569. |
[13] | Meijaard E, Garcia-Ulloa J, Sheil D, et al. (2018) Oil palm and biodiversity: A situation analysis by the IUCN Oil Palm Task Force. 2018. |
[14] | Kadandale S, Marten R, Smith R (2019) The palm oil industry and noncommunicable diseases. Bull World Healt Organ 97: 118-128. |
[15] | Gibb R, Redding DW, Chin KQ, et al. (2020). Zoonotic host diversity increases in human-dominated ecosystems. Nature 584: 398-402. |
[16] | Paterson RRM, Sariah M, Lima N (2013) How will climate change affect oil palm fungal diseases? Crop Prot 46: 113-120. |
[17] | Zhou LW, Cao Y, Wu SH, et al. (2015) Global diversity of the Ganoderma lucidum complex (Ganodermataceae, Polyporales) inferred from morphology and multilocus phylogeny. Phytochem 114: 7-15. |
[18] | Merciere M, Boulord R, Carasco-Lacombe C, et al. (2017) About Ganoderma boninense in oil palm plantations of Sumatra and peninsular Malaysia: Ancient population expansion, extensive gene flow and large scale dispersion ability. Fungal Biol 121: 529-540. |
[19] | Paterson R R M, Lima N. Ecology and biotechnology of thermophilic fungi on crops under global warming[M]//Fungi in Extreme Environments: Ecological Role and Biotechnological Significance. Springer, Cham. 2019: 81-96. |
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 |