A predator–prey model with disease in the predator species only

Nonlinear Analysis: Real World Applications 11 (2010) 2224–2236

Contents lists available at ScienceDirect

Nonlinear Analysis: Real World Applications journal homepage: www.elsevier.com/locate/nonrwa

A predator–prey model with disease in the predator species only Mainul Haque ∗ School of Mathematical Science, University of Nottingham, University Park, Nottingham, NG7 2RD, UK

article

info

Article history: Received 9 November 2008 Accepted 15 June 2009 Keywords: Eco-epidemiological model Infected predator SIS epidemic model Stability and bifurcation analysis

abstract In this paper we propose a predator–prey model with logistic growth in the prey population that includes an SIS parasitic infection in the predator population, with the assumption that the predator has an alternative source of food. For simplicity we initially work with a model involving the fractions of the predator which are susceptible and those infected and then translate the results back to the model with absolute numbers. Important thresholds R10 , R20 , R30 and R40 are identified and their implications have been explained. Our theoretical study indicates that the absence of prey may be beneficial for predator when a transmissible disease runs among the predator population. One important conclusion is that infection in the predator species may save the prey from extinction even if R20 , the basic reproduction number for the prey to be able to invade the predator-only equilibrium, is less than one. Therefore infection in the predator species may be taken as biological control. Finally, analytical results have been supported by numerical simulations with the help of experimental data. © 2009 Elsevier Ltd. All rights reserved.

1. Introduction Mathematical and statistical models can help us to identify key parameters which determine the rich dynamics of an ecological or epidemiological system. In the development of quantitative theory for interaction of predator and prey, mathematical ecology is also an important factor along with the experimental ecology, see [1] and references therein. Transmissible disease in an ecological situation is fast becoming a major field of study in its own right. The first mathematical description of contagious diseases has been formulated by Kermack and McKendric [2]. Both theoretical and experimental investigations in these two fields namely ecology and epidemiology progressed independently along the years, until the late eighties and early nineties. The ecological literature has increasingly emphasized the importance of parasites in shaping the dynamics of both plant and animal communities, [3–5]. Some models merging features of the two phenomena, i.e. the demographics of interacting species and an epidemic evolution in such a composite environment, were then considered, see [6–17] and the references therein. Most of the above models are based on the assumption that the infected prey is more vulnerable to predation except for those developed by Venturino [11], and Haque and Venturino [15]. In these two papers the authors considered the case when disease is spreading among the predator population. The model considered here also studies the situation where an epidemic runs among the predator population but differs from previous models; it includes the recovery rate and also takes into account also the stability of the positive interior equilibrium point. In addition, we assume that the predator has a logistic growth rate since it has sufficient resources for alternative foods; and it is argued that alternative food sources may have an important role in promoting the persistence of predator–prey systems. van Baalen et al. [18] have observed that the switching of the predator species has a significant contribution to the persistence of species in the predator–prey system. They also have analyzed the conditions for stability as well as long-term behavior of the system under bounded

∗

Tel.: +44 0 115 951 4949; fax: +44 0 115 951 4951. E-mail address: [email protected].

1468-1218/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.nonrwa.2009.06.012

M. Haque / Nonlinear Analysis: Real World Applications 11 (2010) 2224–2236

2225

population fluctuations. This leads to specific hypotheses about which types of alternative food (in terms of nutritional quality, availability, and handling time) promote persistence and by which mechanism. 2. Predator–prey model with disease in the predator species only Let us consider a predator–prey model, where P (t ) represents the number of prey and H (t ) denotes the number of predators at the time t. The classical well known model is given by dP dt

= aP 1 −

!

P

dH

− cPH ,

M

dt

= rH 1 −

H

!

K

+ ecPH .

(2.1)

Here a is the growth rate of the prey and M is the prey carrying capacity of the environment. The total predation rate is cHP and the conversion rate is denoted by e. Again r represents the growth rate of the predator and K is the predator carrying capacity of the environment. Now we consider the case when a transmissible disease spreads among the predator population. There are several examples in Table 2 of [19]; an example is common seal (Phoca vitulina) and striped dolphin (Stenella coeruleoaba) when they are infected by Phocine Distemper Virus (PDV). The epidemic of PDV infections are found among harbor seals in 1988. PDV was first identified in 1988 as the cause of death of 18,000 harbor seals (P. vitulina) and 300 gray seals (Halichoerus grypus) along the northern European coast; in these case fish are considered as prey. Therefore we make the strong assumption that the disease does not cross the species barrier. Consideration of the case where the infection crosses the species barrier would make the model less amenable to mathematical analysis. The following basic assumptions have been made to construct our model: (A1) In the absence of prey, the predator population grows logistically with per capita birth rate b − (θ rH /K ) and per capita death rate d + ((1 − θ )rH /K ), where r = b − d > 0 and 0 ≤ θ ≤ 1. The restricted growth in the logistic equation is due to a density-dependent per capita death rate when θ = 0, to a density-dependent per capita birth rate when θ = 1, and to a combination of these two if 0 < θ < 1. (A2) In the presence of disease the predator population is divided into two classes, susceptible predator, X (t ) and infected predator, Y (t ). The total population at time t is H (t ) = X (t )+ Y (t ). Then the fraction of predators that is susceptible is S = HX and the fraction of predators that is infectious is I = HY . (A3) The number of new infected predators per unit time increases due to the infection of susceptible predators through their contacts with infectious predators. If the average number of contacts of a susceptible predator per unit time is denoted βY by I = HY , then H = β I is the per capita contact rate of a susceptible predator with infectious predator per unit time, β XY

therefore H = β ISH is the total number of new cases per unit time. This type of disease incidence function is called standard incidence and is formulated by Anderson and May [20]. (A4) We assume that reproduction takes place also among infected predators, possibly at a different birth rate h than for healthy ones, and the newborns are free from disease. On the basis of the above assumptions we propose the following models dP dt dH dt dX dt dY dt

= aP 1 −

M H

= rH 1 −

= b−

!

P

+ ecXP + ecqYP ,

(2.2b)

!

K

= hY 1 −

(2.2a)

!

K

θ rH

− cXP − cqYP ,

Y

H−

! +

K

d+

β XY H

(1 − θ )rH K

! X−

− γY − d +

β XY H

+ γ Y + ecXP ,

(1 − θ )rH K

(2.2c)

! Y + ecqYP .

(2.2d)

Here note that, b − d = r, H = X + Y , X = susceptible predator, Y = infected predator, c is the searching coefficient, e is the conversion factor. The parameter q ≤ 1 reflects the situation where the infected predators are less efficient to catch prey compared to the healthy predators. Hence the predation rate for the sick predators measured here by the the parameters cq which is smaller than that of healthy individuals c. For sake of simplicity, we put h = 0 in (2.2) to give dP dt

= aP 1 −

P M

! − cXP − cqYP ,

(2.3a)

2226

M. Haque / Nonlinear Analysis: Real World Applications 11 (2010) 2224–2236

dH dt dX dt dY dt

= rH 1 −

= b− β XY

=

H

H

! + ecXP + ecqYP ,

K

θ rH

! H−

K

− d+

d+

(2.3b)

(1 − θ )rH

! X−

K

(1 − θ )rH

β XY

+ ecXP ,

H

(2.3c)

! Y + ecqYP .

K

(2.3d)

In this article, we discuss the dynamical behavior of the system (2.3). For this purpose we put H = X + Y , I = in the model (2.3) that leads to the following set of differential equations: dP dt dH dt dI dt

=P a 1−

=H r 1−

P

!

!

(2.4a)

! + ecP (1 − I + qI ) ,

K

(2.4b)

! β + ecP (q − 1) (1 − I ) − γ + b −

=I

and S = 1 − I

! − cH (1 − I + qI ) ,

M H

Y H

θ rH

!!

K

.

(2.4c)

The model is mathematically well posed in the positively invariant region

Ω = {(P , H , I ) : 0 ≤ P ≤ M , 0 ≤ H ≤ K1 , 0 ≤ I ≤ 1}, where K1 = K (1 + ecM (1 + q)/r ). Moreover any solution with positive initial values P0 H0 , and I0 ultimately approaches this region (Appendix). 2.1. Equilibrium points and their stability The non-negative equilibria of the system (2.4) are (i) E0 = (0, 0, 0), (ii) E1 = (M , 0, 0), (iii) E2 = (0, K , 0), (iv) E3 = β rM (a−cK ) aK (r +ecM ) , (v) E4 = (P4 , H4 , 0) where P4 = ar +ec 2 MK and H4 = ar +ec 2 MK , (vi) E5 = (0, K , I5 ), γ +b

(0, 0, I3 ), where I3 = 1 − R11 and R10 = 0

where I5 = 1 − where

1 R30

β

and R30 = b+γ −r θ , (vii) E6 = (M , 0, I6 ) where I6 = 1 −

I ∗ = 1 − ψ ∗,

H∗ =

aKr + MKeca(ψ + q − qψ) ar +

Kec 2 M

(ψ + q − qψ)

2

,

P∗ =

1 R40

and R40 =

β+ecM (q−1) β+γ

and (viii) E7 = (P ∗ , H ∗ , I ∗ )

Mra − crKM (ψ + q − qψ) ar + Kec 2 M (ψ + q − qψ)2

and ψ is the unique real positive root in ∗

(I)

(II)

0, min 1,

a − cqK

!! if q < 1 with

(1 − q)cK

cqK − a aKr + ecaqMK , min 1, (q − 1)cK eca(q − 1)MK

(III) (0, 1)

if q = 1 with

a cqK

a cqK

> 1,

!! if q > 1 with

a cqK

< 1 and

> 1;

of the cubic equation Aψ 3 + 3Bψ 2 + 3C ψ + D = 0,

(2.5)

where A = β Kec 2 M (1 − q)2 , 3B = Kec 2 M (1 − q)[−b − γ + q(b + γ + 2β)], 3C = Kec 2 Mq[(1 − q)(−2b − 2γ + r ) + β q] + ar [ecM (θ − 1)(1 − q) + β], D = −ar (b + γ − r θ ) − ecqM (cqbK + cqγ K − ar θ ). The above equation has exactly one real positive root if G2 + 4H 3 > 0, where G = A2 D − 3ABC + 2B3 , H = AC − B2 , and

√

1

by Cardano’s method, the root is A1 (w ∗ − wH∗ − B), where w ∗ denotes one of the three values of [ 12 (−G + G2 + 4H 3 ) 3 ]. Now the equilibrium points E0 , E1 and E2 always exist. E3 is admissible if and only if R10 > 1 and E4 is admissible if and only a if R20 = cK > 1. The equilibrium point E5 admissible if and only if R30 > 1. E6 is admissible if and only if R40 > 1. 2.2. Ecological interpretation dI R10 > 1 is the condition for the disease to spread ( dt > 0) if a small infected fraction of the predator is introduced into the predator–prey extinction equilibrium E0 . The expected number of secondary cases produced by each newly infected predator β during its entire infected lifetime is d+γ . However in the neighborhood of E0 , the predator population is exponentially

M. Haque / Nonlinear Analysis: Real World Applications 11 (2010) 2224–2236

2227

increasing at rate r = b − d. R10 is this expected number adjusted for the exponential increase [21]. We need to make this adjustment as I is a fraction with an increasing denominator so that each infected predator must do more than just replace itself for the infection fraction to persist. The threshold R20 is the expected absolute number of offspring of each prey during its lifetime introduced into the predator-only equilibrium E2 . R20 > 1 is the condition for the prey population to take off (increase from zero) in this situation. Similarly R30 is the absolute number of secondary infected predators caused by a single newly infected predator entering the susceptible predator-only equilibrium E2 . R30 > 1 is the condition for the predator population to take off (increase from zero) in this situation. β β R40 does not have a natural similar interpretation; but we can write I6 = β+ecM (q−1) (1 − 14∗ ) where R40∗ = γ +b−ecM (q−1) . R0

The conditions R40 > 1 and R40∗ > 1 are equivalent for the existence of the endemic equilibrium I6 . Now consider a single newly infected predator entering the prey-only equilibrium E1 . R40∗ is the expected number of secondary cases produced by each newly infected predator during its entire infected lifetime adjusted for the exponential increase (or decrease) of the predator population in the neighborhood of this equilibrium. 2.3. Behavior of the system when either prey, predator or disease are absent Note that the paths on the P-axis with P (0) > 0 tend to equilibrium E1 as time becomes large i.e. the prey goes to its carrying capacity when predator is absent. Similarly, the paths on the H-axis with H (0) > 0 tend to the equilibrium E2 as time becomes large, i.e. the predator goes to its carrying capacity when both the disease and the prey are absent. The solution paths on the I-axis with I (0) > 0 go to E3 if R10 > 1 and E0 if R10 ≤ 1. Note that it is natural to interpret E3 and E0 as the same equilibrium for P = H = X = Y = 0 in the absolute system (2.3), although the direction in which this equilibrium is approached is different in the two cases. The solution paths in Ω with P (0) = 0, H (0) > 0 and I (0) > 0 tend to E5 if R30 > 1 and E2 when R30 ≤ 1, i.e. in the absence of the prey the predator goes to its carrying capacity, and the disease persists if and only if its basic reproduction number at E2 exceeds unity. Similarly if P (0) > 0, H (0) = 0 and I (0) > 0 when R40 > 1, the solution paths of the system go to E6 , but to E1 if R40 ≤ 1. Note that again E1 and E6 naturally correspond to the same equilibrium P = M , H = X = Y = 0 in the absolute system (2.3). Finally if P (0) > 0, H (0) > 0 and I (0) = 0 then the solution paths go to E2 if R20 ≤ 1. We conjecture that they go to E4 if R20 > 1 but have not been able to prove this. 3. Stability analysis Here we perform the stability and bifurcation analysis of the system (2.4). System (2.4) can be written in the form ˙ = F (X) where X = (P , H , I )T , and the Jacobian matrix of the system is X j11  ecH (1 − I + qI )

 J=

ecI (1 − I )(q − 1)

−cP (1 − I + qI ) j22 rθI

 −cPH (q − 1) ecPH (q − 1)  ,

K

j33

where j11 = a(1 − 2P /M ) − cH (1 − I + qI ), j22 = r (1 − 2H /K ) + ecP (1 − I + qI ) and j33 = (β + ecP (q − 1))(1 − 2I ) − (γ + b − θ rH /K ). 3.1. Dynamics of the system around E0 , E1 , E3 , E6 , and E2 Straightforward calculations of eigenvalues at equilibrium points E0 , E1 , E3 , and E6 show that at least one of them is strictly positive in each case; therefore they are always unstable. The behavior of the system around E2 can be summarized in the following lemma. Lemma 1. (i) The hyperbolic equilibrium point E2 is locally asymptotically stable (LAS) if R20 < 1 and R30 < 1 and unstable if either R20 > 1 or R30 > 1. (ii) If R20 < 1 and R30 < 1, E2 is globally asymptotically stable (GAS) provided that H (0) > 0 and ξ1 > 0, where ξ1 = (b + γ − r θ) − (a + β + cK1 + ecM (θ + 1)(q + 1)). Proof. (i) The eigenvalues at the equilibrium points E2 are cK (R20 − 1), −r and (b + γ − r θ )(R30 − 1). Hence the proof. (ii) Consider d dt

dI +P , dt dt ≤ PI (a + cK1 + β + ecM (q + 1) − (b + γ ) + r θ + eθ cM (q + 1)),

(PI ) = I

dP

= −ξ1 PI , where ξ1 = (b + γ − r θ ) − (a + β + cK1 + ecM (θ + 1)(q + 1)) > 0. Hence PI → 0 as t → ∞.

2228

M. Haque / Nonlinear Analysis: Real World Applications 11 (2010) 2224–2236

So using the results of the Appendix and considering H1 dH we deduce that given 1 > 0 there exists t0 such that for dt t ≥ t0 , IP ≤ 1 , P ≤ M , I ≤ 1 and K1 ≥ H ≥ K − 1 . Hence for t ≥ t0 , dP

≤ (aP + cHIP − c (K − 1 )P ),

dt d dt

≤ P (a − cK ) + c 1 (K1 + M ), Pe(cK −a)t ≤ c 1 (K1 + M )e(cK −a)t , c 1 (K1 + M )

Pe(cK −a)t − P (t0 )e(cK −a)t0 ≤ P (t ) ≤ P (t0 )e−(cK −a)(t −t0 ) +

(e(cK −a)t − e(cK −a)t0 ),

cK − a c 1 (K1 + M ) cK − a

(1 − e−(cK −a)(t −t0 ) ).

Therefore there exists t1 ≥ t0 such that for t ≥ t1 , 0 ≤ P (t ) ≤ t → ∞. So given 2 > 0 there exists t2 such that for t ≥ t2 , 1 dH H dt

≤r 1−

H

!

H

+ ecP + ecqPI ≤ r 1 −

K

We deduce that lim supt →∞ H ≤ K 1 +

ec (q1 +2 ) r

2c 1 (K1 +M ) . cK −a

But 1 > 0 is arbitrary. Hence P (t ) → 0 as

!

K

+ ec 2 + ecq1 .

! . So H → K as t → ∞.

Hence given 3 > 0 with b + γ > β + ec 3 + r θ 1 +

3

!

K

there exists t3 such that for t ≥ t3 , IP ≤ 3 and H ≤ K + 3 .

For t ≥ t3 , dI dt



≤ I β + ec 3 + r θ 1 +

3  K

! − (b + γ ) + ecq3 .

A similar argument to that used to show that P → 0 as t → ∞ now shows that I → 0 as t → ∞.



3.2. Dynamics of the system around E4 = (P4 , H4 , 0) Lemma 2. (i) E4 is LAS if β < β [4] (R50 < 1) and unstable if β > β [4] (R50 > 1), where β [4] = γ + b + ec (1 − q)P4 − R50

β . β [4 ]

r θ H4 K

and

= (ii) The hyperbolic equilibrium point E4 = (P4 , H4 , 0) transforms into a non-hyperbolic equilibrium point at β = β [4] and the system (2.4) experiences a transcritical bifurcation as β passes through β [4] i.e for R50 = 1. Proof. (i) At E4 the eigenvalues of the Jacobian matrix are ω1 = β − (γ + b + ec (1 − q)P4 − roots of

ω + 2

rH4 K

+

aP4 M

! ω+

ra KM

r θ H4 K

) = β [4] (R50 − 1) and the

! + ec

2

H4 P4 = 0.

(3.1)

Hence as both the roots of (3.1) have negative real parts the result follows. (ii) By following a similar technique to [17], we checked that the system undergoes a transcritical bifurcation as β passes through β [4] .  3.3. Dynamics of the system around E5 = (0, K , I5 ) Suppose that R30 > 1 so that E5 exists. Lemma 3. Define a5 = cK (1 + (q − 1)I5 ). (i) If a < a5 then the system (2.4) will be LAS about E5 and if a > a5 then the system (2.4) will be unstable about E5 . (ii) If a < a5 then the system (2.4) will be GAS about E5 if H (0) > 0, I (0) > 0 and ξ1 > 0. (iii) The hyperbolic equilibrium point E5 = (0, K , I5 ) transforms into a non-hyperbolic equilibrium point at a = a5 and the system (2.4) undergoes a transcritical bifurcation as a passes through a5 . Proof. (i) At E5 , the eigenvalues of the Jacobian matrix are a − a5 , −r and −β I5 , and. The result follows. d (ii) Arguing as in the proof of Lemma 1(ii), we get dt (PI ) = −ξ1 PI, where ξ1 > 0.

M. Haque / Nonlinear Analysis: Real World Applications 11 (2010) 2224–2236

2229

So PI → 0 as t → ∞. Again arguing as in the proof of Lemma 1(ii), P → 0 and H → K as t → ∞. As t → ∞, 1 dI I dt

1

→ β − (γ + b − r θ ) − β I = β 1 −

R30

! −I .

Given 1 > 0 there exists t1 such that for t ≥ t1 , 1 dI I dt and

1

≤β 1− 1 dI I dt

i.e. 1 −

1 R30

R30

+

≥β 1−

!

1

−I ,

2 1

−

R30

1 2

! −I ,

−  1 ≤ I∞ ≤ I ∞ ≤ 1 −

1

+ 1 .

R30

But 1 is arbitrary so I ∞ = I∞ = 1 −

1 R30

,

I →1−

1

as t → ∞.

R30

This completes the proof. Note that this proof works for R30 = 1 too. Hence E2 is also GAS in the case R20 < 1, R30 = 1 if H (0) > 0 and ξ1 > 0. (iii) A similar technique to [17], we verified that the system undergoes a transcritical bifurcation at a = a5 .  3.4. Dynamics of the system around E7 = (P ∗ , H ∗ , I ∗ ) Lemma 4. (i) The necessary and sufficient conditions for the system (2.4) around the positive interior equilibrium E7 are locally asymptotically stable (LAS) if the following inequalities hold true: (I)

 (II)

aP ∗

ecP ∗ I ∗ (q − 1) + cH ∗ I ∗ (q − 1) ≤ H

rθ K



− c + c (q − 1)H I ≤

(III) ecP q − 2 − 2(q − 1)I ∗

(IV) (1 − q) (1 − θ )

≤ βI

aP ∗

∗ ∗

!

arP ∗ MK

≤

arecI ∗ P ∗ 2 MK

M M

rH ∗ K

+

rH ∗ K

, !

+ β + ecP (q − 1) I ∗ , ∗

! + β + ecP ∗ (q − 1) I ∗ ,

+ ec P H (1 − (1 − q)I ) ecq + 2 ∗

∗

∗

r θ H ∗I ∗

!!

K

! + ec 2 P ∗ H ∗ (1 − I ∗ + qI ∗ )2 + rec 2 (q − 1)2 P ∗ I ∗ (1 − I ∗ )H ∗ .

Proof. (i) To show the asymptotic stability of the equilibrium E7 (P ∗ , H ∗ , I ∗ ), we use the method of first approximation. In this method we have to show that (i) The second compound matrix J [2] (E7 ) of J ∗ is stable and (ii) det J ∗ < 0. Now for E7 (P ∗ , H ∗ , I ∗ ) and the diagonal matrix D = diag (I ∗ , H ∗ , P ∗ ), the matrix J ∗ is similar to DJ ∗ (E7 )D−1 = (ϑij )3×3 , where ∗ ∗ ∗ ∗ ∗ ∗ ϑ11 = σ¯ 11 ,ϑ12 = σ¯ 12 HI ∗ , ϑ13 = σ¯ 13 PI ∗ , ϑ21 = σ¯ 21 HI ∗ , ϑ22 = σ¯ 22 , ϑ23 = σ¯ 23 HP ∗ , ϑ31 = σ¯ 31 PI ∗ , ϑ32 = σ¯ 32 HP ∗ , ϑ33 = σ¯ 33 , and J [2] = (σ¯ ij )3×3 , the 2nd computed matrix of J. The matrix J [2] (E7 ) is stable if and only if DJ [2] (E7 )D−1 is stable. Since the diagonal elements of the matrix DJ [7] (E7 )D−1 are negative, therefore the matrix will be stable if it is diagonally dominant in rows. Let µ∗ = max{Ψ1 , Ψ2 , Ψ3 }, where

! Ψ1 = ϑ11 + ϑ12 + ϑ13 =

ecP I (q − 1) + cH I (q − 1)

Ψ2 = ϑ21 + ϑ22 + ϑ23 = H

∗

∗ ∗



rθ K

∗ ∗

− c + c (q − 1)I

∗

! Ψ3 = ϑ31 + ϑ32 + ϑ33 = ecP q − 2 − 2(q − 1)I

−

!



∗

−

aP ∗

− rH ∗ K

M

aP ∗ M

+

rH ∗ K

! , ! !

+ β + ecP (q − 1) I ∗ , ∗

! ! + β + ecP (q − 1) I ∗ . ∗

2230

M. Haque / Nonlinear Analysis: Real World Applications 11 (2010) 2224–2236

When the above conditions hold true then µ∗ < 0, which implies the diagonal dominance as claimed and thus verifies the first condition. Again det J7 = H (1 − q) (1 − θ ) ∗

− H∗ βI

arP ∗ MK

arecI ∗ P ∗ 2 MK

+ ec P H (1 − (1 − q)I ) ecq + 2 ∗

∗

∗

r θ H ∗I ∗

!

!!

K

!

+ ec 2 P ∗ H ∗ (1 − I ∗ + qI ∗ )2 + rec 2 (q − 1)2 P ∗ I ∗ (1 − I ∗ )H ∗ .

Therefore if the above-mentioned conditions hold true then det J7 < 0, this verifies the second condition and completes the proof.  Remark 1. A set of suitable sufficient conditions for the local stability of the system around E7 can be obtained, but we did not include these results due to lack of space. 4. Behavior of the proportional system (2.4) and the absolute system (2.3) The behavior of the system when one or more of the susceptible predator, infected predator and prey populations are initially absent has already been discussed in Section 3. For starting values for which all three types of individuals are initially present note that for all parameter values there is exactly one LAS equilibrium. The following remarks can be done. Remark 2. For the proportional system (2.4) with P (0) > 0, H (0) > 0 and I (0) > 0 we have the following results: (i) For R20 < 1, R30 ≤ 1 and ξ1 > 0 the system approaches E2 , the healthy predator-only equilibrium. (ii) For a < a5 , ξ1 > 0 the system approaches E5 , the other predator-only equilibrium. Now we propose the following conjectures that will be verified numerically in Section 7. Conjecture 1. For R20 > 1 and R50 < 1 the system approaches E4 , the prey- susceptible predator equilibrium. Conjecture 2. For the conditions stated in Lemma 4, the system approaches E7 , the unique positive equilibrium where all three types of individuals are present. Therefore we see that R20 , R30 and R50 are the important parameters which determine the long-term dynamics of the system. R10 and R40 are not so important. The next question is what are the implications of our results for the original system given by (2.3). 4.1. Implications for the absolute system β XY

When we consider the absolute system (2.3) the natural interpretation is that the undefined term H is zero when H = 0. With this interpretation we find that there are six possible distinct equilibria: E 0 = (0, 0, 0, 0). The predators and prey are extinct. This is always admissible. E 1 = (M , 0, 0, 0). The prey only equilibrium. This is always admissible. E 2 = (0, K , K , 0). The prey and the infected predators are extinct. This is always admissible. E 4 = (P4 , H4 , H4 , 0). The infected predator is extinct. This is admissible if and only if R20 > 1. E 5 = (0, K , K (1 − I5 ), KI5 ). The prey is extinct. This is admissible if and only if R30 > 1. E 7 = (P ∗ , H ∗ , H ∗ (1 − I ∗ ), H ∗ I ∗ ). Both the susceptible and infected predator and the prey are present. This is possible if G2 + 4H 3 > 0 and 1 − I ∗ = ψ ∗ in (i)

(ii)

0, min 1, cqK − a

(q − 1)cK

a − cqK

!! if q < 1 with

(1 − q)cK , min 1,

(iii) (0, 1) if q = 1 with

aKr + ecaqMK eca(q − 1)MK a cqK

a cqK

> 1,

!! if q > 1 with

a cqK

< 1 and

> 1.

The equilibrium point E 0 corresponds to E0 and E3 ; E 1 to E1 and E6 in the proportional system. It is straightforward to show that the local and global stability properties of those equilibria E 2 , E 4 , E 5 and E 7 where H 6= 0 follow from those of the corresponding equilibria in the proportional model. However when we try to evaluate the Jacobian of (2.3) at E 0 and E 1 we have a problem as it contains undefined terms.

M. Haque / Nonlinear Analysis: Real World Applications 11 (2010) 2224–2236

2231

It is straightforward to show that global asymptotic stability of E2 or E5 in the proportional system implies that of E 2 and E 5 respectively. Instability does not follow similarly but by considering the second equation in (2.3) it is straightforward to show that E 0 is always unstable. When R20 > 1 by a similar technique it can be shown that E 2 is unstable. We make the following remarks about the absolute system Remark 3. For the absolute system (2.3) with P (0) > 0, H (0) > 0, and X (0) > 0, Y (0) > 0: (i) For R20 < 1, R30 < 1 and ξ1 > 0 the system approaches E 2 , the susceptible predator-only equilibrium. (ii) For R20 > 1 and R50 < 1 the system approaches E 4 , the prey-susceptible predator equilibrium. Similar to the proportional system (2.4), we eventually propose the following conjectures for the absolute system (2.3) Conjecture 3. For R30 > 1, a < a5 and ξ1 > 0, the system approaches E 5 , the predator only equilibrium. Conjecture 4. For the conditions stated in Lemma 4, then the system approaches E 7 , the unique equilibrium where all three types of individuals are present. 5. The role of infection on the predator–prey system To observe the influence of infection in the predator–prey system we set I = 0 in (2.4), and obtain a simple Lotka–Volterra type predator–prey system, with density-dependent logistic growth for the prey population, where the predator has an alternative food source. However it is analyzed by Bazykin et al. [22], but here we report the results in terms of our system parameters for comparison purpose. This system has four biologically relevant equilibria namely : (i) Eˆ 0 = (0, 0), (ii) rM (a−cK ) aK (r +ecM ) Eˆ 1 = (M , 0), (iii) Eˆ 2 = (0, K ) and (iv) Eˆ 4 = (P4 , H4 ) where P4 = ar +ec 2 MK and H4 = ar +ec 2 MK . The first three are always admissible, the last is admissible if and only if R20 > 1. It can easily be verified that Eˆ 0 is always unstable with eigenvalues a

and r, Eˆ 1 is unstable in the H-axis direction with eigenvalue r + ecM, Eˆ 2 will be GAS for H (0) > 0 if R20 < 1 and unstable in the P-axis direction if R20 > 1. It is easily verified that the positive interior equilibrium Eˆ 4 = (H4 , P4 ) is always LAS when it exists. Note that if R20 < 1 then Eˆ 2 will be GAS and Eˆ 4 does not exist. As R20 increases above one Eˆ 2 changes to being unstable and Eˆ 4 comes into existence. 6. The effect of prey removal on the eco-epidemiological system To observe the role of the predator in the eco-epidemiological system we put P = 0 in our system (2.4) and get a reduced system which is well known and has been analyzed extensively by Diekmann and Heesterbeek [23], Hethcote [24] and Brauer and Castillo-Chavez [25]. Here the local stability results have been reported in terms of our system parameters for comparison purposes. This system has four equilibria, namely (i) E˘ 0 = (0, 0) which is always admissible, (ii) E˘ 2 = (K , 0), which is always admissible, (iii) E˘ 3 = (0, I2 ), which is admissible if and only if R10 > 1, and (iv) E˘ 5 = (K , I5 ), which is admissible if and only if R30 > 1. It is easy to show that both the trivial equilibrium E˘ 0 and the axial equilibrium E˘ 3 are

unstable. The axial equilibrium point E˘ 2 = (K , 0) will be GAS for H (0) > 0 if R30 < 1 and unstable if R30 > 1. Arguing as in the proof of Lemma 3(ii) we deduce that if R30 > 1, H (0) > 0 and I (0) > 0 then E˘ 5 is GAS.

The absolute system (2.3) when P ≡ 0 has three distinct equilibria: (i) E˜ 0 = (0, 0, 0), (ii) E˜ 2 = (K , K , 0) and E˜ 5 = (K , K (1 − I5 ), KI5 ). The first two always exist. The third exists if and only if R30 > 1. E˜ 0 is always unstable. E˜ 2 is GAS for H (0) > 0 and X (0) > 0 if R30 < 1 and unstable if R30 > 1. If R30 > 1, H (0) > 0, X (0) > 0 and Y (0) > 0 then E˜ 5 is GAS. 7. Numerical investigations Many biological examples are given in the Table 2 of [19], here we take the example of rabbit/fox predator–prey system as in the UK, the European rabbit plays a vital role in maintaining bio-diversity and the fox (Vulpus vulpus) has become the main carrier of rabies in Europe. Now let the fox be infected by rabies. The disease is transmitted mostly as a result of an infected animal biting a healthy animal. The birth rate of foxes b, its natural death rate d, transmission coefficient of rabies β , per capita rate of recovery from the disease γ and density dependent parameter values θ are taken from the paper of [26] (see Table 3). They calculated these parameter values per day here we convert them per month and we obtain (b, d, β, γ , θ ) = (0.00274, 0.00137, 0.21833, 0, 0.00069) × 30 = (0.0822, 0.0411, 6.5499, 0, 0.0207). We assume that γ is very low and would be 0.00015/month however they calculated 0 per day. The birth rate of rabbit and fox vary in time and they are higher in the summer than in the winter. Here we take the growth rate of rabbits from the same source. Therefore we take the growth rate of rabbits (a) in between 0.2466 (in New Zealand) and 0.6 (in United Kingdom) per month. All other parameter values have been chosen in such a way that they are realistic and at the same time obey the conditions for stability or bifurcation.

2232

M. Haque / Nonlinear Analysis: Real World Applications 11 (2010) 2224–2236

80

60 I 40

20

0 120 250

100

200 150

80

100

H

[2]

Fig. 1. Global stability around E2 . Here R0 tend to E2 .

60

50 0

P

= 0.1233 < 1 and R0[3] = 0.9262 < 1; and all the trajectories starting from any part of the positive octant

1360

1745

0.06

1340 0.04

1740

1320

0.02

1300

1735

1280 P

0 H

1260

I

1730

–0.02

1240 1725

1220

–0.04

1200

1720

–0.06

1180 1160

0

5 beta

10

1715

0

5 beta

10

–0.08

0

5 beta

10

Fig. 2. Transcritical bifurcation around E4 at β = β [4] = 6.6502. The thick line represents the loci of stable equilibria, the dotted line the unstable ones.

To show the global stability of the system around E2 we choose c = 0.02, e = 0.4, q = .3, M = 1000, K = 100 and all other parameters are described before; Fig. 1 shows this fact under different initial conditions and in all cases we get that all the trajectories starting from any part of the positive octant tend to E2 . The yellow circle with a red border represents the corresponding equilibrium point in Figs. 1, 3, 5 and 6. By taking c = .0002, γ = 6.5, M = 3000, K = 500, we get the transcritical bifurcation around E4 when β = β [4] = 6.6502; Fig. 2 describes the fact. In order to show the global stability of the system around E5 (see Fig. 3) we choose same parameter values as for Fig. 1 except e = 0.2, q = 0.3, M = 3000 and K = 500. The system gets transcritical bifurcation with the same parameter values when it passes through a = a5 = 3.089 (see Fig. 4).

M. Haque / Nonlinear Analysis: Real World Applications 11 (2010) 2224–2236

2000

2233

600

1000

400

P

H 0

–1000

200

0

5

10 15 Time in month

0

20

0

5

10 15 Time in month

20

10 8 4 6

I

I 2

4 2 0

1000

2000

500 0

5

10 15 Time in month

20

1000 P

0 0

H

Fig. 3. Global stability around E5 . Here R0 R0 = 1.1023 > 1 ≥ R0 = 0.9684. [2] [4]

7

[2]

600

0.9876

580

0.9876

6 5 4 P

0.9876

560

3

I

H 540

2

0.9876 0.9876

1

520

0.9876

0 500

–1 –2

0

1

2

3

4

5

6

480

7

0.9875 0

a

1

2

3

4

5

6

7

0.9875

a

0

1

2

3

4

5

6

7

a

Fig. 4. Transcritical bifurcation around E5 at a = a5 = 3.089. The thick line represents the loci of stable equilibria, the dotted line the unstable ones.

Now we check our conjectures which have been made in this manuscript. Fig. 6 proofs the Conjecture 2. In the Conjecture 1, we claim that for R20 > 1 and R50 < 1, the system approaches E4 , Fig. 5 shows the fact with all the parameter values same as for Fig. 2 except β < β [4] . 8. Discussion In this article we have proposed and analyzed a mathematical model that consists of three non-linear autonomous differential equations for three different populations namely prey (P), susceptible predator (H), infected predator (I). The conditions for existence and stability of the equilibria of the system have been given. The bifurcation situations have also been observed around important equilibrium points. β a We have investigated that there are five epidemiological threshold quantities for our model: R10 = γ +b ; R20 = cK ; R30 = β ; γ +b −r θ

R40 =

H4 =

β+ecM (q−1) , b+γ

R50 =

β

γ +b+ec (1−q)P4 −

r θ H4 K

and provided that R20 ≥ 1, where

aK (r + ecM ) ar + ec 2 MK

is the equilibrium predator level in the presence of the prey when disease is absent in the prey. Note that R30 ≥ R10 .

2234

M. Haque / Nonlinear Analysis: Real World Applications 11 (2010) 2224–2236

0.65 0.6 0.55 I

0.5 0.45 0.4 0.35 1500 300

1400

280 260

1300

240 220

1200

H

P

Fig. 5. All the trajectories that arise from any part of the positive octant tend to E7 .

0.4

0.3 I 0.2

0.1

0 2000 1400

1800 1300 1600

1200 1400

H

1100

P

Fig. 6. System approaches to E4 with R0 = 6 > 1 and R0 = 0.9849 < 1. Note that all the trajectories starting from any part of the positive octant tend to E4 . [2]

[5]

To show the influence of transmissible disease in the simple predator–prey system we compare the system (2.3) with the system considered in Section 6. In the absence of infection, the only important threshold quantity is R20 which is the basic reproduction number for the susceptible predator invading the infection free equilibrium. If predator and prey are both initially present then we expect that the prey dies out if R20 < 1 and takes off if R20 > 1. Therefore the situation is simpler when disease is absent. Note that if the predator is initially present it does not go to extinction whether infection is present or absent. Ecologically this is because the alternative food source for the predator species saves the predator population from global extinction. Again note that the predator-only equilibrium is E 2 = (0, K , K , 0) in system (2.3) corresponding to Eˆ 2 = (0, K ) in system considered in Section 5. For H (0) > 0, Eˆ 2 is GAS if R20 < 1 and unstable if R20 > 1. On the other hand for H (0) > 0, and X (0) > 0 E 2 is GAS if R20 < 1 and R30 < 1 and unstable if either R20 > 1 or R30 > 1. Hence global asymptotic stability of E 2 in the presence of disease requires additional conditions. So if R20 < 1 but R30 > 1 then E 2 is unstable. In these circumstances the prey population goes to extinction if no infection is present, but does not go to extinction if infection is

M. Haque / Nonlinear Analysis: Real World Applications 11 (2010) 2224–2236

2235

initially present. Thus infection in the predator species may save the prey from extinction even if R20 , the basic reproduction number for the prey to be able to invade the predator-only equilibrium, is less than one. In nature there are many situations where extinction of the prey species is undesirable and the introduction of a non-fatal disease into the predator is a potential method of biological control to ensure persistence of the prey. On the other hand if the prey is a pest which we wish to eliminate then the criteria for global stability of the predator-only equilibria E 2 and Eˆ 2 tell us when this is admissible. Similarly the predator-prey only equilibrium E 4 = (P4 , H4 , H4 , 0) in the absolute full model (2.3) corresponds to the equilibrium Eˆ 4 = (P4 , H4 ) in the full model where infection is not present i.e. the model considered in Section 6. The conditions for both equilibria to exist are the same (R20 > 1) but Eˆ 4 is always LAS when it exists, whereas the condition for the existence of an equilibrium with prey, susceptible predator and infected predator all present (R50 > 1) which destabilizes E 4 in the full model. This conclusion agrees fairly with the conclusion of [8], where the authors have shown that a small amount of infectious agent can destabilize the otherwise stable trophic configuration between a prey species (phytoplankton) and the predator (its grazer) where a disease is spreading among the prey species only. The results of their model are compared to experimental data by use of a combination of nonlinear forecasting techniques. We would also like to point out that in another two eco-epidemiological models where the predator consumes the infected prey, Haque and Venturino [14], and Greenhalgh and Haque [16] have shown that if a similar threshold condition is satisfied one can destabilize an otherwise stable prey only-predator equilibrium where disease is spreading among the prey species only. Thus theoretical studies indicate that infection in either prey species can destabilize an otherwise stable prey-predator equilibrium. For the full model with absolute numbers we conjecture that if R20 > 1 and R50 < 1, so E 4 exists and is LAS then E 4 is actually GAS provided P (0) > 0, H (0) > 0 and X (0) > 0 so the disease dies out but the prey and susceptible predator populations persist. In Section 6, we observed that the basic reproduction number R30 is much more important in determining the behavior of the predator–prey system when the prey is removed than when it is present. In the full model system (2.3) the predatoronly equilibrium E 2 = (0, K , K , 0) is LAS if R20 < 1 and R30 < 1. On the other hand, the corresponding prey free equilibrium

E˜ 2 = (K , K , 0) is GAS if R30 < 1. Therefore, for R20 > 1 and R30 < 1, E 2 = (0, K , K , 0) become unstable and E˜ 2 = (K , K , 0) is GAS. Thus in the absence of prey, the disease will die out if H (0) > 0, and X (0) > 0 and R30 < 1, which does not happen when prey is present. Therefore, the absence of prey benefits the predator when a transmissible disease runs among the predator species. If we consider the total predator population size then provided that H (0) > 0 the predator increases to its carrying capacity, K , when the prey is absent. However the presence of the prey can increase the predator population by driving to an upper equilibrium level E 4 (P4 , H4 , X4 , 0) in the absence of disease. In conclusion, this article emphasizes that the value of the per capita disease contact rate β and the existence of the prey population have important implications for predator species persistence in simple prey-susceptible predator-infected predator models where the predators are infected by transmissible disease. Acknowledgments I am very grateful to the referees for their comments. I am greatly indebted to John R. King (editor of IMA Journal of Mathematics in Medicine and Biology, associate editor of SIAM Journal of Applied Mathematics), from whom I have learned applications of mathematics to biology and medicine. Finally, I am also thankful to Professor Philip Maini, Mathematical Institute(CMB), University of Oxford, who saw the earlier version of the manuscript. Finally, thanks go to the UK-IERI for financial support. Appendix The first equation of the system (2.4) implies that dP

≤ aP (1 − P /M ). dt Therefore by a simple but standard argument we have lim supt →∞ P (t ) ≤ M. So given  > 0 there exists t0 such that for t ≥ t0 , P (t ) ≤ M +  . Clearly lim supt →∞ I (t ) ≤ 1. Again the 2nd equation of (2.4) implies that: dH

≤ (r (1 − H /K ) + ec (1 + q)(M + ))H , for t ≥ t0 . dt Thus lim supt →∞ H (t ) ≤ K (1 + ec (1 + q)(M + )/r ). As  is arbitrary lim supt →∞ H (t ) ≤ K (1 + ec (1 + q)M /r ) = K1 (say). References [1] M. Haque, Ratio-dependent predator-prey models of interacting populations, Bulletin of Mathematical Biology 71 (2) (2009) 430–452. [2] W. Kermack, A. McKendrick, A contribution to the mathematical theory of epidemics, Proceedings of the Royal Society of London 115 (1927) 700–721. [3] D. Minchella, M. Scott, Parasitism: A cryptic determinant of animal community structure, Trends in Ecology & Evolution 6 (8) (1991) 250–254.

2236

M. Haque / Nonlinear Analysis: Real World Applications 11 (2010) 2224–2236

[4] C. Gibson, A. Watkinson, The role of the hemiparasitic annual Rhinanthus minor in determining grassland community structure, Oecologia 89 (1) (1992) 62–68. [5] A. Dobson, M. Crawley, Pathogens and the structure of plant communities, Trends in Ecology & Evolution 9 (10) (1994) 393–399. [6] K. Hadeler, H. Freedman, Predator–prey populations with parasitic infection, Journal of Mathematical Biology 27 (6) (1989) 609–631. [7] H. Freedman, A model of predator–prey dynamics as modified by the action of a parasite, Mathematical Biosciences 99 (2) (1990) 143–155. [8] E. Beltrami, T. Carroll, Modeling the role of viral disease in recurrent phytoplankton blooms, Journal of Mathematical Biology 32 (8) (1994) 857–863. [9] E. Venturino, The influence of diseases on Lotka–Volterra systems, Rocky Mountain Journal of Mathematics 24 (1) (1993) 381–402. [10] E. Venturino, The effects of diseases on competing species, Mathematical Biosciences 174 (2) (2001) 111–131. [11] E. Venturino, Epidemics in predator–prey models: Disease in the predators, Mathematical Medicine and Biology 19 (3) (2002) 185–205. [12] Y. Xiao, L. Chen, Modeling and analysis of a predator–prey model with disease in the prey, Mathematical Biosciences 171 (1) (2001) 59–82. [13] H. Hethcote, W. Wang, L. Han, Z. Ma, A predator–prey model with infected prey, Theoretical Population Biology 66 (3) (2004) 259–268. [14] M. Haque, E. Venturino, The role of transmissible diseases in the Holling–Tanner predator–prey model, Theoretical Population Biology 70 (3) (2006) 273–288. [15] M. Haque, E. Venturino, An ecoepidemiological model with disease in predator: The ratio-dependent case, Mathematical Methods in the Applied Sciences 30 (14) (2007). [16] D. Greenhalgh, M. Haque, A predator–prey model with disease in the prey species only, Mathematical Methods in the Applied Sciences 30 (8) (2007). [17] M. Haque, J. Zhen, E. Venturino, An ecoepidemiological predator–prey model with standard disease incidence, Mathematical Methods in the Applied Sciences 32 (7) (2009). [18] M. van Baalen, V. Krivan, P. van Rijn, M. Sabelis, Alternative food, switching predators, and the persistence of predator–prey systems, The American Naturalist 157 (5) (2001) 512–524. [19] F. Gulland, The impact of infectious diseases on wild animal populations: A review, Ecology of Infectious Diseases in Natural Populations 1 (1995) 20–51. [20] R. Anderson, R. May, Population biology of infectious diseases, Report of the Dahlem workshop on population biology of infectious disease agents, 1982. [21] H. Thieme, Epidemic and demographic interaction in the spread of potentially fatal diseases in growing populations, Mathematical Biosciences 111 (1) (1992) 99–130. [22] A. Bazykin, A. Khibnik, B. Krauskopf, Nonlinear Dynamics of Interacting Populations, World Scientific, 1998. [23] O. Diekmann, J. Heesterbeek, Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation, Wiley, 2000. [24] H. Hethcote, The mathematics of infectious diseases, SIAM Review 42 (4) (2000) 599–653. [25] F. Brauer, C. Castillo-Chavez, Mathematical Models in Population Biology and Epidemiology, Springer, 2001. [26] J. Ireland, R. Norman, J. Greenman, The effect of seasonal host birth rates on population dynamics: The importance of resonance, Journal of Theoretical Biology 231 (2) (2004) 229–238.