# Analysis of an Ecoepidemiological Model with Prey Refuges

Authors: Shufan Wang; Zhihui Ma
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

An ecoepidemiological system with prey refuges and disease in prey is proposed. Bilinear incidence and Holling III functional response are used to model the contact process and the predation process, respectively. We will study the stability behavior of the basic system from a local to a global perspective. Permanence of the considered system is also investigated.


---

## Body

## 1. Introduction

Ecoepidemiology is the branch of biomathematics that understands the dynamics of disease spread on the predator-prey system. Modeling researches on such ecoepidemiological issues have received much attention recently [1–13]. Anderson and May [1] investigated a prey-predator model with prey infection and observed destabilization due to the spread of infectious diseases within animal and plant communities. Chattopadhyay analyzed predator-prey system with disease in the prey [5] and applied the ecoepidemiological study to the Pelicans at risk in the Salton Sea [6]. Bairagi et al. [13] made a comparative study on the role of prey infection in the stability aspects of a predator-prey system with several functional responses. An ecoepidemiological model with prey harvesting and predator switching was investigated by Bhattacharyya and Mukhopadhyay [14]. Kooi et al. [15] studied stabilization and complex dynamics in a predator-prey system with disease in predator. Most of the above-mentioned studies focused on the role of disease in regulating the dynamical consequences of the interacting populations concerned, such as disease-induced stabilization and destabilization of population states [13].In fact, the dynamical consequences of the predator-prey model can be determined by much ecological effect, such as the Allee effect and prey refuge. Theoretical research and field observations on population dynamics of prey refuges lead to the conclusion that prey refuges have two influences (stabilizing and destabilizing effect) on predator-prey models and prey extinction can be prevented by the addition of prey refuges [13–34]. Here, stabilization (destabilization) of stability refers to cases where an equilibrium point changes from an attractor (a repeller) to a repeller (an attractor) due to increase in the value of a control parameter [17]. Ruxton [16] proposed a continuous-time predator-prey model under the assumption that the rate of prey moving to refuges is proportional to predator density and the results showed that the hiding behavior of prey has a stabilizing effect. The stabilizing effect was also observed in a simple predator-prey system by González-Olivares and Ramos-Jiliberto [17]. Ma et al. [23] formulated a predator-prey model with a class of functional response incorporating the effect of prey refuges and observed the stabilizing and destabilizing effect due to the increases in the prey refuges.In the present research, we formulate a mathematical model of prey-predator interaction with prey refuges and disease in prey. We mainly study the positivity and boundedness, the stability behavior of the disease-free equilibrium point, and the permanence of the basic model.


## 2. Model Formulation

The basic model comprises two population subclasses—(i) prey population with densityN(t) and (ii) predators with density Y(t). To formulate our model, we make the following assumptions.(1)
The prey population increases logistically with intrinsic growth rater and environmental carrying capacity K.(2)
The prey population is divided into two subclasses—the susceptible prey(S(t)) and the infected prey (I(t)) due to infectious disease. We also assumed that N(t)=S(t)+I(t) at any instant of time t.(3)
The susceptible prey is capable of reproducing only and the infected prey is removed by death at a rated1.(4)
The disease is spread only among the prey population and the disease is not genetically inherited. The infected prey does not become immune.(5)
Susceptible prey becomes infected with the simple mass action lawβSI, where β measures the force of infection.(6)
The predators vanish due to natural death at a constant rated2. They consume susceptible and infected prey following the Holling III functional response with predation coefficients c1 and c2, respectively. The consumed prey is converted into predator with efficiency e.(7)
It is assumed that there is a quantityγX of (the susceptible and the infected) prey population incorporating refuges. γ  (0<γ<1) denotes that a constant proportion of (the susceptible and the infected) prey use refuges. All the above-mentioned parameters are assumed to be positive.With the previously mentioned assumptions, the generalized predator-prey system with prey refuges and disease in prey can be represented by the following equations:(2.1)dSdt=rS(1-S+IK)-βSI-c1(1-γ)2S2Ym+(1-γ)2S2,dIdt=βSI-c2(1-γ)2I2Ym+(1-γ)2I2-d1I,dYdt=ec1(1-γ)2S2Ym+(1-γ)2S2+ec2(1-γ)2I2Ym+(1-γ)2I2-d2Y.DefiningG={(S,I,Y)∈R+3∣S+I≤K,S≥0,I≥0,Y≥0}, it is easy to show that the set G is the positively invariable set of system (2.1).


## 3. The Positivity and Boundedness

Theorem 3.1.
All solutions of system (2.1) initiating R+3 are positive and ultimately bounded.Proof.
Let(S(t),I(t),Y(t)) be one of the solutions of system (2.1).
Integrating (2.1) with initial conditions (S0,I0,Y0), we have
(3.1)S(t)=S0exp[∫0t‍(r(1-S(ξ)+I(ξ)K)-βI(ξ)-c1(1-γ)2S(ξ)Y(ξ)m+(1-γ)2S2(ξ))dξ]>0,I(t)=I0exp[∫0t‍(βS(ξ)-c2(1-γ)2I(ξ)Y(ξ)m+(1-γ)2I2(ξ)-d1)dξ]>0,Y(t)=Y0exp[∫0t‍(ec1(1-γ)2S2(ξ)m+(1-γ)2S2(ξ)+ec2(1-γ)2I2(ξ)m(1-γ)2I2(ξ)-d2)dξ]>0.
Hence all solutions starting inR+3 remain in R+3 for all t≥0.Next, we will prove the boundedness of the solutions.BecausedS/dt≤rS(1-S/K), then we have
(3.2)limt→+∞supS(t)≤K.LetW=eS+eI+Y, then we obtain that
(3.3)dWdt=reS(1-SK)-reKSI-ed1I-d2Y≤reS-ed1I-d2Y=(r+1)eS-(eS+d1eI+d2Y)≤e(r+1)K-δW,
where δ=min{1,d1,d2}.Hence, we have(3.4)dWdt+δW≤e(r+1)K.
That is
(3.5)0≤W(S,I,Y)≤e(r+1)Kδ+W(S(0),I(0),Y(0))eδt,0≤W≤e(r+1)Kδ(t→+∞).Thus, all curves of system (2.1) will enter the following region:
(3.6)B={(S,I,Y):0≤W≤e(r+1)Kδ+ε,∀ε>0}.


## 4. The Equilibrium Point

All equilibrium points of system (2.1) can be obtained by solving the following equations:
(4.1)rS(1-S+IK)-βSI-c1(1-γ)2S2Ym+(1-γ)2S2=0,βSI-c2(1-γ)2I2Ym+(1-γ)2I2-d1I=0,ec1(1-γ)2S2Ym+(1-γ)2S2+ec2(1-γ)2I2Ym+(1-γ)2I2-d2Y=0,
These points are as follows:(1)
the trivial equilibrium pointE0(0,0,0),(2)
the equilibrium pointEK(K,0,0),(3)
the predator-extinction equilibrium pointE-(S-,I-,0),(4)
the disease-free equilibrium pointE*(S*,0,Y*), where
(4.2)S-=d1β=˙SM,I-=r(βK-d1)β(r+βK)=˙IM,S*=11-γd2mec1-d2,Y*=erS*d2(1-S*K).LetR1=(ec1(1-γ)2K2/(m+(1-γ)2K2))/d2. It is clear to show that the disease-free equilibrium point E*(S*,0,Y*) has its ecological meaning when R1>1.LetR0=βK/d1, the predator-extinction equilibrium point E-(S-,I-,0) is nonnegative when R0>1.


## 5. The Stability Property

In this section, we will study the local and global stability of the equilibrium points of system (2.1).Theorem 5.1.
LetR-1=((ec1(1-γ)2SM2/(m+(1-γ)2SM2))+(ec2(1-γ)2IM2/(m+(1-γ)2IM2)))/d2, R-2=max(R1,R-1), then one has the following. (1)
The trivial equilibrium pointE0(0,0,0) is always unstable.(2)
IfR0<1, R1<1, then the predator-extinction equilibrium point E-(S-,I-,0) is globally asymptotically stable. If R0>1 or R1>1, the equilibrium point EK(K,0,0) is unstable.(3)
IfR0>1 and R-1<1, the predator-extinction equilibrium point E-(S-,I-,0) is locally asymptotically stable. If R-2<1, the predator-extinction equilibrium point E-(S-,I-,0) is locally asymptotically stable. If R-1>1, the predator-extinction equilibrium point E-(S-,I-,0) is unstable.Proof.
The Jacobian matrix of system (2.1) at the trivial equilibrium point E0(0,0,0) is
(5.1)(r000-d1000-d2).Clearly, the trivial equilibrium pointE0(0,0,0) is unstable.The Jacobian matrix of system (2.1) at the equilibrium point EK(K,0,0) is
(5.2)(-r-r-βK-c1(1-γ)2K2m+(1-γ)2K20βK-d1000ec1(1-γ)2K2m+(1-γ)2K2-d2).According to the theorem about the local stability, the local stability of the equilibrium pointEK(K,0,0) is determined only by the sign of (ec2(1-γ)2K2/(m+(1-γ)2K2))-d2 and βK-d1.Therefore, ifR0<1 and R1<1, the equilibrium point EK(K,0,0) is locally asymptotically stable.Next, we will prove the global stability of the equilibrium pointEK(K,0,0).Defining the Lyapunov functionV=eI+Y, then we obtain that
(5.3)dVdt=e(βS-d1)I+(ec1(1-γ)2S2m+(1-γ)2S2-d2)Y≤e(βK-d1)I+(ec1(1-γ)2K2m+(1-γ)2K2-d2)Y.Hence, ifR0<1 and R1<1, dV/dt≤0.HavingE={(S,I,Y)∈G∣dV/dt=0}={I=0,Y=0}, the maximum invariable set of system (2.1) is M=E={I=0,Y=0}.According to the LaSalle invariable set theorem,limt→+∞I(t)=0, limt→+∞Y(t)=0.Thus, the limit equation of system (2.1) is
(5.4)dSdt=rS(1-SK).Clearly, the equilibrium pointS=K is globally asymptotically stable.According to the limit system theorem, ifR0<1 and R1<1, the equilibrium point EK(K,0,0) is globally asymptotically stable.SettingR0>1, the Jacobian matrix of system (2.1) at the predator-extinction equilibrium point E-(S-,I-,0) is
(5.5)(-rS-K-rS-K-βS--c1(1-γ)2S-2m+(1-γ)2S-2βI-0-c2(1-γ)2I-2m+(1-γ)2I-200ec1(1-γ)2S-2m+(1-γ)2S-2+ec2(1-γ)2I-2m+(1-γ)2I-2-d2).The characteristic equation of system (2.1) at the predator-extinction equilibrium point E-(S-,I-,0) is
(5.6)(λ-A)(λ2-Bλ+C)=0,
where
(5.7)A=c1(1-γ)2SM2m+(1-γ)2SM2+ec2(1-γ)2IM2m+(1-γ)2IM2-d2,B=-2SMK<0,C=βIM(SMK+βSM)>0.According to the Routh-Hurwitz rule, the predator-extinction equilibrium pointE-(S-,I-,0) is locally asymptotically stable when R-1<1.Next, we will prove the global stability of the predator-extinction equilibrium pointE-(S-,I-,0).Defining the Lyapunov functionV=Y, then we have
(5.8)dVdt=(ec1(1-γ)2S2m+(1-γ)2S2+ec2(1-γ)2I2m+(1-γ)2I2-d2)Y.LetΦ(S,I)=(c1(1-γ)2S2/(m+(1-γ)2S2))+(ec2(1-γ)2I2/(m+(1-γ)2I2)) in which (S,I)∈D={(S,I)∣S+I≤K,S≥0,I≥0}.Thus(5.9)Φ′(t)=2c1em(1-γ)2S(m+(1-γ)2S2)2dSdt+2c2em(1-γ)2I(m+(1-γ)2I2)2dIdt.Clearly,(2c1em(1-γ)2S/(m+(1-γ)2S2)2)>0 and (2c2em(1-γ)2I/(m+(1-γ)2I2)2)>0 on the set D'=D-{S=0∣I=0}.Thus,Φ'(t)=0 on the set D', if and only if dS/dt=0 and dI/dt=0, then S=SM, I=IM.IfS=I=0, then Φ(S,I)=0. If S=SM, then I=IM. If S<SM, then dI/dt<0; that is, I<IM.Therefore, the maximum value of the functionΦ(t) is obtained at the point (K,0) or (SM,IM).Hence, ifR-2<1, dV/dt≤0; that is, E={(S,I,Y)∈G∣dV/dt=0}={Y=0}.According to the LaSalle invariable set theorem,limt→+∞Y(t)=0. Thus, the limit equation of system (2.1) is
(5.10)dSdt=rS(1-SK)-βSI,dIdt=βSI-d1I.According to the results of the appendix section, ifR0>1, the equilibrium point (SM,IM) is globally asymptotically stable.According to the limit system theorem, ifR0>1 and R-2<1, the predator-extinction equilibrium point E-(S-,I-,0) is globally asymptotically stable.Theorem 5.2.
LetR1*=(β/d1(1-γ))(d2m/(ec1-d2)), R2*=(((2d2-ec1)(1-γ)K)/2d2)(ec1-d2)/d2m. If R1>1, then the disease-free equilibrium point E*(S*,0,I*) is nonnegative, (1)
ifR1*<1 and R2*<1, then the disease-free equilibrium point E*(S*,0,I*) is locally asymptotically stable;(2)
ifR0<1 and R2*<1, then the disease-free equilibrium point E*(S*,0,I*) is globally asymptotically stable;(3)
ifR1*>1 or R2*>1, then the disease-free equilibrium point E*(S*,0,I*) is unstable.Proof.
AssumingR1>1, the Jacobian matrix of system (2.1) at the disease-free equilibrium point E*(S*,0,Y*) is
(5.11)(a11-βS*-rS*K-c1(1-γ)2(S*)2m+(1-γ)2(S*)20a2202ec1m(1-γ)2S*Y*(m+(1-γ)2(S*)2)200),
where
(5.12)a11=r(1-2S*K)-2c1m(1-γ)2S*Y*(m+(1-γ)2(S*)2)2,a22=βS*-d1.The characteristic equation of system (2.1) at the disease-free equilibrium point E*(S*,0,Y*) is
(5.13)(λ-a22)(λ2-a11λ+2c12em(1-γ)4(S*)3Y*(m+(1-γ)2(S*)2)3)=0.Clearly, ifR1*<1, a22=βS*-d1<0.Again, we have(5.14)a11=r(1-2S*K)-2c1m(1-γ)2S*Y*(m+(1-γ)2(S*)2)2=r(1-2S*K)-2r(1-d2ec1)(1-S*K)=r[1-2S*K-2+2S*K+2d2ec1-2d2ec1S*K]=rec1[(2d2-ec1)-2d2(1-γ)Kd2mec1-d2].Hence, ifR2*<1, then a11<0.According to the Routh-Hurwitz rule, ifR1*<1 and R2*<1, then the disease-free equilibrium point E*(S*,0,Y*) is locally asymptotically stable.Next, we will study the global stability of the disease-free equilibrium pointE*(S*,0,Y*).Defining the Lyapunov functionV=I, then we obtain
(5.15)dVdt=(βS-c2(1-γ)2Ym+(1-γ)2I-d1)I≤(βK-d1)I.Thus, ifR0<1, dV/dt≤0 and dV/dt=0 if and only if I=0; that is, E={(S,I,Y)∈G∣dV/dt=0}={I=0}.According to the LaSalle invariable set theorem,I(t)→0 when t→+∞. The limit system of system (2.1) is
(5.16)dSdt=rS(1-SK)-c1(1-γ)2S2Ym+(1-γ)2S2,dYdt=ec1(1-γ)2S2Ym+(1-γ)2S2-d2Y.Clearly, the equilibrium points of the system (5.16) are
(5.17)E0′(0,0),EK′(K,0),E*′(S*,Y*)
in which S* and Y* are similar as the equilibria expression of system (2.1).It is easy to show that the equilibrium pointE*′ has its ecological meaning when R1>1.According to the Routh-Hurwitz rule, the equilibrium pointE0′ is unstable. If R1>1, then EK′(K,0) is unstable.Again, the Jacobian matrix of system (5.16) at the equilibrium point E*′ is
(5.18)(a11-c1(1-γ)2(S*)2m+(1-γ)2(S*)22ec1m(1-γ)2S*Y*(m+(1-γ)2(S*)2)20),
where
(5.19)a11=r(1-2S*K)-2c1m(1-γ)2S*Y*(m+(1-γ)2(S*)2)2.The characteristic equation of system (5.16) at the equilibrium point E*′ is
(5.20)λ2-a11λ+2c12em(1-γ)4(S*)3Y*(m+(1-γ)2(S*)2)3=0.According to the above study, ifR2*<1, then a11<0.Hence, ifR2*<1, then the equilibrium point E*′ is locally asymptotically stable in the region D by the Routh-Hurwitz rule.It is easy to note that the globally asymptotically stability of the equilibrium pointE*′ implies that there is no close orbit in the region D for the considered system.Let(5.21)S=d2mec1-d2S~,Y=1c1d2m(ec1-d2)Y~,dt=d2(1-γ)2S~2+ec1-d2d2(ec1-d2)dt~,
and rewrite S~, Y~, and t~ into S, Y, and t, then system (5.16) becomes as follows:
(5.22)dSdt=S(A0+A1S+A2S2+A3S3)-(1-γ)2S2Y≡P(S,Y),dYdt=(1-γ)2S2Y-Y≡Q(S,Y),
where A0=r/d2>0, A1=-r/(d2K)d2m/(ec1-d2)<0, A2=(1-γ)2r/(ec1-d2)>0, A3=-(1-γ)2r/(ec1-d2)d2m/(ec1-d2)<0.Thus, the positive equilibrium pointE*′(S*,Y*) of system (5.16) becomes the positive equilibrium point EP(SP,YP) of system (5.22), where
(5.23)SP=11-γ,YP=11-γec1d2(ec1-d2)(r-r(1-γ)Kd2mec1-d2)=rec1d2(ec1-d2)(1-γ)(1-1(1-γ)Kd2mec1-d2).Considering the Dulac functionB(S,Y)=S-2Yn-1, then we have
(5.24)∂(BP)∂S+∂(BQ)∂Y=S-2Yn-1[2A3S3+(A2+n(1-γ)2)S2-(A0+n)]≡S-2Yn-1Φ(S,n).In order to prove the global stability, we will prove only that there exists a real numbern such that Φ(S,n)≤0.Clearly, ifS=0 or S=-(A2+n(1-γ)2)/3A3, then 2S[3A3S+A2+n(1-γ)2]=0 since ∂Φ(S,n)/∂S=2S[3A3S+A2+n(1-γ)2].Let-(A2+n(1-γ)2)/3A3>0  (A2>0, A3<0), that is, A2+n(1-γ)2>0, then we obtain that if 0<S<-(A2+n(1-γ)2)/3A3, then ∂Φ(S,n)/∂S>0. Otherwise, ∂Φ(S,n)/∂S<0.Therefore, the functionΦ(S,n) has the maximum value at the point S=-(A2+n(1-γ)2)/3A3 and Φ(0,n)=-(A0+n)<0.Hence, there exists only one real numbern, such that A2+n(1-γ)2>0 and
(5.25)max0<S<+∞Φ(S,n)=Φ(S,n)∣S=-(A2+n(1-γ)2)/3A3=[A2+n(1-γ)2]327A32-A2+n(1-γ)2(1-γ)2+A2(1-γ)2-A0=0.Thus, we will prove only that there existsA2+n(1-γ)2>0, such that
(5.26)[A2+n(1-γ)2]3+M[A2+n(1-γ)2]+N=0,
where M=-27A32/(1-γ)2<0, N=27A32[A2/(1-γ)2-A0]>0.The discriminant of the cubic equation is(5.27)Δ=N24+M227=27A344[(A2(1-γ)2-A0)2-4A32(1-γ)6]=[A2Sp2-A0+2A3Sp3][A2Sp2-A0-2A3Sp3].It is easy to show that ifd2<ec1<2d2, then A2Sp2-A0=r(2d2-ec1)/d2(ec1-d2)>0.Again, ifR1>1 and R2*<1, then A2Sp2-A0-2A3Sp3>0 and 2A3Sp3+A2Sp2-A0≤0. If 2A3Sp3+A2Sp2-A0=0, then Δ=0.According to Shengjin's distinguishing means, the cubic equation has one negative real root and two positively real roots.If2A3Sp3+A2Sp2-A0<0, the cubic equation has three roots which are not equal.According to the Descartes rule of signs, the cubic equation has two positively real roots and one negatively real root at most.Therefore, the cubic equation has at least one positively real root. That is to say, there exists a numbern, such that A2+n(1-γ)2>0.Furthermore, we obtain that(5.28)∂(BP)∂S+∂(BQ)∂Y≤0.According to the Bendixson-Dulac theorem, there does not exist the limit cycle for the limit system.Hence, the equilibrium pointE*′ is globally asymptotically stable.Therefore, ifR0<1 and R2*<1, then the equilibrium point E*′ is globally asymptotically stable according to the limit system theorem.


## 6. Permanence

Theorem 6.1.
IfR0>1, R1>1, R-1>1, and R1*>1, then system (2.1) is permanent.Proof.
Considering the average Lyapunov functionV(S,I,Y)=Sα1Iα2Yα3, where αi(i=1,2,3) is positive, then in the region R+3, we have
(6.1)V˙V=α1[r(1-S+IK)-βI-c1(1-γ)2SYm+(1-γ)2S2]+α2[βS-c2(1-γ)2IYm+(1-γ)2I2-d1]+α3[ec1(1-γ)2S2m+(1-γ)2S2+ec2(1-γ)2I2m+(1-γ)2I2-d2]≡Ψ(S,I,Y).In order to prove the permanence of system (2.1), we only indicate the following results: the function Ψ(S,I,Y)>0 for all boundary equilibrium points.Letα1>(d1α2+d2α3)/r, then
(6.2)Ψ(E0)=α1r-α2d1-α3d2>0,Ψ(EK)=α2(βK-d1)+α3(ec1(1-γ)2K2m+(1-γ)2K2-d2),Ψ(E-)=α3(ec1(1-γ)2SM2m+(1-γ)2SM2+ec2(1-γ)2IM2m+(1-γ)2IM2-d2),Ψ(E*)=α2(βS*-d1).Hence, ifR0>1 and R1>1, then Ψ(EK)>0. If R-2>1, then Ψ(E~)>0. If R1*>1, then Ψ(E*)>0.Therefore, system (2.1) is permanent when R0>1, R1>1, R-1>1, and R1*>1 by the average Lyapunov function theorem [34].By simple computation,R1>1(<1) is equivalent to (1-γ)K(ec1-d2)/d2m>1(<1).Hence, ifR1>1, then R1*>1 implies R0>1.



---

*Source: 371685-2012-11-21.xml*