A Predator-prey Model with Disease in Prey

The present investigation deals with the disease in the prey population having significant role in curbing the dynamical behaviour of the system of prey-predator interactions from both ecological and mathematical point of view. The predator-prey model introduced by Cosner et al. [1] has been wisely modified in the present work based on the biological point of considerations. Here one introduces the disease which may spread among the prey species only. Following the formulation of the model, all the equilibria are systematically analyzed and the existence of a Hopf bifurcation at the interior equilibrium has been duly carried out through their graphical representations with appropriate discussion in order to validate the applicability of the system under consideration.


Introduction
Around 1800, the British Economist Malthus formulated a single species model [2] and subsequently modified by Verhulst [3].Lotka and Volterra [4,5], considered two populations and analysed the model.They formulated the model after considering predatorprey type of situations.Many researchers have studied the techniques as predator-prey interactions, mutualisms and competitive mechanisms and made an attempt to develop a more biologically realistic model.Prey dependent predator-prey models have also been studied extensively in several investigations (cf.[6][7][8][9][10][11][12]). The deterministic prey dependent predator-prey model exhibits not only the "paradox of enrichment", formulated by Hariston et al. [13] and Rosenzweig [14] but also the "biological control paradox" which was taken by Luck [15].
Mathematical epidemiology has become a interesting subject of research work since the seminal model of  on SIRS (susceptible-infected-removedsusceptible) systems, in which a disease is transmitted upon contact has been thoroughly investigated.More recently epidemic models with demographic evolution have been introduced.Ecoepidemiology is much young subject in the branch of science.Here epidemic and demographic aspects are merged within one model.Gao and Hethcote [17] and Mena-Lorca and Hethcote [18] consider the dynamics of a reproducing population, which is also a subject of epidemics.A disease spreading in interacting populations was first studied by Hadeler and Freedman [19].The Lotka-Volterra model has been used as the demographic basis in which the influence of a disease propagating in one of the two species has been investigated.It is known that viral, bacterial, fungal and metazoan parasites can intervene host vulnerable to predation (cf.[20][21][22][23][24]).
The simplest models contain a bilinear mass action term, quadratic in both the interacting populations, called also Holling type I.This term appears due to the fact that an individual can in principle interact with the whole other population, the product of the two populations is the obvious outcome.We consider the fact that in general a single individual can feed only until the stomach is full, a saturation function indicate the intake of food.The latter can be modeled using the concept of the "law of diminishing returns" or technically speaking Michaelis-Menten or Holling type II term.The present model is a modification of the previous model studied by Cosner et al. [1], allowing a disease to spread among the prey species only.Based upon thorough analysis of the problem under study, all the equilibria are adequately characterized and their nature of stability are properly discussed.The standard approach has finally been used to establish the local stability, Hopf bifurcations and limit cycles.

Basic assumption and the mathematical model
A general type predator-prey model will have the following structure: dX dt = f (X)X − g(X, P )P, dP dt = eg(X, P )P − dP, where X, P are the population size of prey and predator species, respectively, f (X) the per capita prey growth rate in the absence of predators and g(X, P ) the rate at which an individual predator consumes prey.Also g(X, P ) represents the functional response of the model.The parameter e is the conversion factor and d is the natural death rate of predator population.It has been observed in [1] that if the predator searches prey with line formation and moving in a direction perpendicular or transverse to the line then the rate at which an individual predator consumes prey has been used by g(X, P ) = (c 1 X)/(a+X) as in the case of traditional Holling type II prey dependent response function.Such scenario assumed that encounters involved only an individual predator and a single prey item and that while one predator was handling prey, others would continue their searching strategy.Furthermore, the prey form patches, herds that are large enough that the predators can aggregate before all the prey are consumed or escape.This type of behaviour is limited by the requirement that the line of foragers must be short enough to permit transmission of a signal.So the corresponding response function is probably only accurate at low to moderate group size of predator.In this situation the number of encounters between prey and individual predator is proportional to the number of predators due to aggregation upon contact with prey and thus each encounter between the prey and a single predator is converted rapidly into an encounter between the prey and all the predators.As a result the total number of encounters depends upon another factor of P .Also suppose that the prey are gathered in large cluster that the predator do not entirely consume with the action of an encounter due to disperseness within some short period after the cluster is attacked.Then the predation rate per predator per encounter depends on that time and also on the handling time rather than depending on P .This is why we choose the better predator dependent response function as g(X, P ) = c 1 XP/(a+XP ).Now we deal with a predator-prey model with Cosner et al. type functional response [1] given by with X(0) = X 0 > 0 and P (0 where r, k, c 1 , a, e, d all belong to R + and r, k, c 1 and a represent growth rate of the prey, carrying capacity, search rate for prey and half saturation constant, respectively.The first equation characterizes that the prey population grows logistically with carrying capacity k and intrinsic growth rate r in the absence of predator population.System of equations ( 1) has the following equilibria: (i) trivial equilibrium Ê0 (0, 0), (ii) axial equilibrium Ê1 (k, 0) and (iii) positive interior equilibrium Ê * (X * , P * ), where X * = ad/((ec 1 − d)P * ) and P * is the real positive root of the cubic equation Z 3 + 3hZ + g = 0, where h = −rea/(3(ec 1 − d)) and g = a 2 der/(k(ec 1 − d) 2 ).This equation has exactly one real positive root if g 2 + 4h 3 > 0 i.e., 27ad 2 > 4k 2 er(ec 1 − d).Using Cardano's method, we obtain the root as m − h/m where m denotes one of the three values of [(1/2)(−g + (g 2 + 4h 3 )) 1/3 ].
The predators consume prey population according to the Cosner et al. [1] type of functional response.This type of functional response is different from others who have derived the fact that it increases with P .Logically, the functional response could reflect a higher rate of predation of prey per predator than would be possible if predator acted for foraging individually.It is very important to point out that similar "foraging" techniques were stated by German and American submarines in attacks on enemy convoys during second world war (cf.[25] and [26]).Now the above model is modified by introducing transmissible disease in the prey species only for the purpose of curbing the dynamical behaviour of the system.We introduce the following facts: (i) In the presence of disease, the prey population X is divided among susceptibles S and infected I individuals.Therefore, the total prey population is X(t) = S(t)+I(t).
(ii) Assume that only susceptible prey can compete for resources.
(iii) The disease spreads with bilinear mass action incidence rate λSI.
(iv) The infected prey population neither recover from the disease nor is capable of reproducing (cf.[27]).
(v) Further, we assume that the epidemics cannot be transmitted to the predator population.Predators hunt both sound and sick prey with different search rates depending on the various kind of parasitism (cf.[28][29][30]).
With the above modification we propose and analyze the following model: where r, k, λ, c 2 , c 1 , a, γ, e, d all belong to R + and represent growth rate of the prey, carrying capacity, force of infection, search rate for infected prey, search rate for susceptible prey, half saturation constant, total death rate of infected prey (natural death rate + death rate due to infection), conversion factor and the natural death rate of the predator population, respectively.The right hand side of each of equations ( 2)-( 4) is a smooth function of variables (S, I, P ) and all the parameters are non-negative.So local existence and uniqueness properties hold in the positive octant Π = {(S, I, P ): S > 0, I > 0, P > 0}.
From equation (2) it follows that S = 0 is an invariant subset that is S ≡ 0 for some t.Thus S(t) > 0 for all t, if S 0 > 0. The same argument follows for equations (3) and (4).Hence Π is an invariant set.If J be the Jacobian matrix of the system of equations (2)-(4) which is diagonalizable and hence system (2)-( 4) is obviously not conservative.
3 Boundedness of the system Proposition 1.The sound prey population is bounded.
Proof.From equation (2) we have By using simple but standard argument we have lim sup t→+∞ S(t) < k.Hence the proof.Now there exists a T 1 > 0 such that for all t > T 1 , we have S(t) < k + = W (say).
This characteristic is related to a similar threshold phenomenon in epidemic theory (cf.[32][33][34]).Now if the "basic reproductive ratio R 0 " in epidemic theory exceeds the value 1 then the epidemics outbreak will occur.This number, R 0 can be thought of as the number of all individuals who will get the disease in all time following "successful" contact with a typical sick individual, that is, the expected number of secondary cases produced by an infective during its whole infectious period.

Stability analysis
The system of equations ( 2)-( 4) can be written as in the form Ẋ = F (X), where X = (S, I, P ) T .The Jacobian matrix J ≡ DF (X) of system (2)-(4) at any arbitrary point www.mii.lt/NA(S, I, P ) is given by J(S, I, P ) = (δ ij ) 3×3 , where, assuming ζ = P/(a + (S + I)P ), the elements are defined by We denote J k = J, the Jacobian evaluated at E k , where k = 1, 2, 3, 4, 5.The determinant of J k = det(J), the trace of J k = tr(J) and , where M k denotes the sum of the second order principal minors of J k .

Local stability analysis for the simplest equilibria
It is easy to verify that the eigenvalues of the Jacobian matrix at E 0 are r, −γ, −d.Therefore, it is an unstable hyperbolic saddle critical point in the direction orthogonal to IP -plane.The eigenvalues of the Jacobian matrix at E 1 are −r, λk − γ and −d.Therefore, E 1 will be stable if the condition λk < γ holds.When λk < γ, then the system E 1 will be saddle with SP -plane its stable manifold and I-axis as unstable manifold.
Remark.Since the boundary equilibrium point E 2 is locally asymptotically stable for all positive values of the parameter of system ( 2)-( 4), therefore, the system around the positive interior equilibrium E 5 (S 5 , I 5 , P 5 ) is not persistent.
7 Non-existence periodic solutions around E 5 (S 5 , I 5 , P 5 ) In this section, we would like to prove that under some suitable conditions, there is no periodic solution of system ( 2)-( 4) around the positive interior equilibrium E 5 (S 5 , I 5 , P 5 ).
To prove this, the following criteria by Li and Muldowney [35] and Arino et al. [36] can be applied.Consider the general autonomous ordinary differential equation where F is a C 1 function in some open subset of R N .Denoting by J = (∂F/∂X) the Jacobian matrix of (14) and J [2] the N 2 × N 2 matrix which is the second additive compound matrix associated the Jacobian matrix J.The definition of the second additive compound matrix can be established in the paper of Li and Muldowney [35].Let J = (ξ ij ) Nonlinear Anal.Model.Control, 2013, Vol. 18, No. 2, 191-209 be an n × n matrix.The second additive compound matrix A [2] is the N 2 × N 2 matrix defined as follows, for any integer i = 1, 2, . . .N  2 , let (i) = (i 1 , i 2 ) be the ith number in the lexicographic ordering of integer pairs (i 1 , i 2 ) such that 1 i 1 < i 2 n.Then, the element in the ith row and jth column of J [2] is ) r+s ξ irjs if exactly one entry i r of (i) does not occur in (j) and j s does not occur in (i), 0 if neither entry from (i) occurs in (j).
A simple closed rectifiable curve that is invariant with respect to system (14) cannot exist if any of the following conditions is satisfied on R n : where λ 1 λ 2 • • • λ n are the eigenvalues of (1/2)((∂F/∂x) * +∂F/∂x) and where ∂F/∂x is the Jacobian matrix of F while the asterisk denotes transposition.
If X ∈ R N then the corresponding logarithmic norm of J [2] (that we denote by µ ∞ (J [2] )) endowed by the vector norm where µ ∞ (J [2] ) < 0 implies the diagonal dominance by row matrix J [2] .Then, the following result holds.
Before we arrive at a condition in which there is no closed rectified curve (no periodic solution), we carry out some changes of coordinates to lower the number the parameters in system ( 2)-( 4).Let us make the following parameter non-dimensional: System ( 2)-( 4) preserves the same form but with r = 1 and k = 1.The new parameters with prime ( ) mentioned in (15) are introduced in the foregoing analysis by dropping the prime notations.
Li and Muldowney's criteria has been adopted here in the revised coordinates for the non existence of periodic solutions of system ( 2)-( 4).The logarithmic norm µ ∞ , endowed by the vector norm |X| ∞ of the second additive compound matrix J [2] associated with the Jacobian matrix J, calculated on E 5 , is negative iff the suprema of the following functions satisfy: Sufficient conditions to satisfy ( 16), ( 17) and ( 18) are, respectively, (i) A direct application of Li and Muldowney's criteria shows that under the conditions (i), (ii) and (iii) there is no periodic solution for the present system (2)-( 4), under consideration.

Simulations
We have performed numerical simulation for the positive equilibrium of the updated system ( 2)-( 4) under consideration.The local stability characteristics of the present system around the equilibria E 2 , E 3 and E 5 are shown in Figs.1(a), 1(b) and 1(c).It is very interesting that the dynamical system enters into Hopf bifurcation at the interior equilibrium E 5 .Here we have investigated the phase portrait of Hopf bifurcation exhibited in Fig. 2 when the parameters have values r = 0.2 per month, k = 11000 metric tons, a = 1500 metric tons, c 1 = 0.00003 per month, c 2 = 0.00001 per month, λ = 0.0066 per month, parameter.This feature corresponds to the case when c 1 > c 2 , that is, when the consumption rate of susceptible prey population is higher than that of infected prey population.On the other hand when c 1 < c 2 , the parameters have values r = 0.8 per month, k = 50 metric tons, a = 1 metric tons, c 1 = 0.03 per month, c 2 = 1.5 per month, λ = 0.12 per month, e = 0.4, γ = 0.91 per month, d = 0.163 per month, then the system also enters into another Hopf bifurcation at the interior equilibrium E 5 .This happens when the consumption rate of infected prey population is greater than the rate of the susceptible prey population.In this context, one may observe the existence of a limit cycle.

Discussion
The dynamical behaviour investigated in our present model is an important issue of ecoepidemiological interactions based upon predator-prey model with disease in prey species only.Here we have considered three nonlinear autonomous ordinary differential eqations for three different classes of populations, namely, the susceptible prey, the infected prey and the predator.In this investigation, the boundedness of the solutions, the existence and stability of different equilibria have been thoroughly examined.The present system yields mainly five equilibria E 0 , E 1 , E 2 , E 3 and E 5 .We have summarized the sufficient conditions for the stability of all possible equilibria of model system ( 2)-( 4) in Table 1.It has been pointed out in Proposition 2 that the infection will spread only when the "basic reproductive ratio R 0 " is greater than 1.On the other hand when R 0 < 1 then the disease will naturally die out.The axial equilibrium position E 1 is found to be stable when λk < γ and unstable when γ < λk − r − d.The boundary equilibrium E 2 be always stable results in the system being non-persistent.The boundary equilibrium E 3 is locally asymptotically stable if the conditions (i), (ii) and (iii) of (7) in Proposition 3(I) hold and on the other hand when λ > λ [1] then E 3 is unstable.Also using the Routh-Hurwitz criteria it can be shown that Ê0 (0, 0) is always unstable, Ê1 (k, 0) is always stable and Ê * (X * , P * ) is stable if (i) k < 2X * and (ii) ec 1 X * P * × (2a + X * P * )/(a + X * P * ) 2 < d.The dynamical behaviour of the systems represented  by ( 1) and ( 2)-( 4) about Ê0 and E 0 , respectively, appears to be the same.Here Ê1 is always stable but E 1 is stable when λk < γ.Thus one may conclude that the dynamical behaviour of system (2)-( 4) is influenced by the force of infection.But the system is unstable when γ < (λk − r − d).Thus inclusion of the disease in the disease-free predator-prey system is adequate to destabilize the otherwise stable equilibrium.Again the equilibrium position E 3 is analogous with Ê * .But their dynamical behavior differs in several ways.Both the equilibria have the same feasibility condition, but for the stability of E 3 an additional condition λS 3 < γ is required.Again the system E 3 is unstable when λ > λ [2] .Thus the infection has an important influence on the ecosystem.Now one may conclude that the conditions (I) of the Proposition 3 are true then the inclusion of disease to the disease-free predator-prey system can help preventing total extinction and behave as a biological control.Next, further attention was focused on the interior equilibrium E 5 into two cases: E * 5 with condition c 1 > c 2 and E 5 * * with condition c 1 < c 2 in Propositions 4 and 5, respectively.In first cases we have noticed locally asymptotically stable behaviour in Proposition 4. In this situation the predator population wants to predict more susceptible prey population than infected prey population.On the other hand we have observed the condition when predators try to predict susceptible prey but infected prey is caught easily which is shown in the Proposition 5.Both the cases as mentioned above are usually found in our real life on the interior equilibrium E 5 into two cases: E * 5 with condition c 1 > c 2 and E * * 5 with condition c 1 < c 2 in Propositions 4 and 5, respectively.In first cases we have noticed locally asymptotically stable behaviour in Proposition 4. In this situation the predator population wants to capture more susceptible prey population than infected prey population.On the other hand we have observed the condition when predators try to hunt susceptible prey but infected prey is caught easily which is shown in the Proposition 5.Both the cases as mentioned above are usually found in real life.Lastly, we describe on both cases when the system enters Hopf bifurcation with bifurcation parameter λ.The stable limit cycles justify a behaviour similar to the one exhibited by the demographic model [37].Also we apply Bendixson's criterion in R n (according to [35]) to find the condition of non-existence of periodic solutions around interior equilibrium E 5 .
Finally Venturino (cf.[38][39][40]) considered recovery from the disease, a step closer to the real situation, which is not assumed herein.A future direction of this work can be well extended by keeping this factor in mind.