# The Dynamics of an Eco-Epidemiological Model with Nonlinear Incidence Rate

Authors: Raid Kamel Naji; Arkan N. Mustafa
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

This paper treats the dynamical behavior of eco-epidemiological model with nonlinear incidence rate. A Holling type II prey-predator model withSI-type of disease in prey has been proposed and analyzed. The existence, uniqueness, and boundedness of the solution of the system are studied. The local and global dynamical behaviors are investigated. The conditions, which guarantee the occurring of Hopf bifurcation of the system, are established. Finally, further investigations for the global dynamics of the proposed system are carried out with the help of numerical simulations.


---

## Body

## 1. Introduction

In the beginning of the twentieth century a number of attempts have been made to predict the evolution and existence of species mathematically. Indeed, the first major attempt in this direction was due to the well-known classical Lotka-Volterra model in 1927, since then many complicated model for two or more interacting species has been proposed according to the Lotka-Volterra model by taking into account the effect of competition, time delay, functional response, and so forth, see for example [1, 2] and the references theirin. On the other hand, over the last few decades, mathematics has been used to understand and predict the spread of disease relating important public-health questions to basic transmission parameters, for the detailed history of mathematical epidemiology and basics for SIR epidemic models (or Kermack-McKendrick model) may be found in the classical books [3–5]. However, recently Haque and Venturino [6] have been discussed mathematical models of diseases spreading in symbiotic communities.During the last three decades, there has been growing interest in the study of infectious disease coupled with prey-predator interaction model. In many ecological studies of prey-predator systems with disease, it is reported that the predators take a disproportionately high number of parasite-infected prey. Some studies have even shown that parasite could change the external features or behavior of the prey so that infected prey are more vulnerable to predation, see [7–13] and the references theirin.Later on, many authors have been proposed and studied eco-epidemiological mathematical models incorporating ratio-dependent functional response, toxicant, external sources of disease, predator switching, and infected prey refuge [14–23]. In all these models, the authors assumed that the infection affects the prey only and the disease transmission follows the simple law of mass action with a constant rate of transmission, also they assumed that the predator consumes either the susceptible prey or the infected prey. Recently, Haque and Greenhalgh in 2010 [24] have proposed and studied a predator-prey model with logistic growth in the prey population, where a disease spreads among the prey according to a susceptible-infected-susceptible (SIS) epidemic model. They assumed that the predator do not consume infected prey.On the other hand, there is another category papers in literature, in which the authors consider the eco-epidemiological models where the disease spreads in predator population [25–27]. The authors in [25] have proposed a ratio-dependent demographic predator-prey model in which a disease spreads among predators via homogeneous mixing, while the author in [26] proposed 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. Finally, in [27] the authors considered a system of delay differential equations modeling the predator-prey eco-epidemic dynamics with a transmissible disease in the predator population.Keeping the above in view, we combine a prey-predator model with an epidemiological model where the disease in prey is modeled by a susceptible-infected (SI) epidemic system. The eco-epidemiological prey-predator model proposed here differs from previous models; it uses the nonlinear incidence rate [28] and assume that the predator consumes the susceptible as well as the infected prey according to the modified Holling type II functional response. Finally, the formulation of this model and the local as well as the global stability analysis of the proposed model are described in the following sections.


## 2. The Mathematical Model

An eco-epidemiological model consisting of prey, which is divided into two classes: susceptible prey and infected prey, interacting with a predator is proposed according to the following assumptions.(1)
In the absence of disease, the prey population grows logistically with carrying capacityK>0 and an intrinsic growth rate constant r>0.(2)
In the presence of disease, the prey population is divided into two classes, namely the susceptible preyS(t) and infected prey I(t), and hence the total prey population at time t will be S(t)+I(t). Further, it is assumed that only the susceptible prey can reproduce reaching to its carrying capacity. However, the infected prey does not grow, recover, reproduce, or compete.(3)
It is assumed that the disease is transmitted from infected prey to susceptible prey by contact according to the nonlinear incidence rate of the formλIS/(1+I), which was proposed by Gumel and Moghadas 2003 [28] and used by many authors, where λI measures the infection force of the disease and 1/(1+I) measures the inhibition effect from the crowding effect of the infected individuals. This incidence rate seems more reasonable than that depends on simple law of mass action λSI because it includes the crowding effect of the infected individuals and prevents the unboundedness of the contact rate by choosing suitable parameters.(4)
Finally, it is assumed that the predator species consumes the prey species (susceptible as well as infected) according to the modified Holling type II functional response. However, in the absence of prey population the predator population decays exponentially.Then the dynamics of such a model can be represented in the following set of differential equations:(2.1)dSdt=rS(1-SK)-λI1+IS-α1SYβ+S+mI=f(S,I,Y);S(0)≥0,dIdt=λI1+IS-α2IYβ+S+mI-μ1I=g(S,I,Y);I(0)≥0,dYdt=θ1S+θ2Iβ+S+mIY-μ2Y=h(S,I,Y);Y(0)≥0.
Here, all the parameters are assumed to be positive constants. Moreover, the parameter λ represents the infected rate; α1 and α2 represent the maximum predation rates of S and I, respectively; β is the half saturation constant; the parameter m represents the predator preference rate between S and I. The parameters θ1 and θ2 are the conversion rates of S and I, respectively; μ1 represents the death rate if I due to the disease, while μ2 represents the natural death rate of Y. In addition, since the density of population cannot be negative then, the state space of the system (2.1) is R+3={(S,I,Y)∈R3:S≥0,I≥0,Y≥0}.Note that system (2.1) can be separated into two independent subsystems. The first system is obtained by assuming the absence of the predators and can be written in the following form:(2.2a)dSdt=rS(1-SK)-λI1+IS=f1(S,I),dIdt=λI1+IS-μ1I=g1(S,I).
However, the second subsystem is obtained in the absence of the infected prey and takes the form:
(2.2b)dSdt=rS(1-SK)-α1SYβ+S=f2(S,Y),dYdt=θ1Sβ+SY-μ2Y=h2(S,Y).Obviously, the interaction functionsf, g, and h of the system (2.1) are continuous and have continuous partial derivatives on the state space R+3; therefore these functions are Lipschitzian on R+3 and then the solution of the system (2.1) with non negative initial condition exists and is unique. In addition, all the solutions of the system (2.1) which initiate in nonnegative octant are uniformly bounded as shown in the following theorem.Theorem 2.1.
All solutions of system (2.1) that initiate in the state space R+3 are uniformly bounded in the region Ψ={(S,I,Y)∈R+3:0≤S+I+Y≤θ/μ}, where μ=min(μ1,μ2) and θ=(r+μ)K.Proof.
Let(S(t),I(t),Y(t)) be any solution of the system with the nonnegative initial conditions. From the first equation, we get
(2.3)dSdt≤rS(1-SK).
Then by solving the above differential inequality, we obtain
(2.4)Limt→∞SupS(t)≤K.
Let W=S+I+Y, then from the model we get
(2.5)dWdt=rS(1-SK)-(α1-θ1)SYβ+S+mI-(α2-θ2)SYβ+S+mI-μ1I-μ2Y.
Now, since the conversion rate constant from prey population to predator population cannot be exceeding the maximum predation rate constant of predator population to prey population. Hence from biological a point of view, we have always θi≤αi; i=1,2. Hence, we obtain that
(2.6)dWdt≤rS-μ1I-μ2Y≤θ-μW.
So again by solving the above linear differential inequality, we get that
(2.7)Limt→∞W(t)≤θμ.
So the proof is completed.Now since the dynamical system (2.1) is said to be dissipative if all populations initiating in R+3 are uniformly limited by their environment. Accordingly, system (2.1) is a dissipative system.


## 3. Stability Analysis of Subsystems

In this section, the local stability of the subsystems (2.2a) and (2.2b) is discussed and the obtained results are summarized below.The first subsystem (2.2a) has at most three nonnegative equilibrium points. The equilibrium points p1=(0,0) and p2=(K,0) always exist. However, the positive equilibrium point p3=(S^,I^) where
(3.1)S^=μ1λ(1+I^);I^=-A22A1+12A1A22-4A1A3
with A1=rμ1>0, A2=rμ1+λ2K-r(λK-μ1) and A3=-r(λK-μ1) exists in the interior of the positive quadrant of the SI-plane if and only if
(3.2)λK>μ1,
otherwise p3 does not exist. Further, it is observed that since
(3.3)dIdt=I(λ1+IS-μ1)≤I(λK-μ1).
So if λK<μ1, then the right-hand side will be negative and then I species become extinct. Therefore, condition (3.2) represents the necessary condition for survival of infected species.Moreover, the variational matrices of subsystem (2.2a) about the equilibrium points p1 and p2 can be written as
(3.4)V(p1)=(r00-μ1),V(p2)=(-r-Kλ0Kλ-μ1).
Thus, the eigenvalues of V(p1) are  λ11=r>0 and λ21=-μ1<0, then p1 is a saddle point with locally stable manifold in the I-direction and with locally unstable manifold in the S-direction. While the eigenvalues of V(p2) are λ12=-r<0 and λ22=λK-μ1. So, p2 is locally asymptotically stable point if and only if
(3.5)λK<μ1,
and it is a saddle point with locally stable manifold in the S-direction and with locally unstable manifold in the I-direction provided that condition (3.2) holds.In addition to the above, the variational matrix about the positive equilibrium pointp3 can be written as follows V(p3)=(aij)2×2, where
(3.6)a11=-rS^K<0,a12=-λS^(1+I^)2<0,a21=λI^1+I^>0,a22=-λS^I^(1+I^)2<0.
Obviously, the characteristic equation of V(p3) can be written in the following form:
(3.7)λ2+Aλ+B=0
with A=-(a11+a22)=(rS^)/K+(λS^I^)/(1+I^)2>0(3.8)B=a11a22-a12a21=rλS^2I^K(1+I^)2+λ2S^I^(1+I^)3>0,
Consequently, the eigenvalues of V(p3) are given by
(3.9)λ13=-A2+12A2-4B;λ23=-A2-12A2-4B.
Clearly both eigenvalues λ13 and λ23 have negative real parts, and hence p3=(S^,I^) is locally asymptotically stable in the interior of SI-plane whenever it exists.Similarly, the second subsystem (2.2b) has three nonnegative equilibrium points q1=(0,0), q2=(K,0), and q3=(S-,Y-), where
(3.10)S-=βμ2θ1-μ2;Y-=rα1K(K-S-)(β+S-)=rβθ1α1K(θ1-μ2)2[K(θ1-μ2)-βμ2].
Clearly, q3 exists in the interior of positive quadrant of SY-plane under the following conditions:
(3.11)θ1>μ2,S-<K.
Now the variational matrices of the subsystem (2.2b) at q1 and q2 can be written as
(3.12)V(q1)=(r00-μ2),V(q2)=(-r-λKβ+K0(θ1-μ2)K-βμ2β+K).
So, the eigenvalues of V(q1) are γ11=r>0 and γ21=-μ2<0; hence, q1 is a saddle point with locally stable manifold in Y-direction and with locally unstable manifold in the S-direction.Also the eigenvalues ofV(q2) are γ12=-r<0 and γ22=(K(θ1-μ2)-βμ2)/(β+K). So, q2 is locally asymptotically stable if and only if
(3.13)θ1<μ2orK<βμ2θ1-μ2.
However, it is a saddle point with locally stable manifold in the S-direction and with locally unstable manifold in the Y-direction provided that condition (3.11) holds.Finally, the variational matrix of the subsystem (2.2b) at the positive equilibrium point q3 can be written in the form V(q3)=(bij)2×2, where
(3.14)b11=rKS-[-1+K-S-β+S-],b12=-α1S-β+S-<0,b21=βθ1Y-(β+S-)2>0,b22=0.
Therefore, it is easy to verify that the eigenvalues of V(q3) satisfy the following relations:
(3.15)γ13+γ23=rKS-[-1+K-S-β+S-];γ13γ23=βα1θ1Y-S-(β+S-)3>0.
Thus, if the following condition holds: then both eigenvalues have negative real parts and hence q3 is locally asymptotically stable in the interior of the positive quadrant of SY-plane:(3.16a)[-1+K-S-β+S-]<0.
Further, both eigenvalues have positive real parts and hence q3 will be unstable point if the following condition holds:
(3.16b)[-1+K-S-β+S-]>0.


## 4. Stability Analysis of System (2.1)

In this section, system (2.1) has been analyzed locally as well as globally. An investigation of system (2.1) shows that there are at most five possible nonnegative equilibrium points. The existence conditions for these points are discussed in the following.The trivial equilibrium pointE0=(0,0,0) and the axial equilibrium  point E1=(K,0,0) always exist.The disease-free equilibrium pointE2=(S-,0,Y-), where S- and Y- are given by (3.10), exists under condition (3.11).The predator-free equilibrium pointE3=(S^,I^,0), where S^ and I^ are given by (3.1), exists under condition (3.2).The positive-equilibrium pointE4=(S~,I~,Y~) exists in the Int.R+3 if and only if S~, I~, and Y~ represent the positive solution of the following algebraic set of nonlinear equations: (4.1a)rα1(1-SK)-λIα1(1+I)=Yβ+S+mI,(4.1b)λα2(1+I)S-μ1α2=Yβ+S+mI,(4.1c)θ1S+θ2Iβ+S+mI-μ2=0.
Now, from (4.1c) we get
(4.1d)S=A~1+B~1IC~1
with A~1=μ2β>0; B~1=mμ2-θ2 and C~1=θ1-μ2. So, by combining (4.1a) with (4.1b) and then substituting (4.1d) in the resulting equation we get that
(4.1e)D~1I2+D~2I+D~3=0,
where D~1=rα2B~1, D~2=rα2A~1+(rα2+λα1K)B~1+(λα2-rα2-α1μ1)KC~1, and D~3=(rα2+λα1K)A~1-(rα2+α1μ1)KC~1.Obviously, (4.1e) has a unique positive root, say I~, if and only if D~1 and D~3 has opposite signs and hence S~=S(I~) follows directly from (4.1d). Finally, substituting the value of (S~,I~) in (4.1a) gives
(4.1f)Y~=(β+S~+mI~)[r(K-S~)(1+I~)-λKI~α1K(1+I~)].Clearly, Y~ is positive if and only if the following condition holds:(4.2a)r>λKI~(K-S~)(1+I~).
Therefore, the positive equilibrium point E4=(S~,I~,Y~) exists uniquely in the Int.R+3 if and only if in addition to condition (4.2a) one set of the following sets of conditions is satisfied.
(4.2b)C~1>0,B~1>0,D~3<0(4.2c)C~1>0,B~1<0,D~3>0I~<βμ2θ2-mμ2(4.2d)C~1<0,B~1<0,I~>βμ2θ2-mμ2.Now in order to study the local stability of system (2.1), the variational matrix of system (2.1) is computed at each of the above equilibrium points and then the eigenvalues are determined as shown in the following.The Variational matrix at the trivial equilibrium point is determined as(4.3)V(E0)=(r000-μ1000-μ2).
So, the eigenvalues of V(E0) are λ1=r>0, λ2=-μ1<0 and λ3=-μ2<0; hence, E0=(0,0,0) is a saddle point with locally stable manifold in IY-plane and with locally unstable manifold in the S-direction.The variational matrix at the axial equilibrium point can be written as(4.4)V(E1)=(-r-λK-α1Kβ+K0λK-μ1000(θ1-μ2)K-βμ2β+K).
Hence, the eigenvalues of V(E1) are λ1=-r<0, λ2=λK-μ1, and λ3=((θ1-μ2)K-μ2β)/(β+K). Therefore, E1=(K,0,0) is locally asymptotically stable if and only if (4.5a)λK<μ1withθ1<μ2orK<min{μ1λ,μ2βθ1-μ2}.
However, it is a saddle point provided that
(4.5b)λK>μ1and/orK>μ2βθ1-μ2>0.Now, the variational matrix of system (2.1) at the disease-free equilibrium point E2=(S-,0,Y-) is written as
(4.6)V(E2)=(rKS-[-1+K-S-β+S-]S-[-λ+mr(K-S-)K(β+S-)]-α1S-β+S-0λS--α2Y-β+S--μ10θ1βY-(β+S-)2Y-[θ2β+S-(θ2-mθ1)(β+S-)2]0).
Then, the eigenvalues of V(E2) satisfy the following relations λ-S+λ-Y=γ13+γ23 and λ-Sλ-Y=γ13·γ23, where γ13 and γ23 represent the eigenvalues of V(q3) and satisfy (3.15). However,
(4.7)λ-I=λS--α2Y-β+S--μ1,
here λ-S, λ-I, and λ-Y represent the eigenvalues of V(E2) in the S-, I-, and Y-directions, respectively. Therefore, it is clear that the eigenvalues λ-S and λ-Y have negative real parts if and only if condition (3.16a) holds. However the eigenvalue λ-I is negative or positive if and only if the following conditions hold respectively:(4.8a)S-<K(α1μ1+α2r)α1Kλ+α2r,(4.8b)S->K(α1μ1+α2r)α1Kλ+α2r.Consequently, E2=(S-,0,Y-) is locally asymptotically stable if and only if conditions (3.16a) and (4.8a) hold. However it is a saddle point with non empty stable and unstable manifolds if at least one of conditions (3.16b) and (4.8b) hold.The variational matrix of system (2.1) at the predator-free equilibrium point E3=(S^,I^,0) can be written as
(4.9)V(E3)=(-rKS^-λS^(1+I^)2-α1S^β+S^+mI^λI^1+I^-λS^I^(1+I^)2-α2I^β+S^+mI^00S^(θ1-μ2)+I^(θ2-mμ2)-βμ2β+S^+mI^).
Straightforward computation shows that, the eigenvalues of V(E3) can be written as λ^S=λ13; λ^I=λ23, where λ13 and λ23 are given by (3.9), and the eigenvalue in the Y-direction is written as
(4.10)λ^Y=S^(θ1-μ2)+I^(θ2-mμ2)-βμ2β+S^+mI^.
According to (3.9), both eigenvalues λ^S and λ^I have negative real parts, while the eigenvalue λ^Y will be negative or positive if and only if the following conditions hold respectively:(4.11a)S^(θ1-μ2)+I^(θ2-mμ2)<μ2β,(4.11b)S^(θ1-μ2)+I^(θ2-mμ2)>μ2β.Therefore, E3=(S^,I^,0) is locally asymptotically stable provided that condition (4.11a) holds. However, it will be a saddle point with locally stable manifold in SI-plane and with locally unstable manifold in the Y-direction provided that condition (4.11b) holds.Finally, the local stability conditions for the positive equilibrium point are established in the following theorem.Theorem 4.1.
The positive equilibrium pointE4=(S~,I~,Y~) is locally asymptotically stable in the Int.R+3 provided that the following conditions hold: (4.12a)Y~<min{rα1KH~2,λS~mα2G~2H~2,λmα1G~2H~2},(4.12b)mθ1S~β+S~<θ2<mθ1orθ2I~β+mI~<θ1<θ2m,(4.12c)Γ>KG~2(λH~2-mα1G~2Y~rH~2-α1KY~),(4.12d)λH~2S~-mα2G~2Y~G~(λH~2+α2G~Y~)>α1S~α2I~,where  G~=1+I~,  H~=β+S~+mI~, and Γ=(βθ2+(θ2-mθ1)S~)/(βθ1+(mθ1-θ2)I~).Proof.
Straightforward computation shows that the variational matrix nearE4=(S~,I~,Y~) can be written as V(E4)=(cij)3×3 with i,j=1,2,3 where
(4.13)c11=-rS~K+α1S~Y~H~2,c12=-λS~G~2+α1mY~S~H~2,c13=-α1S~H~<0,c21=λI~G~+α2I~Y~H~2>0,c22=-λS~I~G~2+mα2I~Y~H~2,c23=-α2I~H~<0,c31=Y~βθ1+Y~I~(mθ1-θ2)H~2,c32=Y~βθ2+Y~S~(θ2-mθ1)H~2,c33=0.
Then, the characteristic equation of V(E4) can be written as follows:
(4.14)λ~3+D1λ~2+D2λ~+D3=0,
where
(4.15)D1=-(c11+c22),D2=(c11c22-c12c21)-(c13c31+c23c32),D3=c23(c11c32-c12c31)+c13(c22c31-c21c32),Δ=D1D2-D3=D1(c11c22-c21c12)+c31(c11c13+c12c23)+c32(c22c23+c21c13).
Note that according to Routh-Hurwitz criterion E4 is locally asymptotically stable if D1>0,  D3>0, and D1D2-D3>0.
Clearly from the condition (4.12a) we obtain that c11<0, c22<0 with c12>0, and hence D1>0. Also due to condition (4.12b), we obtain that c31>0 and c32>0. Therefore, the second term of D3 will be positive provided that conditions (4.12a) and (4.12b) hold. Further, it is easy to verify that the first term of D3 will be positive and hence D3>0 if the (sufficient) condition (4.12c) holds.
Now, according to the above set of conditions we have that the first and second terms ofΔ are positive. However the third term of Δ will be positive too if the following (sufficient) condition holds c22c23+c13c21>0, which is satisfied provided that condition (4.12d) holds. This completes the proof.Furthermore, in the following theorem the global stability conditions ofE1 are established.Theorem 4.2.
Assume thatE1=(K,0,0) is locally asymptotically stable in R+3. Then, if the following condition holds
(4.16)α1βK<μ2,
the equilibrium point E1 is globally asymptotically stable.Proof.
Consider the following function:(4.17)V(S,I,Y)=[S-K-Kln(SK)]+I+Y,
it is easy to see that V(S,I,Y)∈C1(R+3,R), in addition V(K,0,0)=0, while V(S,I,Y)>0 for all (S,I,Y)∈R+3 and (S,I,Y)≠(K,0,0). Further,
(4.18)dVdt=(S-K)[r(1-SK)-λI1+I-α1Yβ+S+mI]+[λS1+II-α2YIβ+S+mI-μ1I]+[θ1S+θ2Iβ+S+mIY-μ2Y].
Now due to the fact θi≤αi; i=1,2 that is mentioned in Theorem 2.1, then we get
(4.19)dVdt≤-rK(S-K)2-[μ1-λK1+I]I-[μ2-Kα1β+S+mI]Y.
Clearly, μ1-λK/(1+I)>0 under the local stability condition (4.5a). However, μ2-Kα1/(β+S+mI)>0 under the condition (4.16). Therefore, dV/dt is negative definite, and hence the proof is complete.Theorem 4.3.
Suppose that the disease-free equilibrium pointE2=(S-,0,Y-) is locally asymptotically stable in the Int.R+2 of the SY-plane and let the following condition holds:
(4.20)α1β2<rKYmax.
Then, the equilibrium point E2=(S-,0,Y-) is globally asymptotically stable in the Int.R+2 of the SY-plane.Proof.
Obviously, system (2.1) will be reduced to disease-free subsystem (2.2b) in case of absence of the infected species (I=0). Also, it has been shown that (S-,Y-) represents the unique positive equilibrium point of the subsystem (2.2b).
Consider now the functionF2(S,Y)=1/SY, clearly F2:Int.R+2→R which is a continuously differentiable function. Further, since
(4.21)Δ1(S,Y)=∂∂S(F2f2)+∂∂Y(F2h2)=-rKY+α1(β+S)2<-rKYmax+α1β2.
Here, f2 and h2 are given in subsystem (2.2b). Also, since the system (2.1) is bounded as shown in Theorem 2.1, Ymax represents the upper bound constant for the variable Y. Therefore, under the condition (4.20) it is observed that Δ1(S,Y) does not change sign and is not identically zero. So, by Bendixson-Dulac criterion, there is no periodic curve in the Int.R+2 of the SY-plane. Hence, the equilibrium point (S-,Y-) of subsystem (2.2b) and then the associated disease-free equilibrium point (S-,0,Y-) of system (2.1) is globally asymptotically in the Int.R+2 of the SI-plane.Theorem 4.4.
Assume that the positive equilibrium pointE4=(S~,I~,Y~) of system (2.1) is locally asymptotically stable in the Int.R+3. Then, it is a globally asymptotically stable on the subregion of the Int.R+3 that satisfies the following conditions:(4.22a)Y~<min{rHH~α1K,λS~HH~mGG~}(4.22b)[(C1mα1+C2α2)Y~1HH~+(C2G~-C1)λGG~]2<4C1C2[rK-α1Y~HH~][λS~GG~-mα2Y~HH~],where G=(1+I) and H=(β+S+mI), while G~ and H~ are given in Theorem 4.2. Further, C1 and C2 are positive constants to be determined later on.Proof.
Consider the following function:(4.23)V1(S,I,Y)=C1[S-S~-S~ln(SS~)]+C2[I-I~-I~ln(II~)]+[Y-Y~-Y~ln(YY~)],
where C1 and C2 are positive constants to be determined. Note that it is easy to see that V1(S,I,Y)∈C1(R+3,R) with V1(S~,I~,Y~)=0 and V1(S,I,Y)>0, for all (S,I,Y)∈R+3 with (S,I,Y)≠(S~,I~,Y~). Further,
(4.24)dV1dt=C1(S-S~)[r(1-SK)-λIG-α1YH]+C2(I-I~)[λSG-α2YH-μ1]+(Y-Y~)[Q1S+Q2IH-μ2]dV1dt=-C1[rK-α1Y~HH~](S-S~)2-C2[λS~GG~-mα2Y~HH~](I-I~)2+[(C1mα1+C2α2)Y~1HH~+(C2G~-C1)λGG~](S-S~)(I-I~)+[-C1α1H~+(βθ1-(θ2-mθ1)I~)HH~](S-S~)(Y-Y~)+[-C2α2H~+(βθ2+(θ2-mθ1)S~)HH~](I-I~)(Y-Y~).
So by choosing the constants Ci; i=1,2 as follows:
(4.25)C1=βθ1-(θ2-mθ1)I~α1H~,C2=βθ2-(θ2-mθ1)S~α2H~.
Then, we get
(4.26)dV1dt=-C1[rK-α1Y~HH~](S-S~)2-C2[λS~GG~-mα2Y~HH~](I-I~)2+[(C1mα1+C2α2)Y~1HH~+(C2G~-C1)λGG~](S-S~)(I-I~).
Clearly, Ci; i=1,2 are positive under the local stability condition (4.12b). Further, r/K-α1Y~/HH~>0 and λS~/GG~-mα2Y~/HH~>0 under the condition (4.22a). Therefore using condition (4.22b), we obtain that
(4.27)dV1dt<-[C1[rK-α1Y~HH~](S-S~)+C2[λS~GG~-mα2Y~HH~](I-I~)]2.
Therefore, dV1/dt is negative and hence V1 is Lyapunov function with respect to E4, so E4 is globally asymptotically stable and the proof is complete.


## 5. Hopf Bifurcation

In this section, the occurrence of Hopf bifurcation in system (2.1) near the equilibrium points is studied.Theorem 5.1.
Assume that condition (4.8a) holds, then system (2.1) has a Hopf bifurcation near the disease-free equilibrium point E2=(S-,0,Y-) as the parameter value K passes through the value K-=β+2S-.Proof.
According to the variational matrix of system (2.1) at the disease-free equilibrium point E2=(S-,0,Y-), it is easy to verify that the eigenvalues can be written as
(5.1)λ-S,Y=T2±12T2-4D,λ-I=λS--α2Y-β+S--μ1,
where T=(r/K)S-[-1+(K-S-)/(β+S-)] and D=α1θ1Y-βS-/(β+S-)3>0.
Clearly, as shown above, the eigenvalueλ-I<0 if and only if condition (4.8a) holds. However, the eigenvalues λ-S and λ-Y are pure imaginary complex numbers for T=0 or K=K-, so there is a neighborhood around K-=β+2S- such that λ-S and λ-Y can be written as
(5.2)λ-S=ω(K)+iϖ(K),λ-Y=ω(K)-iϖ(K),
where  ω(K)=(r/2K)S-[-1+(K-S-)/(β+S-)] represents the real part of λ-S and λ-Y.
Now, since(5.3)[dω(K)dK]K=K-=rS-2(β+2S-)2+rS-22(β+2S-)2(β+S-)≠0.
Therefore, system (2.1) has a Hopf bifurcation near the disease-free equilibrium point at K=K-, and hence the proof is complete.Now since the variational matrix around the predator-free equilibrium point has two eigenvalues, given by (3.9), with negative real parts while the third that given by (4.10) is real and is negative or positive depending on conditions (4.11a) or (4.11b), respectively, then there is no possibility to have a Hopf bifurcation near this point.Finally, the conditions that guarantee the occurring of Hopf bifurcation near the positive equilibrium point are established in the following theorem.Theorem 5.2.
Assume that the conditions (4.12a)–(4.12c) hold. Then, system (2.1) exhibits a Hopf bifurcation near the positive equilibrium point E4 as the parameter K passes through the value
(5.4)K1=-2rc22S~H~2H~2[N1+N12+4c22N2]-2α1c22S~Y~,
where N1=c13c31+c12c21-c222, N2=c12(c21c22+c23c31)+c32(c22c23+c21c13), and cij;i,j=1,2,3 represent the elements of V(E4): provided that the following condition holds
(5.5)[λ2S~2I~2G~4+α1α2m2S~I~2Y~2H~4]Π1+α2λS~I~Y~(1+I~)G~2H~3Π2<(α2+α1S~)mλS~I~2Y~G~2H~2Π1+α1α2mS~I~Y~2H~5Π2+I~Y~H~3[α2Y~(mI~+α1S~)H~2+α1λS~G~]Π3.
Here, Π1=λ/G~+(α2Y~)/H~2, Π2=βθ1+(mθ1-θ2)I~, and Π3=βθ2+(θ2-mθ1)S~.Proof.
Recall that the characteristic equation ofV(E4) is given by(5.6a)λ~3+D1λ~2+D2λ~+D3=0,
where Dj; j=1,2,3 are given in Theorem 4.1. It had been observed that the conditions (4.12a)–(4.12c) guarantee that Dj>0; j=1,2,3 for all values of K.
Now since the Hopf bifurcation near the positive equilibrium pointE4 of system (2.1) occurs if and only if V(E4) have two complex conjugate eigenvalues with the third eigen value real and negative such that there exists a constant parameter value, say K1, satisfying:(1)Re(λ~(K))|K=K1=0.(2)(d/dK)(Re(λ~(K)))|K=K1≠0; where λ~ is a complex eigenvalue of V(E4).
So, by simplifyingΔ=D1D2-D3 and then equating to the zero we get c22(c11)2-N1c11-N2=0.
Obviously, conditions (4.12a)-(4.12b) guarantee that N1<0, while condition (5.5) guarantees that N2<0. Therefore, by solving the above second order equation we get the c11=N1/2c22+(1/2c22)N12+4c22N2 (The other root is omitted because c11<0).
Substituting the value ofc11 in this equation and then solving for K we get that K=K1. Accordingly, for K=K1, we have D1D2=D3 and then the above characteristic equation (5.6a) can be written as
(5.6b)(λ~2+D2)(λ~+D1)=0,which has the following three roots:
(5.7)λ~1=iD2(K1),λ~2=-iD2(K1),λ~3=-D1(K1).
However, for all values of K in the neighborhood of K1, these roots can be written in general as
(5.8)λ~1(K)=a(K)+ib(K),λ~2(K)=a(K)-ib(K),λ~3(K)=-D1(K).
Clearly, a(K1)=0 which means that the first condition of Hopf bifurcation holds. Now, the proof will follows, if we can verify the above second condition (known as transversality condition) of Hopf bifurcation when Re(λ~(K))=a(K).
Thus, by substitutingλ~1(K)=a(K)+ib(K) in (5.6b) and calculating the derivative with respect to the K, it is obtained that
(5.9)A(K)a'(K)-B(K)b'(K)+E(K)=0,B(K)a'(K)+A(K)b'(K)+R(K)=0,
where
(5.10)A(K)=3[a(K)]2+2D1(K)a(K)+D2(K)-3[b(K)]2,B(K)=6a(K)b(K)+2D1(K)b(K),E(K)=[a(K)]2D1′(K)+D2′(K)a(K)+D3′(K)-D1′(K)[b(K)]2,R(K)=2a(K)b(K)D1′(K)+D2′(K)b(K).
Thus, by solving the linear system (5.9) for the unknown a'(K) we get that
(5.11)a′(K)=ddKRe(λ~(K))=-B(K)R(K)+A(K)E(K)[A(K)]2+[B(K)]2.
So for K=K1, it is easy to verify that
(5.12)ddKRe(λ~(K))|K=K1=-B(K1)R(K1)+A(K1)E(K1)[A(K1)]2+[B(K1)]2≠0.
Also, we have that λ~3(K1)=-D1(K1)<0. Hence the proof is complete.


## 6. Numerical Simulation

To visualize the above analytical findings and understand the effect of varying the parameters on the global dynamics of system (2.1), numerical simulation is done in this section. The system (2.1) is integrated using six-order Runge-Kutta method along with predictor-corrector method under different sets of initial values and different sets of parameter values.For the following set of hypothetical, biologically feasible set of parameters, system (2.1) is solved numerically starting at different initial points as shown in the following figures:
(6.1)r=1,K=500,λ=1,α1=1,α2=1,β=50,m=0.75,θ1=0.75,θ2=0.5,μ1=0.3,μ2=0.1.
According to Figures 1 and 2, system (2.1) has a unique globally asymptotically stable positive equilibrium point in the Int.R+3. For the parameters values (6.1) with the infective rates λ=0.5, λ=0.01, and λ=1.25, system (2.1) is solved numerically and the trajectories are drawn in Figures 3, 4, and 5 respectively.Figure 1
System (2.1) approaches globally asymptotically stable point under the parameters values (6.1).Time series for the point attractor in Figure1. (a) The trajectories starting at (10,10,10). (b) The trajectories starting at (1,1,1). (c) The trajectories starting at (20,20,20).
(a)
(b)
(c)(a) System (2.1) approaches periodic attractor in Int.R+3 for λ=0.5 with the rest of parameters are given by (6.1). (b) Time series of the trajectory in (a).
(a)
(b)(a) System (2.1) approaches periodic attractor in Int.R+2 of the SY-plane for λ=0.01 with the rest of parameters are given by (6.1). (b) Time series of the trajectory in (a).
(a)
(b)(a) System (2.1) approaches stable point in Int.R+2 of the SI-plane for λ=1.25 with the rest of parameters are given by (6.1). (b) Time series of the trajectory in (a).
(a)
(b)Clearly, Figure3 shows that system (2.1) losses its stability and approaches to periodic attractor as the infective rate decreases to λ=0.5, keeping the rest of parameters as in (6.1). While the disease will disappear and system (2.1) approaches periodic attractor in the Int.R+2 of the SY-plane as the infective rate decreases further reaching λ=0.01 as shown in Figure 4. Moreover, the predator will faces extinction as the infective rate increases reaching to λ=1.25, and hence the system will approach asymptotically stable point in the Int.R+2 of the SI-plane. In fact, it is observed that for 0<λ≤0.02 the trajectories of system (2.1) approach periodic attractor in the Int.R+2 of the SY-plane and for 0.02<λ<0.62 they approach the periodic attractor in the Int.R+3. However, increasing the infective rate farther, 0.62≤λ<1.12, stabilizes the system (2.1) at positive equilibrium point. Finally, for 1.12≤λ, the system approaches asymptotically stable point in the Int.R+2 of the SI-plane.Now, the effect of varying the intrinsic growth rate of the susceptible prey on the global dynamics of system (2.1) is investigated. It is observed that for r≤0.89 system (2.1) approaches asymptotically stable point in the Int.R+2 of the SI-plane; however, for 0.89<r<1.55, system (2.1) has a globally stable point in the Int.R+3. Finally, for 1.55≤r system (2.1) has a periodic attractor in the Int.R+3 as shown in Figure 6.(a) System (2.1) approaches periodic attractor in Int.R+3 for r=1.55 with the rest of parameters are given by (6.1). (b) Time series of the trajectory in (a).
(a)
(b)Now, for the following set of parameters values, it is observed that the system (2.1) has also periodic attractor in the Int.R+3. 						(6.2)r=1,K=500,λ=0.5,α1=0.5,α2=1,β=50,m=0.75,θ1=0.25,θ2=0.5,μ1=0.3,μ2=0.1.
However, the dynamics of system (2.1) is studied further at different values of the carrying capacity. It is observed that the system (2.1) undergoes a Hopf bifurcation as the carrying capacity parameter passes through K≅270, which confirm our result in Theorem 5.2, as shown in Figure 7.System (2.1) undergoes a Hopf bifurcation as the carrying capacity pass through K≅270 keeping other parameters fixed as in (6.2).
(a)
(b)
(c)
(d)
(e)
(f)Clearly, System (2.1) approaches asymptotically to the positive equilibrium point in the Int.R+3 when K=260 with the rest of parameters given by (6.2) as shown in Figures 7(a)-7(b). However, system (2.1) approaches to small periodic attractor in the Int.R+3 when K=280 with the rest of parameters given by (6.2) as shown in Figures 7(c) and 7(d). Finally, system (2.1) approaches larger periodic attractor in the Int.R+3 when K=300 with the rest of parameters given by (6.2) as plotted in Figures 7(e) and 7(f).In addition to the above, the effect of other parameters on the dynamics of system (2.1) is also studied and the obtained results can be summarized as follows.For the parameters values given by (6.2) with β≤160, system (2.1) has a periodic attractor, while as the half-saturation parameter β increases, β>160, it is observed that system (2.1) approaches asymptotically positive equilibrium point. On the other hand, for the parameters values in (6.2) with the natural death rate for the predator in the range μ2<0.3, system (2.1) has a periodic dynamic. However, for 0.3≤μ2≤0.43, the solution of system (2.1) approaches asymptotically stable positive equilibrium point. Further increase the natural death rate, that is, 0.43<μ2 will lead to extinction in the predator species and then the solution of system (2.1) approaches asymptotically the predator-free equilibrium point in the Int.R+2 of the SI-plane.


## 7. Discussions and Conclusions

In this paper, an eco-epidemiological prey-predator system with nonlinear incidence rate has been proposed and analyzed. The existence, uniqueness, and boundedness of the solution of the system are discussed. The existence and stability analysis of the equilibria of the system are investigated. The global stability analysis of the proposed model is also studied using suitable Lyapunov function. The conditions that guarantee the accruing of the Hopf bifurcation in system (2.1) are presented. Moreover, numerical simulations are used to study the global dynamics of our system and confirm our analytical results. It is observed the following: (1)
The infective rate constantλ has a vital role on the dynamics of system (2.1). In fact decreasing the value of λ more than a specific value, say λ1, causes destabilizing the system first and then leads to extinction in the infected individuals; however increasing λ more than a specific value, say λ2, leads to extinction in the predator species. Finally, for λ1<λ<λ2, system (2.1) has asymptotically stable point.(2)
Decreasing the intrinsic growth rate of the susceptible prey below a specific value, sayr1, leads to extinction in the predator species; however, increasing this parameter more than a specific value, say r2, causes coexistence of all species and the system approaches periodic dynamic. Finally for the values of intrinsic growth rate in the range r1<r<r2, the system approaches asymptotically stable point. Consequently, increasing this parameter has a destabilizing effect on the system.(3)
System (2.1) persists and has asymptotically stable point for relatively small values of carrying capacity; however, increasing the value of carrying capacity more than a specific value, say K>K1, the positive equilibrium point losses its stability and then the system approaches periodic attractor in the Int.R+3. Indeed, K1 represents a Hopf bifurcation point of system (2.1).Finally, the effects of other parameters of system (2.1) are also studied and similar conclusions are obtained as those mentioned in Section 6.



---

*Source: 852631-2012-10-11.xml*