Stability Analysis of an Eco-epidemiological Model Incorporating a Prey Refuge

The present paper deals with the problem of a predator-prey model incorporating a prey refuge with disease in the prey-population. We assume the predator population will prefer only infected population for their diet as those are more vulnerable. Dynamical behaviours such as boundedness, permanence, local and global stabilities are addressed. We have also studied the effect of discrete time delay on the model. The length of delay preserving the stability is also estimated. Computer simulations are carried out to illustrate our analytical findings.


Introduction
The dynamic relationship between predator and their prey has long been and will continue to be one of the dominant topics in both applied mathematics and theoretical ecology due to its universal existence and importance.These problems may appear to be simple mathematically at first sight, they are, in fact, often very challenging and complicated.
The formation of classical Lotka-Volterra [1,2] model was a milestone progress in the study of predator-prey interactions.Similarly, after the pioneering work of , epidemiological studies also received special attention to the researchers.Disease in ecological systems is an important issue.Anderson and May [4] were the pioneers for formulating the mathematical model on this topic.A lot of research articles have already been appeared on this subject [4,5].Most of the studies mainly focused on parasite infection in prey population only [6][7][8][9].The dynamics of predator-prey system with infection in prey population is an important study from modelling point of view.Most of predators preferentially consume diseased prey [10].Predators harm parasite directly by consuming infected prey and also harms parasites indirectly by reducing the density of susceptible hosts.Predators behaviour can magnify this direct effect because predators often prefer infected prey over uninfected prey [11].So that infectious disease can be a factor to regulate human and animal population size.
The research of the hiding behaviour of preys has been incorporated as a new ingredient of prey-predator models and its consequences on the dynamics of prey-predator interactions can be recognized as one of the major issues in both applied mathematics and theoretical ecology.In nature, prey populations often access to areas where they are safe from their predators.Such refugia are usually playing two significant roles, serving both to reduce the chance of extinction due to predation and to damp preypredator oscillations.These are therefore a potentially important means of increasing species richness in natural communities and of stabilizing population sizes, biomass and productivity.It is well known that many more attentions have paid on the effects of a prey refuge for predator-prey system.Predator-prey interactions often exhibit spatial refuge which afford the prey some degree of protection from predation and reduce the chance of extinction due to predation [3,12].Hassel [13] showed that adding a large refuge to a model, which exhibited divergent oscillations in the absence of refuge, replaced the oscillatory behaviour with a stable equilibrium.These mathematical models and a number of experiments indicate that refuge have a stabilizing effect on predator-prey interactions.
Time delays of one type or another have been incorporated into epidemiological models by many authors [1,2,14].In general, delay-differential equations exhibit much more complicated dynamics than ordinary differential equations since a time-delay could cause a stable equilibrium to become unstable and cause the populations to fluctuate.
In this paper, we have investigated the dynamical behaviour of a ratio-dependent predator-prey systems with infection in prey population [9], and the effect of refuge in the infected prey.Here we have studied the boundedness, permanence, local and global stabilities of the non-equilibrium points of this system.We have also considered a discrete time-delay in the interaction term of the predator equation.
The rest of the paper is structured as follows: In Section 2, we present a brief sketch of the construction of the model, which may indicate the epidemiological relevance of it.In Section 3, boundedness of the basic deterministic model (3) is discussed.Section 4 deals with the boundary equilibrium points and their stability.In Section 5, we find out the condition for which system (3) is permanent.In Section 6, we find the necessary and sufficient condition for the existence of the interior equilibrium point E * (s * , i * , y * ) and study its stability.Computer simulations of some solutions of the system (3) are presented in Section 7. The effect of discrete time-delay on the system (3) is studied in Section 8.In Section 9, we calculate the length of delay for which the system preserves stability.Section 10 contains the general discussions of the paper.
They make the following assumptions in formulating the mathematical model of a predator-prey system with disease in the prey population: 1.In the absence of disease, the prey population grows logistically with carrying capacity K ∈ + and intrinsic birth rate r 1 ∈ + .
2. In the presence of virus, the prey population is divided into two groups, namely susceptible prey denoted by S(T ) and infected prey denoted by I(T ).Therefore at time T , the total population is N (T ) = S(T ) + I(T ).
3. The disease is not genetically inherited.The infected populations do not recover or become immune.We assume that the disease transmission follows the simple law of mass action βS(T )I(T ) with β as the transmission rate.
4. The infected prey I(T ) is removed by death (say, its death rate is positive constant c) or by predation before having the possibility of reproducing.However, the infected prey population I(T ) still contribute with S(T ) towards the carrying capacity of the system.
5. The infected prey is more vulnerable than susceptible prey.We assume that the predator population consumes only infected prey with ratio-dependent Michaelis-Menten functional response function It is assumed that the predator has the death rate constant d (d > 0), and the predation coefficient b (b > 0).The coefficient in conversing prey into predator is p (0 < p ≤ 1).
This paper extends the above model by incorporating a refuge protecting mI of the infected prey, where m ∈ [0, 1) is constant.This leaves (1 − m)I of the infected prey available to the predator, and modifying system (1) accordingly yields the system: with initial data S(0) ≥ 0, I(0) ≥ 0, Y (0) ≥ 0.
The model we have just specified has nine parameters, which makes the analysis difficult.To reduce the number of parameters and to determine which combinations of parameters control the behaviour of the system, we nondimensionalize system (2).We choose , and t = βKT.
Then system (2) takes the form (after some simplification) 3 Boundedness In theoretical eco-epidemiology, boundedness of a system implies that the system is biologically well behaved.The following theorem ensures the boundedness of system (3): Theorem 1.All solutions of system (3) that start in 3 + are uniformly bounded.
Proof.Let, (s(t), i(t), y(t)) be any solution of the system (3).Since, ds dt ≤ rs(1 − s).We have, Applying a theorem on differential inequalities [15], we obtain and for t → ∞ Thus, all the solutions of (3) enter into the region Hence the theorem.

Boundary equilibria and their stability
In this section, we study the stability of the boundary equilibrium points of the system (3).
In the following lemma we have mentioned the boundary equilibria of the system (3) and the condition of their existence.
Lemma 1. System (3) always has two boundary equilibrium points, namely the trivial equilibrium E 0 (0, 0, 0) and the axial equilibrium E 1 (1, 0, 0).The predator-free equilibrium point E 2 (ŝ, î, 0) exists if and only if b 2 < 1.When this condition is satisfied, ŝ, î are given by In terms of original parameters of the system, the condition b 2 < 1 becomes c < βK.This implies that if the ratio of the death rate of the infected prey to the carrying capacity (c/K) is less than the transmissiom rate (β), then the predator becomes extinct and conversely.
The system (3) cannot be linearized at E 0 (0, 0, 0) and E 1 (1, 0, 0) and therefore local stability of E 0 and E 1 cannot be studied [16].Therefore, we are only interested in the stability in the predator-free equilibrium point E 2 (ŝ, î, 0).
The variational matrix V (E 2 ) at the equilibrium point E 2 is given by The eigen values are Since, B > 0 and C > 0, therefore, the signs of the real parts of λ 1 , λ 2 are negative.This implies that E 2 is locally asymptotically stable in the si-plane.Now E 2 is asymptotically stable in the y-direction if and only if pl − b 1 < 0, i.e., b 1 > pl.

Permanence of the system
To prove the permanence of the system (3), we shall use the "Average Liapunov" functions [17].
. Proof.We consider the average Liapunov function of the form V (s, i, y) = s α1 i α2 y α3 where each α i (i − 1, 2, 3) is assumed positive.
In the interior of 3 + , we have To prove the permanence of the system we shall have to show that ψ(s, i, y) > 0 for all boundary equilibria of the system.The following condition should be satisfied for the equilibrium point E 2 , After some simple calculation, we can get Hence the theorem.
6 The interior equilibrium point: its existence and stability First we consider the existence and uniqueness of the interior equilibrium point E * (s * , i * , y * ).
Lemma 2. The interior equilibrium point E * (s * , i * , y * ) of the system (3) exists if and only if the following two conditions are satisfied: Furthermore, s * , i * , y * are given by In terms of original parameters of the system, the conditions (i) and (ii) respectively become pb > d and ap(βK −c) > (1−m)(pb−d), which are the necessary and sufficient conditions for the co-existence of the susceptible prey, infected prey and the predator.
From Lemma 2, we can observe that the interior equilibrium point E * (s * , i * , y * ) exists if and only if both the conditions (i) and (ii) are satisfied.If any one of the condition is violated then E * (s * , i * , y * ) does not exists.Now from condition (ii) of the Lemma 2, we have, .
Hence, to exists the interior equilibrium point E * (s * , i * , y * ) the refuge constant m should lies in the interval Remark 1.It is to be noted that the existence of E * destabilizes E 2 .

Local stability analysis of E *
The variational matrix of ( 3) at E * is given by where The characteristic equation is where Theorem 3. E * is locally asymptotically stable if and only if D > 0 and ∆ > 0.
Proof.We notice that (i) D > 0 ⇔ A 1 > 0, (ii) A 3 > 0 for all values of the parameters, and Hence, the theorem follows from Routh Hurwitz criterion.

Global stability analysis of E *
Now, we shall study the global dynamics of the system (3) around the positive equilibrium E * (s * , i * , y * ).We use Liapunov function to prove the global result.
Theorem 4. Existence of positive interior equilibrium of the system of equations (3) implies its global stability around the positive interior equilibrium if the following two conditions: (ii) y > max{y * , iy * i * } or y < min{y * , iy * i * } holds true.Proof.Let us consider the following positive definite function about E * : Differentiating J with respect to t along the solution of (3), we have (after some calculation) Therefore, dJ dt is negative definite if the above conditions of the theorem are satisfied and consequently, J is a Liapunov function with respect to all solutions in the interior of the positive orthant, proving the theorem.

Numerical simulation
Analytical studies can never be completed without numerical verification of the results.In this section we present computer simulation of some solutions of the system (3).Beside verification of our analytical findings, these numerical solutions are very important from practical point of view.

Model with discrete delay
It is already mentioned that time-delay is an important factor in biological system.It is also reasonable to assume that the effect of the infected prey on the predator population will not be instantaneous, but mediated by some discrete time lag τ required for incubation.As a starting point of this section, we consider the following generalization of the model ( 3) involving discrete time-delay: All parameters are the same as in system (3) except that the positive constant τ represents the reaction time or gestation period of the predator y.
The system (5) has the same equilibria as in the previous case.The main purpose of this section is to study the stability behaviour of E * (s * , i * , y * ) in the presence of discrete delay (τ = 0).We linearize system (5) by using the following transformation: Then linear system is given by where and We look for the solution of the model ( 6) of the form U (t) = ρe λt , 0 = ρ ∈ .This leads to the following characteristic equation: where It is well known that the signs of the real parts of the solutions of ( 7) characterize the stability behaviour of E * .Therefore, substituting λ = ξ + iη in (7) we obtain real and imaginary parts, respectively as and A necessary condition for a stability change of E * is that the characteristic equation ( 7) should have purely imaginary solutions.Hence to obtain the stability criterion, we set ξ = 0 in ( 8) and ( 9).Then we have, Eliminating τ by squaring and adding ( 10) and ( 11), we get the equation for determining η as where Substituting η 2 = σ in (12), we get a cubic equation given by Since d 3 = a 2 3 − a 2 6 > 0 for the parameter values given in previous case, we assume that d 3 ≥ 0 and have the following claim.
then equation ( 13) has no positive real roots.
In fact, notice that Set, Then the roots of equation ( 15) can be expressed as If Hence neither σ 1 nor σ 2 is positive.Thus equation (8.12) does not have positive roots.Since h(0) = d 3 ≥ 0, it follows that the equation (8.9) has no positive roots.
Then the equilibrium point E * of the delay model (5) is absolutely stable; that is E * is asymptotically stable for all τ ≥ 0.
Remark.Proposition 1 indicates that if the parameters satisfy the conditions (i) and (ii), then the steady state of the delay model ( 5) is asymptotically stable for all delay values; that is, independent of the delay.However, we should point out that if the conditions (condition (ii)) in Proposition 1 are not satisfied, then the stability of the steady state depends on the delay value and the delay could even induce oscillation.
Proposition 2. Suppose that If either (ii) d 3 < 0 or (iii) d 3 ≥ 0 and d 2 < 0 is satisfied, then the steady state E * of the delay model ( 5) is asymptotically stable when 0 ≤ τ < τ 0 and unstable when τ > τ 0 , where when τ = τ 0 , a Hopf bifurcation occurs; that is, a family of periodic solutions bifurcates from E * as τ passes through the critical value τ 0 .
Proposition 2 indicates that the delay model could exhibit Hopf bifurcation at certain value of the delay if the parameters satisfy the conditions in (ii) and (iii).However, for the parameter values given in Section 7, neither (ii) nor (iii) holds.

Estimation of the length of delay to preserve stability
We consider the system (3) and the space of all real valued continuous functions defined on [−τ, ∞] satisfying the initial conditions on [−τ, 0].We linearize the system (3) about its interior equilibrium E * (s * , i * , y * ) and get Taking Laplace transform of the system given by ( 18), we get, and s 1 (α), i 1 (α) and y 1 (α) are the laplace transform of s 1 (t), i 1 (t) and y 1 (t) respectively.
We have already shown that E * (s * , i * , y * ) is stable in absence of delay.Hence, by continuity, all eigenvalues will continue to have negative real parts for sufficiently small τ > 0 provided one can gurantee that no eigenvalues with positive real parts bifurcates from infinity as τ increases from zero.This can be proved by using Butler's lemma [19].