# Stability and Bifurcation Analysis on an Ecoepidemiological  Model with Stage Structure and Time Delay

Authors: Lingshu Wang; Guanghui Feng
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

An ecoepidemiological predator-prey model with stage structure for the predator and time delay due to the gestation of the predator is investigated. The effects of a prey refuge with disease in the prey population are concerned. By analyzing the corresponding characteristic equations, the local stability of each of the feasible equilibria of the model is discussed. Further, it is proved that the model undergoes a Hopf bifurcation at the positive equilibrium. By means of appropriate Lyapunov functions and LaSalle’s invariance principle, sufficient conditions are obtained for the global stability of the semitrivial boundary equilibria. By using an iteration technique, sufficient conditions are derived for the global attractiveness of the positive equilibrium.


---

## Body

## 1. Introduction

In the natural world, species does not exist alone. While species spreads the disease, it also competes with the other species for space or food, or it is predated by other species. The construction and study of models for the population dynamics of predator-prey systems have been an important topic in theoretical ecology. Following Anderso and May [1], who were the first to propose an ecoepidemiological model by merging the ecological predator-prey model introduced by Lotka and Volterra, the effect of disease in ecological system is an important issue from mathematical and ecological point of view. Ecoepidemiology which is a relatively new branch of study in theoretical biology tackles such situations by dealing with both ecological and epidemiological issues.The research of the hiding behaviour of preys has been incorporated as a new ingredient of predator-prey models. In nature, prey populations often have access to areas where they are safe from their predators. Such refugia are usually playing two significant roles, serving both to reduce the chance of extinction due to predation and to damp predator-prey oscillations. It is well known that many more attentions have been paid on the effects of a prey refuge for predator-prey model. In [2], Wang considered an ecoepidemiological model incorporating a prey refuge with disease in the prey population
(1)dSdt=rS(t)(1-S(t)+I(t)K)-βS(t)I(t)-b1(1-m)S(t)Y(t),dIdt=βS(t)I(t)-b2(1-m)I(t)Y(t)-dI(t),dYdt=pb1(1-m)S(t)Y(t)+pb2(1-m)I(t)Y(t)-cY(t),
where S(t) and I(t) represent the densities of susceptible and infected prey population at time t, respectively, and Y(t) represents the density of the predator population at time t. The parameters r, K, β, b1, b2, c, d, and p are positive constants in which r and K represent the prey intrinsic growth rate and the carrying capacity, respectively. β is the transmission rate of the susceptible prey into the infected prey. b1 and b2 are the capturing rates of the susceptible prey and the infected prey, respectively. p describes the efficiency of the predator in converting consumed prey into predator offspring. The constant proportion infected prey refuge is (1-m)I, where m∈[0,1) is a constant. By means of appropriate Lyapunov functions and limit theory, sufficient conditions are obtained for the global stability of the semitrivial boundary equilibria of model (1).We note that it is assumed in system (1) that each individual predator admits the same ability to feed on prey. This assumption seems to be not realistic for many animals. In the natural world, there are many species whose individuals pass through an immature stage during which they are raised by their parents, and the rate at which they attack prey can be ignored. Moreover, it can be assumed that their reproductive rate during this stage is zero. Stage structure is a natural phenomenon and represents, for example, the division of a population into immature and mature individuals. Stage-structured models have received great attention in recent years (see, e.g., [3–5]).Time delays of one type or another have been incorporated into biological models by many researchers (see, e.g., [5–7]). In general, delay differential equations exhibit much more complicated dynamics than ordinary differential equations since a time delay could cause the population to fluctuate. Time delay due to gestation is a common example, because generally the consumption of prey by the predator throughout its past history governs the present birth rate of the predator. Therefore, more realistic models of population interactions should take into account the effect of time delays.Based on the above discussions, in this paper, we incorporate a stage structure for the predator and time delay due to the gestation of predator into the model (1). To this end, we study the following differential equations:
(2)dSdt=rS(t)(1-S(t)+I(t)K)-βS(t)I(t),dIdt=βS(t)I(t)-dI(t)-b(1-m)I(t)Y2(t),dY1dt=pb(1-m)I(t-τ)Y2(t-τ)-(r1+d1)Y1(t),dY2dt=r1Y1(t)-d2Y2(t)-aY22(t),
where Y1(t) and Y2(t) represent the densities of the immature and the mature predator population at time t, respectively. The parameters d1, d2, and r1 are positive constants in which d1 and d2 are the death rates of the immature and the mature predator, respectively. r1 denotes the rate of immature predator becoming mature predator. τ≥0 is a constant delay due to the gestation of the predator.The initial conditions for system (2) take the form
(3)S(θ)=ϕ1(θ)≥0,I(θ)=ϕ2(θ)≥0,Y1(θ)=φ1(θ)≥0,Y2(θ)=φ2(θ)≥0,θ∈[-τ,0),ϕ1(0)>0,ϕ2(0)>0,φ1(0)>0,φ2(0)>0,(ϕ1(θ),ϕ2(θ),φ1(θ),φ2(θ))∈C([-τ,0],R+04),
where R+04={(x1,x2,x3,x4):xi≥0,i=1,2,3,4}.It is well known by the fundamental theory of functional differential equations [8] that model (2) has a unique solution (S(t),I(t),Y1(t),Y2(t)) satisfying initial conditions (3).The organization of this paper is as follows. In the next section, we show the positivity and the boundedness of solutions of model (2) with initial conditions (3). In Section 3, we investigate the stability of the semitrivial equilibria of the model (2). In Section 4, we discuss the stability of the positive equilibrium of the model (2). Further, we study the existence of Hopf bifurcation at the positive equilibrium. A brief discussion is given in Section 5 to conclude this work.


## 2. Preliminaries

In this section, we show the positivity and the boundedness of solutions of model (2) with initial conditions (3).Theorem 1.
Solutions of model (2) with initial conditions (3) are positive for all t≥0.Proof.
Let(S(t),I(t),Y1(t),Y2(t)) be a solution of model (2) with initial conditions (3). It follows from the first and the second equations of model (2) that
(4)S(t)=S(0)exp⁡{∫0t‍[r-rKS(s)-(rK+β)I(s)]ds}>0,I(t)=I(0)exp⁡{∫0t‍[βS(s)-d-b(1-m)Y2(s)]ds}>0.
Let us considerY1(t) and Y2(t) for t∈[0,τ]. Since ϕ2(θ)≥0,φ2(θ)≥0, for θ∈[-τ,0], we derive from the third equation of model (2) that
(5)dY1dt≥-(r1+d1)Y1(t).
Since φ1(0)>0, a standard comparison argument shows that
(6)Y1(t)≥Y1(0)e-(r1+d1)t>0;
that is, Y1(t)>0, for t∈[0,τ]. For t∈[0,τ], it follows from the fourth equation of (2) that
(7)dY2dt≥-d2Y2(t)-aY22(t).
Since φ2(0)>0, a standard comparison argument shows that
(8)Y2(t)≥Y2(0)exp⁡{∫0t‍(-d2-aY2(s))ds}>0;
that is, Y2(t)>0, for t∈[0,τ]. In a similar way, we treat the intervals [τ,2τ],…,[nτ,(n+1)τ],n∈N. Thus, S(t)>0,I(t)>0,Y1(t)>0, and Y2(t)>0 for all t≥0. This completes the proof.Theorem 2.
Positive solutions of model (2) with initial conditions (3) are ultimately bounded.Proof.
Let(S(t),I(t),Y1(t),Y2(t)) be a positive solution of model (2) with initial conditions (3). Denote d^=min⁡{d,d1,d2}. Define
(9)W(t)=pS(t-τ)+pI(t-τ)+Y1(t)+Y2(t).
Calculating the derivative of W(t) along the positive solutions of (2), it follows that
(10)dWdt=prS(t-τ)-prKS2(t-τ)-prKS(t-τ)I(t-τ)-pdI(t-τ)-d1Y1(t)-d2Y2(t)-aY22(t)≤-d^W(t)+p(r+d^)S(t-τ)-prKS2(t-τ)-prKS(t-τ)I(t-τ)≤-d^W(t)-prK[S(t-τ)-K(r+d^)2r]2+pK(r+d^)24r≤-d^W(t)+pK(r+d^)24r,
which yields
(11)lim⁡sup⁡t→∞W(t)≤pK(r+d^)24rd^.
If we choose M1=K(r+d^)2/4rd^, M2=pK(r+d^)2/4rd^, then
(12)lim⁡sup⁡t→∞S(t)≤M1,lim⁡sup⁡t→∞I(t)≤M1,lim⁡sup⁡t→∞Yi(t)≤M2,(i=1,2).
This completes the proof.


## 3. Boundary Equilibria and Their Stability

In this section, we discuss the stability of the boundary equilibria of model (2).Model (2) always has two boundary equilibria, namely, the trivial equilibrium E0(0,0,0,0) and the axial equilibrium EK(K,0,0,0). It is easy to show that if Kβ>d, model (2) admits a predator-extinction equilibrium E1(S1,I1,0,0), where
(13)S1=dβ,I1=r(Kβ-d)β(Kβ+r).The characteristic equation of model (2) at the equilibrium E0(0,0,0,0) is of the form
(14)(λ-r)(λ+d)(λ+r1+d1)(λ+d2)=0.
Clearly (14) has a positive real root. Accordingly, the equilibrium E0 is unstable.The characteristic equation of model (2) at the equilibrium EK(K,0,0,0) takes the form
(15)(λ+r)(λ+d2)(λ+r1+d1)[λ-(Kβ-d)]=0.
Hence, if Kβ<d, (15) has no positive real root. Accordingly, the equilibrium EK is locally asymptotically stable. If Kβ>d, (15) has a positive real root. Accordingly, the equilibrium EK is unstable.Theorem 3.
IfKβ<d, then the semitrivial equilibrium EK is globally stable.Proof.
Based on the above discussions, we only prove the global attractivity of the equilibriumEK. Let
(16)VK(t)=c1[S(t)-K-Kln⁡S(t)K]+I(t)+1pY1(t)+1pY2(t)+b(1-m)∫t-τt‍I(s)Y2(s)ds,
where c1=Kβ/(Kβ+r). Calculating the derivative of VK(t) along the positive solutions of model (2), it follows that
(17)ddtVK(t)=c1S(t)-KS(t)S˙(t)+I˙(t)+1pY˙1(t)+1pY˙2(t)+b(1-m)I(t)Y2(t)-b(1-m)I(t-τ)Y2(t-τ).=-βrKβ+r[S(t)-K]2-(d-Kβ)I(t)-1pd1Y1(t)-1pd2Y2(t)-apY22(t).
If Kβ<d, then it follows from (17) that V˙K(t)≤0. By Theorem  5.3.1, in [8], solutions are limited to M, the largest invariant subset of {V˙K(t)=0}. Clearly, we see from (17) that V˙K(t)=0 if and only if S(t)=K,I(t)=0,Y1(t)=0,Y2(t)=0. Accordingly, the global asymptotic stability of EK follows from LaSalle’s invariant principle. This completes the proof.The characteristic equation of model (2) at the equilibrium E1 is of the form
(18)(λ2+rKS1λ+β(Kβ+r)KS1I1)×(λ2+g1λ+g0+f0e-λτ)=0,
where g1=r1+d1+d2,g0=d2(r1+d1),f0=-pbr1(1-m)I1. Clearly, the roots of equation λ2+(r/K)S1λ+(β(Kβ+r)/K)S1I1=0 have negative real part. When τ=0, if pbrr1(1-m)(Kβ-d)<βd2(r1+d1)(Kβ+r), then the roots of (18) have negative real part. Accordingly, E1 is locally asymptotically stable. If pbrr1(1-m)(Kβ-d)>βd2(r1+d1)(Kβ+r), then E1 is unstable. It is easily seen that
(19)g12-2g0=(r1+d1)2+d22>0,g02-f02>0.Hence, if0<pbrr1(1-m)(Kβ-d)<βd2(r1+d1)(Kβ+r), by Lemma B in [7], it follows that the equilibrium E1 is locally asymptotically stable for all τ≥0. If pbrr1(1-m)(Kβ-d)>βd2(r1+d1)(Kβ+r), then E1 is unstable for all τ≥0.Theorem 4.
LetKβ>d hold; the predator-extinction equilibrium E1 of model (2) is globally stable provided that
(20)pbrr1(1-m)(Kβ-d)<βd2(r1+d1)(Kβ+r).Proof.
Based on the above discussions, we only prove the global attractivity of the equilibriumE1. Define
(21)V11(t)=c1(S(t)-S1-S1ln⁡S(t)S1)+I(t)-I1-I1ln⁡I(t)I1+k1Y1(t)+k2Y2(t),
where c1=Kβ/(Kβ+r) and k1=1/p,k2=(r1+d1)/(pr1). Calculating the derivative of V11(t) along the positive solutions of (2), it follows that
(22)V˙11(t)=c1S(t)-S1S(t)S˙(t)+I(t)-I1I(t)I˙(t)+k1Y˙1(t)+k2Y˙2(t)=-βrKβ+r[S(t)-S1]2-b(1-m)I(t)Y2(t)+b(1-m)I(t-τ)Y2(t-τ)-[d2(r1+d1)pr1-b(1-m)I1]Y2(t)-a(r1+d1)pr1Y22(t).
Define(23)V1(t)=V11(t)+b(1-m)∫t-τt‍I(s)Y2(s)ds.
We derive from (22) and (23) that
(24)V˙1(t)=-βrKβ+r[S(t)-S1]2-[d2(r1+d1)pr1-b(1-m)I1]Y2(t)-a(r1+d1)pr1Y22(t).
If 0<pbrr1(1-m)(Kβ-d)<βd2(r1+d1)(Kβ+r), it then follows from (24) that V˙1(t)≤0. By Theorem  5.3.1, in [8], solutions are limited to M, the largest invariant subset of {V˙1(t)=0}. Clearly, we see from (24) that V˙1(t)=0, if and only if S(t)=S1,Y2(t)=0. It follows from the first and fourth equations of (2) that 0=S˙(t)=r-(r/K)S1-((Kβ+r)/K)I(t),0=Y˙2(t)=r1Y1(t), which yields I(t)=I1,Y1(t)=0. Using LaSalle’s invariant principle, the global asymptotic stability of E1 follows. This completes the proof.


## 4. Stability of Positive Equilibrium

In this section, we are concerned with the stability of the positive equilibriumE* and the existence of Hopf bifurcations at the positive equilibrium E* of model (2).If the following holds,(H1)
pbrr1(1-m)(Kβ-d)>βd2(r1+d1)(Kβ+r), then model (2) has a unique positive equilibrium E*(S*,I*,Y1*,Y2*), where
(25)S*=(Kprr1b2(1-m)2+ad(r1+d1)(Kβ+r)-b(1-m)d2(r1+d1)(Kβ+r)Kprr1b2)×(prr1b2(1-m)2+βa(r1+d1)(Kβ+r))-1,I*=r(K-S*)Kβ+r,Y1*=d2+aY2*r1Y2*,Y2*=βS*-db(1-m).The characteristic equation of model (2) at the equilibrium E* takes the form
(26)λ4+p3λ3+p2λ2+p1λ+p0+(q2λ2+q1λ+q0)e-λτ=0,
where
(27)p3=r1+d1+d2+2aY2*+rKS*,p2=(r1+d1)(d2+2aY2*)+rKS*(r1+d1+d2+2aY2*)+Kβ+rKβS*I*,p1=Kβ+rKβS*I*(r1+d1+d2+2aY2*)+rKS*(r1+d1)(d2+2aY2*),p0=Kβ+rKβS*I*(r1+d1)(d2+2aY2*),q2=-pbr1(1-m)I*,q1=pbr1(1-m)I*[b(1-m)Y2*-rKS*],q0=pbr1(1-m)S*I*[b(1-m)rKY2*-Kβ+rKβI*].
It is easy to show that
(28)p3>0,p0+q0>0,p1+q1>0,p2+q2>0.Whenτ=0, (26) becomes
(29)λ4+p3λ3+(p2+q2)λ2+(p1+q1)λ+p0+q0=0.
If the following holds, (H2)
(p1+q1)[p3(p2+q2)-(p0+q0)]>p32(p0+q0), then by the Routh-Hurwitz theorem, when τ=0, the coexistence equilibrium E* of model (2) is locally asymptotically stable and E* is unstable if (p1+q1)[p3(p2+q2)-(p0+q0)]<p32(p0+q0).Ifiω(ω>0) is a solution of (26), separating real and imaginary parts, we have
(30)(q2ω2-q0)sin⁡ωτ+q1ωcos⁡⁡ωτ=p3ω3-p1ω,(q2ω2-q0)cos⁡⁡ωτ-q1ωsin⁡ωτ=ω4-p2ω2+p0.
Squaring and adding the two equations of (30), it follows that
(31)ω8+h3ω6+h2ω4+h1ω2+h0=0,
where
(32)h3=p32-2p2,h2=p22+2p0-2p1p3-q22,h1=p12-2p0p2+2q0q2-q12,h0=p02-q02.Assume that the following holds:(H3)
h3>0,h2>0,h1>0.If h0>0, by the general theory on characteristic equations of delay differential equations from [9] (Theorem 4.1), E* remains stable for all τ>0. If h0<0, then (31) has a unique positive root ω0; that is, (26) admits a pair of purely imaginary roots of the form ±iω0. From (30), we see that
(33)τn=2nπω0+1ω0arccos((q2ω02-q0)(ω04-p2ω02+p0)+q1ω0(p3ω03-p1ω0))×((q1ω0)2+(q2ω02-q0)2)-1,n=0,1,2,….
By Theorem  3.4.1, in [9], we see that E* remains stable for τ<τ0.In the following, we claim that(34)d(Re(λ))dτ|τ=τ0>0.
This will show that there exists at least one eigenvalue with a positive real part for τ>τ0. Moreover, the conditions for the existence of a Hopf bifurcation (Theorem  2.9.1 in [9]) are then satisfied yielding a periodic solution. To this end, differentiating equation (26) with respect to τ, it follows that
(35)(dλdτ)-1=4λ3+3p3λ2+2p2λ+p1-λ(λ4+p3λ3+p2λ2+p1λ+p0)+2q2λ+q1λ(q2λ2+q1λ+q0)-τλ.
Hence, a direct calculation shows that
(36)sgn⁡{d(Reλ)dτ}λ=iω0=sgn⁡{Re(dλdτ)-1}λ=iω0=sgn⁡{×((q1ω0)2+(q2ω02-q0)2)-1((3p3ω02-p1)(p3ω02-p1)+2(2ω02-p2)×(ω04-p2ω02+p0))×(ω02(p1-p3ω02)2+(ω04-p2ω02+p0)2)-1+(-q12+2q2(q0-q2ω02))×((q1ω0)2+(q2ω02-q0)2)-1}.
We derive from (30) that
(37)ω02(p1-p3ω02)2+(ω04-p2ω02+p0)2=(q1ω0)2+(q2ω02-q0)2.
Hence, it follows that
(38)sgn⁡{d(Reλ)dτ}λ=iω0=sgn⁡{4ω06+3h3ω04+2h2ω02+h1(q1ω0)2+(q2ω02-q0)2}>0.
Therefore, if (H3) holds, then the transversal condition holds and a Hopf bifurcation occurs at ω=ω0,τ=τ0.In conclusion, we have the following results.Theorem 5.
For model (2), let (H1) hold, and we have the following.(i)
If (H2) and (H3) hold, h0>0, then the positive equilibrium E* is locally asymptotically stable for all τ≥0.(ii)
If (H2) and (H3) hold, h0<0, then there exists a positive number τ0, such that the positive equilibrium E* is locally asymptotically stable if 0≤τ<τ0 and is unstable if τ>τ0. Further, model (2) undergoes a Hopf bifurcation at E* when τ=τ0.(iii)
If(p1+q1)[p3(p2+q2)-(p1+q1)]<p32(p0+q0), then the positive equilibrium E* is unstable for all τ≥0.Now, we are concerned with the global attractiveness of the positive equilibriumE*.Theorem 6.
Let (H1) hold, and then the positive equilibrium E*(S*,I*,Y1*,Y2*) of model (2) is globally attractive provided that
(39)aβ(r1+d1)(Kβ+r)≠prr1b2(1-m)2.Proof.
Let(S(t),I(t),Y1(t),Y2(t)) be any positive solution of model (2) with initial conditions (3). Let
(40)MS=lim⁡sup⁡t→+∞S(t),mS=lim⁡inf⁡t→+∞S(t),MI=lim⁡sup⁡t→+∞I(t),mI=lim⁡inf⁡t→+∞I(t),MYi=lim⁡sup⁡t→+∞Yi(t),mYi=lim⁡inf⁡t→+∞Yi(t),(i=1,2).
We now claim that MS=mS=S*,MI=mI=I*,MYi=mYi=Yi*(i=1,2). The technique of proof is to use an iteration method.
We derive from the first and the second equations of model (2) that
(41)dSdt=rS(1-S+IK)-βSI,dIdt≤βSI-dI.
Consider the following auxiliary equations:
(42)dx1dt=rx1(1-x1+x2K)-βx1x2,dx2dt=βx1x2-dx2.
If Kβ>d, then, by Theorem  3.1 in [2], it follows from (42) that
(43)lim⁡t→+∞x1(t)=dβ,lim⁡t→+∞x2(t)=r(Kβ-d)β(Kβ+r).
By comparison, we obtain that
(44)MS=lim⁡sup⁡t→+∞S(t)≤dβ≔M1S,MI=lim⁡sup⁡t→+∞I(t)≤r(Kβ-d)β(Kβ+r):=M1I.
Hence, for ɛ>0, sufficiently small, there is a T1>0 such that if t>T1, then I(t)≤M1I+ɛ. We therefore derive from the third and the fourth equations of model (2) that, for t>T1+τ,
(45)dY1dt≤pb(1-m)(M1I+ɛ)Y2(t-τ)-(r1+d1)Y1(t),dY2dt=r1Y1(t)-d2Y2(t)-aY22(t).
Consider the following auxiliary equations:
(46)dz1dt=pb(1-m)(M1I+ɛ)z2(t-τ)-(r1+d1)z1(t),dz2dt=r1z1(t)-d2z2(t)-az22(t).
If (H1) holds, then, by Lemma  2.4 in [10], it follows from (46) that
(47)lim⁡t→+∞z1(t)=(pb(1-m)(M1I+ɛ)×[pbr1(1-m)(M1I+ɛ)-d2(r1+d1)])×(a(r1+d1)2)-1,lim⁡t→+∞z2(t)=pbr1(1-m)(M1I+ɛ)-d2(r1+d1)a(r1+d1).
By comparison, for ɛ>0, sufficiently small, we obtain that
(48)MY1=limsup⁡t→+∞Y1(t)≤pb(1-m)M1I[pbr1(1-m)M1I-d2(r1+d1)]a(r1+d1)2:=M1Y1,MY2=limsup⁡t→+∞Y2(t)=pbr1(1-m)M1I-d2(r1+d1)a(r1+d1):=M1Y2.
Hence, for ɛ>0, sufficiently small, there is a T2≥T1+τ such that if t>T2, then Y2(t)≤M1Y2+ɛ.
Forɛ>0, sufficiently small, we derive from the first and the second equations of model (2) that, for t>T2,
(49)dSdt=rS(1-S+IK)-βSI,dIdt≥βSI-dI-b(1-m)(M1Y2+ɛ)I.
Consider the following auxiliary equations:
(50)dx1dt=rx1(1-x1+x2K)-βx1x2,dx2dt=βx1x2-dx2-b(1-m)(M1Y2+ɛ)x2.
If (H1) holds, then, by Theorem  3.1 in [2], it follows from (50) that
(51)lim⁡t→+∞x1(t)=d+b(1-m)(M1Y2+ɛ)β,lim⁡t→+∞x2(t)=rKβ+r[K-d+b(1-m)(M1Y2+ɛ)β].
By comparison, for ɛ>0, sufficiently small, we conclude that
(52)mS=lim⁡inf⁡t→+∞S(t)≥d+b(1-m)M1Y2β≔N1S,mI=lim⁡inf⁡t→+∞≥rKβ+r(K-d+b(1-m)M1Y2β):=N1I.
Hence, for ɛ>0, sufficiently small, there is a T3≥T2 such that if t>T3, then I(t)≥N1I-ɛ. For ɛ>0, sufficiently small, we derive from the third and the fourth equations of model (2) that for t>T3+τ(53)dY1dt≥pb(1-m)(N1I-ɛ)Y2(t-τ)-(d1+r1)Y1(t),dY2dt=r1Y1(t)-d2Y2(t)-aY22(t).
Consider the following auxiliary equations:
(54)dz1dt=pb(1-m)(N1I-ɛ)z2(t-τ)-(r1+d1)z1(t),dz2dt=r1z1(t)-d2z2(t)-az22(t).
Since (H1) holds, by Lemma 2.4 of [10], it follows from (54) that
(55)lim⁡t→+∞z1(t)=(pb(1-m)(N1I-ɛ)×[pbr1(1-m)(N1I-ɛ)-d2(r1+d1)])×(a(r1+d1)2)-1,lim⁡t→+∞z2(t)=pbr1(1-m)(N1I-ɛ)-d2(r1+d1)a(r1+d1).
By comparison, for ɛ>0, sufficiently small, we obtain that
(56)mY1=lim⁡inf⁡t→+∞Y1(t)≥pb(1-m)N1I[pbr1(1-m)N1I-d2(r1+d1)]a(r1+d1)2:=N1Y1,mY2=liminf⁡t→+∞Y2(t)≥pbr1(1-m)N1I-d2(r1+d1)a(r1+d1):=N1Y2.
Hence, for ɛ>0, sufficiently small, there is a T4≥T3+τ, such that if t>T4, Y2(t)≥N1Y2-ɛ.
Forɛ>0, sufficiently small, we derive from the first and the second equations of model (2) that, for t>T4,
(57)dSdt=rS(1-S+IK)-βSI,dIdt≤βSI-dI-b(1-m)(N1Y2-ɛ)I.
Consider the following auxiliary equations:
(58)dx1dt=rx1(1-x1+x2K)-βx1x2,dx2dt=βx1x2-dx2-b(1-m)(N1Y2-ɛ)x2.
If (H1) holds, then, by Theorem  3.1 in [2], it follows from (58) that
(59)lim⁡t→+∞x1(t)=d+b(1-m)(N1Y2-ɛ)β,lim⁡t→+∞x2(t)=rKβ+r[K-d+b(1-m)(N1Y2-ɛ)β].
By comparison, for ɛ>0, sufficiently small, we obtain that
(60)MS=lim⁡sup⁡t→+∞S(t)≤d+b(1-m)N1Y2β≔M2S,MI=lim⁡sup⁡t→+∞≤rKβ+r(K-d+b(1-m)N1Y2β):=M2I.
Therefore, for ɛ>0, sufficiently small, there is a T5≥T4 such that if t>T5, I(t)≤M2I+ɛ.
Forɛ>0, sufficiently small, we derive from the third and the fourth equations of model (2) that, for t>T5+τ,
(61)dY1dt≤pb(1-m)(M2I+ɛ)Y2(t-τ)-(d1+r1)Y1(t),dY2dt=r1Y1(t)-d2Y2(t)-aY22(t).
Consider the following auxiliary equations:
(62)dz1dt=pb(1-m)(M2I+ɛ)z2(t-τ)-(r1+d1)z1(t),dz2dt=r1z1(t)-d2z2(t)-az22(t).
Since (H1) holds, by Lemma 2.4 of [10], it follows from (62) that
(63)lim⁡t→+∞z1(t)=(pb(1-m)(M2I+ɛ)×[pbr1(1-m)(M2I+ɛ)-d2(r1+d1)])×(a(r1+d1)2)-1,lim⁡t→+∞z2(t)=pbr1(1-m)(M2I+ɛ)-d2(r1+d1)a(r1+d1).
By comparison, for ɛ>0, sufficiently small, we conclude that
(64)MY1=limsup⁡t→+∞Y1(t)≥pb(1-m)M2I[pbr1(1-m)M2I-d2(r1+d1)]a(r1+d1)2:=M2Y1,MY2=limsup⁡t→+∞Y2(t)≥pbr1(1-m)M2I-d2(r1+d1)a(r1+d1):=M2Y2.
Therefore, for ɛ>0, sufficiently small, there is a T6≥T5+τ such that if t>T6, y2(t)≤M2Y2+ɛ.
Forɛ>0, sufficiently small, it follows from the first and the second equations of model (2) that for t>T6(65)dSdt=rS(1-S+IK)-βSI,dIdt≥βSI-dI-b(1-m)(M2Y2+ɛ)I(t).
Consider the following auxiliary equations:
(66)dx1dt=rx1(1-x1+x2K)-βx1x2,dx2dt=βx1x2-dx2-b(1-m)(M2Y2+ɛ)x2.
If (H1) holds, then, by Theorem  3.1 in [2], it follows from (66) that
(67)lim⁡t→+∞x1(t)=d+b(1-m)(M2Y2+ɛ)β,lim⁡t→+∞x2(t)=rKβ+r[K-d+b(1-m)(M2Y2+ɛ)β].
By comparison, for ɛ>0, sufficiently small, we obtain that
(68)mS=lim⁡inf⁡t→+∞S(t)≤d+b(1-m)M2Y2β:=N2S,mI=lim⁡inf⁡t→+∞≤rKβ+r(K-d+b(1-m)M2Y2β):=N2I.
Hence, for ɛ>0, sufficiently small, there is a T7≥T6 such that if t>T7, I(t)≥N2I-ɛ. We therefore obtain from the third and the fourth equations of model (2) that for t>T7+τ(69)dY1dt≥pb(1-m)(N2I-ɛ)Y2(t-τ)-(d1+r1)Y1(t),dY2dt=r1Y1(t)-d2Y2(t)-aY22(t).
Consider the following auxiliary equations:
(70)dz1dt=pb(1-m)(N2I-ɛ)z2(t-τ)-(r1+d1)z1(t),dz2dt=r1z1(t)-d2z2(t)-az22(t).
Since (H1) holds, by Lemma 2.4 of [10], it follows from (70) that
(71)lim⁡t→+∞z1(t)=(pb(1-m)(N2I-ɛ)×[pbr1(1-m)(N2I-ɛ)-d2(r1+d1)])×(a(r1+d1)2)-1,lim⁡t→+∞z2(t)=pbr1(1-m)(N2I-ɛ)-d2(r1+d1)a(r1+d1).
By comparison, for ɛ>0, sufficiently small, we obtain that
(72)mY1=lim⁡inf⁡t→+∞Y1(t)≥pb(1-m)N2I[pbr1(1-m)N2I-d2(r1+d1)]a(r1+d1)2:=N2Y1,mY2=lim⁡inf⁡t→+∞Y2(t)≥pbr1(1-m)N2I-d2(r1+d1)a(r1+d1):=N2Y2.
Continuing this process, we derive eight sequences MkS,MkI,MkY1,MkY2,NkS,NkI,NkY1,NkY2(k=1,2,…) such that, for k≥2,
(73)MkS=d+b(1-m)Nk-1Y2β,MkI=rr+Kβ(K-MkS),MkY1=pb(1-m)MkY2MkIr1+d1,MkY2=pbr1(1-m)MkI-d2(r1+d1)a(r1+d1),NkS=d+b(1-m)MkY2β,NkI=rr+Kβ(K-NkS),NkY1=pb(1-m)NkY2NkIr1+d1,NkY2=pbr1(1-m)NkI-d2(r1+d1)a(r1+d1).
It is readily seen that
(74)NkS≤mS≤MS≤MkS,NkI≤mI≤MI≤MkI,NkYi≤mYi≤MYi≤MkYi(i=1,2).
Note that the sequences MkS,MkI,MkY1,MkY2 are nonincreasing and the sequences NkS,NkI,NkY1,NkY2 are nondecreasing. Hence, the limit of each sequence in MkS,MkI,MkY1,MkY2,NkS,NkI,NkY1,NkY2 exists. Denote
(75)lim⁡k→+∞MkS=S¯,lim⁡k→+∞MkI=I¯,lim⁡k→+∞MkYi=Y¯i,(i=1,2),lim⁡k→+∞NkS=S_,lim⁡k→+∞NkI=I_,lim⁡k→+∞NkYi=Y_i,(i=1,2).
From (73), we can obtain
(76)S¯=1β[d+b(1-m)Y_2],I¯=rr+Kβ(K-S¯),Y¯1=pb(1-m)Y2¯I¯r1+d1,Y¯2=pbr1(1-m)I¯-d2(r1+d1)a(r1+d1),S_=1β[d+b(1-m)Y¯2],I_=rr+Kβ(K-S_),Y_1=pb(1-m)Y2_I_r1+d1,Y_2=pbr1(1-m)I_-d2(r1+d1)a(r1+d1).
It follows from (76) that
(77)β(Kβ+r)I¯=r(Kβ-d)+br(1-m)d2a-prr1b2(1-m)2a(r1+d1)I_,(78)β(Kβ+r)I_=r(Kβ-d)+br(1-m)d2a-prr1b2(1-m)2a(r1+d1)I¯,
and (77) minus (78) results in
(79)[β(Kβ+r)-prr1b2(1-m)2a(r1+d1)](I¯-I_)=0.
If aβ(r1+d1)(Kβ+r)≠prr1b2(1-m)2, then we derive from (79) that I¯=I_. It therefore follows from (76) that S¯=S_,Y¯1=Y_1,Y¯2=Y_2. We therefore conclude that E* is globally attractive. The proof is complete.


## 5. Conclusion

In this paper, we have incorporated a prey refuge, stage structure for the predator, and time delay due to the gestation of the predator into an ecoepidemiological predator-prey model. By using Lyapunov functions and the LaSalle invariant principle, the global stability of each of the boundary equilibria of the model is discussed. By using the iteration technique and comparison arguments, sufficient conditions are derived for the global attractivity of the positive equilibrium of the model. By Theorem4, we see that the predator population go to extinction if 0<pbrr1(1-m)(Kβ-d)<βd2(r1+d1)(Kβ+r). By Theorem 6, we see that if pbrr1(1-m)(Kβ-d)>βd2(r1+d1)(Kβ+r) and aβ(r1+d1)(Kβ+r)≠prr1b2(1-m)2, then both the prey and the predator species of model (2) are permanent.



---

*Source: 727818-2014-08-12.xml*