Analysis of stability and Hopf bifurcation for an eco-epidemiological model with distributed delay ∗

In this paper, the dynamical behavior of an eco-epidemiological model with distributed delay is studied. Sufficient conditions forthe asymptotical stability of all the equilibria are obtained. We prove that there exists a threshold value of the infection rate b beyond which the positive equilibrium bifurcates towards a periodic solution. We further analyze the orbital stability of the periodic orbits arising from bi- furcation by applying Poore's condition. Numerical simulation with some hypothetical sets of data has been done to support the analytical findings.


Introduction
The mathematical modelling of epidemics has become a very important subject of research after the seminal model of Kermac-McKendric on SIRS (susceptible-infectedremoved-susceptible) systems [18], in which the evolution of a disease which gets transmitted upon contact is described.Important studies in the following decades have been carried out, with the aim of controlling the effects of diseases and of developing suitable vaccination strategies [10,22,29,9,16].After the seminal models of Vito Volterra [36] and Alfred James Lotka [20] in the mid 1920s for predator-prey interactions, mutualist and competitive mechanisms have been studied extensively in the recent years by researchers [21,25,32].
In the natural world, however, species do not exist alone, it is of more biological significance to study the persistence-extinction threshold of each population in systems of two or more interacting species subjected to parasitism.Mathematical biology, namely predator-prey systems and models for transmissible diseases are major fields of study in their own right.But little attention has been paid so far to merge these two important areas of research (see [2,27,37,38,39]).Eco-epidemiology is rather a young subject of study, which tries to merge the epidemics models with some demographic issues.The first papers along these guidelines were indeed [8,24], where the dynamics of a reproducing population is studied, which is also subject to an epidemics.A model for a disease spreading among interacting populations was first described in [13].In [35], the Lotka-Volterra model is taken as the demographic basis on which to study the influence of a disease propagating in one of the two species.In an acquatic medium, a model of this kind has been first introduced in [3], while in the important domain of viral diseases affecting plankton, it is studied in [31].Also the recent ecological literature has emphasized the importance of parasites in shaping the dynamics of both plant and animal communities [11].Nowadays it has been observed that viral, bacterial, fungal and metazoan parasites can mediate host vulnerable to predation [15].
Studies on ecology and epidemiology share some common features.It is very important from both the ecological and the mathematical points of view to study ecological systems subject to epidemiological factors.A number of studies have been performed in this direction, since transmissible disease in ecological situation cannot be ignored [4,6,14].In order to study the influence of disease on an environment where two or more interacting species are present.In this paper, we shall put emphasis on such an eco-epidemiological system consisting of three populations, namely, the healthy prey (which is susceptible), the infected prey (which becomes infective by direct contact) and the predator population.
We have two populations: 1.The prey, whose total population density is denoted by N(t). 2. The predator, whose population density is denoted by P (t).
We make the following assumptions: (A 1 ) In the absence of transmissible disease the prey population grows according to logistic law with carrying capacity K (K > 0) and an intrinsic birth rate constant r (r > 0) [12,14,30], i.e., (A 2 ) In the presence of infection, the total prey population N(t) are divided into two distinct classes, namely, susceptible populations, S(t), and infected populations, [12,14,30].Therefore, at any time t, the total density of prey population is EJQTDE, 2012 No. 44, p. 2 (A 3 ) We assume the infected prey I(t) are removed by death (say its death rate is a positive constant d 1 ), or by predation before having the possibility of reproducing.(A 4 ) We assume that the disease is spread among the prey population only and the disease is not genetically inherited.The infected populations do not recover or become immune.The incidence is assumed to be the simple mass action incidence bS(t)I(t), where b > 0 is called the transmission coefficient.Hence, the SI model of the infected prey is: Numerous field studies show that infected prey are more vulnerable to predation compared with their non-infected counterpart [17,26].Lafferty and Morris [19] quantified that the predation rates on infected prey may be thirty one times higher compared to that on susceptible prey.Thus, we consider the case when the predator mainly eats the infected prey.And we are assuming in a more realistic fashion that the present level of the predator affects instantaneously the growth of the infected prey, but that the growth of the predator is influenced by the amount of infected prey in the past.More precisely, the number of predators grows depending on the weight-averaged time of the Michaelis-Menten function of I(t) over the past by means of the function Q(t) given by the integral [33] where the exponential weight function satisfies Clearly, this assumption implies that the influence of the past fades away exponentially and the number 1 δ might be interpreted as the measure of the influence of the past.So, the smaller the δ > 0, the longer the interval in the past in which the values of I are taken into account.
From the above assumptions we have the following model: where S(t), I(t) and P (t) denote the quantities of sound prey, infected prey and predator, respectively.d 2 (> 0) and h(> 0) the death rate of predator and conversion rate, respectively.EJQTDE, 2012 No. 44, p. 3 The integro-differential system (1.2) can be transformed [5,23] into the system of differential equations on the interval [0, +∞) (1.3) We understand the relationship between the two systems as follows: If (S, I, P ) : [0, +∞) → R 3 is the solution of (1.2) corresponding to continuous and bounded initial function ( S0 , Ī0 , P0 ) : [0, +∞) → R 3 , then (S, I, Q, P ) : [0, +∞) → R 4 is a solution of (1.3) with S0 = S 0 , Ī0 = I 0 , P0 = P 0 , and Conversely, if (S, I, Q, P ) is any solution of (1.3) defined on the entire real line and bounded on (−∞, 0], then Q is given by (1.1), and so (S, I, P ) satisfies (1.2).System (1.3) will be analyzed with the following initial conditions The organization of the paper is as follows: Section 2 deals with some basic results, e.g., positivity, existence of equilibria, boundedness of solutions.Section 3 is devoted to studying the dynamical behavior of the linearized system around each of the equilibria.In Section 4, we discuss the orbital stability of bifurcating limit cycle using Poore's condition.The globally asymptotical stability of predator-free equilibrium is presented in section 5. We perform a numerical analysis in Section 6 to support our analytical findings.The paper ends with discussion on the results obtained in the previous sections.

Some basic results
It is important to show positivity and boundedness for the system (1.3) as they represent populations.Positivity implies that populations survives and boundedness may be interpreted as a natural restriction to growth as a consequence of limited resources.
In this section, we present some basic results, such as the positive invariance of system (1.3), the existence of equilibria, the boundedness of solutions.

Equilibria
The model (1.3) has the following equilibria: (i) the predator-free equilibrium E 1 (S
Proof.Let us put Eq. (1.3) in a vector form by setting and where G : R 4 → R 4 and G ∈ C ∞ (R 4 ).Then Eq. (1.3) becomes where Due to the well-known theorem by Nagumo [34] and the solution of Eq. (2.3) with + for all t > 0.

Boundedness of solutions
Now, let us prove that the solutions of the system (1.3) are bounded for t ≥ 0.
As is obvious for system (1.3), we have lim sup t→+∞ S(t) ≤ K. Then there is a T > 0 such that for any sufficiently small ǫ > 0 we have S(t) ≤ K + ǫ for t > T .
. Recall that S(t) ≤ K + ǫ for t > T .Then there exists an M 1 , depending only on the parameters of system (1.3), such that V 1 < M 1 for t > T .Then I(t) has an ultimately above bound M 1 .It follows from the third equation of system (1.3) that Q(t) has an ultimately above bound M 1 m .And from the last equation of system (1.3), we can find that P (t) has an ultimately above bound hM

Local stability
In this section, we discuss the local stability of each equilibrium of the system (1.3).
Let Ē( S, Ī, Q, P ) be any arbitrary equilibrium of (1.3).The variational matrix of the system at Ē is given by .
The variational matrix at E 1 (S 1 , I 1 , 0, 0) is given by , which gives the following characteristic equation in λ as: The equation has two roots with negative real part because their two coefficients are positive.Therefore, the stability of E 1 is determined by the sign of the real part of the following equation Now this gives rise to the following cases: 2) has two eigenvalues with real part negative, consequently, the predator-free equilibrium E 1 is stable.
Case 2. When d 2 = h, there is one eigenvalue zero and the other with real part negative.System (1.3) has a saddle-node bifurcation at E 1 .
The variational matrix for the equilibrium E 2 is given by .
The above variational matrix gives the following characteristic equation in λ: where By the Routh-Hurwitz criterion, it follows that all eigenvalues of (3.3) have negative real parts if and only if p 1 > 0, p 3 > 0, p 4 > 0 and The first factor has two pure imaginary roots λ 1 , λ 2 = λ1 at b = b * , where λ 1 = iω 0 = i p 3 p 1 .Also, the following conditions hold: ) ) where The following condition has to be checked: 3) and calculating the derivative, we have Thus the conditions for Hopf bifurcation are verified.We may summarize the above discussion in the form of following proposition: Next we discuss the orbital stability of the limit cycle arising out of Hopf-bifurcation.We apply here the Poore's condition [28] for orbital stability.EJQTDE, 2012 No. 44, p. 9 We apply Poore's condition for verification of the orbital stability of the Hopf bifurcating limit cycle.Let a real, n-dimensional (n ≥ 2), first order system of autonomous differential equations be of the form where ⊤ denotes a vector of m-real parameters.Assume that there exist a combination of the parameters, say, v 0 , a critical point a 0 such that the variational matrix F(a 0 , v 0 ) has exactly two, nonzero, purely imaginary eigenvalues, say ±iω 0 with ω 0 > 0 and other n − 2 eigenvalues with nonzero real parts.To vary one, some or all the m parameters, an m-dimensional vector function f (ε) is introduced with the property that f (0) = 0, and hence we confine our analysis on the following system of ODE: It follows from the definition of F(x, ε) that F(a 0 , 0) = 0 and eigenvalues of F x (a 0 , 0) are the same as those of F(a 0 , v 0 ).It is assumed that where k ≥ 3, D is a domain in R n containing a 0 and ε 0 > 0. As there is no nonzero eigenvalue of the variational matrix, F x (a 0 , 0) = 0 ∈ R n , and hence the implicit function theorem guarantees the existence of a critical point a ε which is k-times continuously differentiable in ε and satisfies F(a ε , ε) = 0 for ε in a small neighborhood of ε = 0. Using this definition of a ε , we introduce a change of variables which is similar to that used by K.O.Friedrichs [7]: This reduces the differential equation (4.2) in the following form where Thus the problem of periodic solutions of (4.2) is reduced to a perturbation problem in the small parameter µ ′ .Now, the stability information of the bifurcating periodic orbits is contained in the following theorem of Poore, coupled with an algebraic expression, which completely reduces the determination of stability to an algebraic problem.By the assumptions of Theorem 2.3 in [28], the differential equation in (4.4) is continuously differentiable in µ ′ , when δ = δ(µ ′ ) and η = η(µ ′ ) and in the function y in a neighborhood of the periodic orbit.Thus the existing periodic orbit will be asymptotically orbitally stable with asymptotic phase if n − 1 of the characteristic multipliers of the variational equation EJQTDE, 2012 No. 44, p. 10 have moduli less than one.The information developed by Poore about the modulus of each of the characteristic multipliers is given by the following lemma: Lemma 4.1.Assume that there exists a nonconstant bifurcated periodic solution, i.e., the hypothesis in Theorem 2.3 in [28] is satisfied.Then the n-characteristic multipliers, denoted by ̺ l (µ ′ ) for 1 ≤ l ≤ n, are continuous functions of µ ′ at µ ′ = 0 and satisfy the relations , where β(µ ′ ) is continuous at µ ′ = 0 and β(0) = 0. Here, T is the period of oscillation of the periodic solution y(s, µ ′ ) of Theorem 2.3 in [28], and the λ 0 l 's denote the n − 2 eigenvalues (continuing multiplicities) of the variational matrix F x (a 0 , 0) = Fx (a 0 , 0) with nonzero real parts.
The following lemma reduces the above condition of stability to an algebraic expression as follows: Lemma 4.2.[28] Let F(x, ε) satisfy the hypotheses in Theorem 2.3 in [28] and let u and v denote the left and right eigenvector, respectively, for the eigenvalues +iω 0 (ω 0 > 0) of the variational matrix A 0 .If u and v are normalized in the sense that uv = 1 (u is row and v is a column vector), then where the complex conjugate of v, and the definition of Φ(s) can be found in (2.6) of [28].
Written out in component form, the above expression reduces to We rewrite our system of ODE (1.3) in the following form: Now all the second-and third-order derivatives of F i (i = 1, 2, 3, 4) are as follows: , (mP +I) 4 .And q 11 q 12 q 13 q 14 q 21 q 22 q 23 q 24 q 31 q 32 q 33 q 34 q 41 q 42 q 43 q 44 where EJQTDE, 2012 No. 44, p. 12 ] and ∆ = −( rS 2 K + 2iω 0 )q 11 − bI 2 q 12 .Putting ω 0 = 0 in the above expressions, we get the components of Q = (J E 2 ) −1 .Now we are to calculate the left and right eigenvector u and v respectively, of the variational matrix J E 2 , i.e., we have to find out a row vector u = (u 1 , u 2 , u 3 , u 4 ), such that uJ E 2 = iω 0 u, and a column vector Proceeding in the above way and solving the set of equations, we find the left eigenvector u = (u 1 , u 2 , u 3 , u 4 ), where where ω 0 is given by Eq. (3.10).Now, writing the expression of (4.5) in detail, we have the following: First term: Second term: In the above terms, Putting the values of F xxx , F xx , u, v and components of the matrix Q in terms of the parameters of the model, the sign of the real part of the resulting expression (4.5) can be deduced.This in turn indicates the orbital stability of the limit cycle arising out of Hopf bifurcation.
This completes the proof.

Numerical results
In the previous sections, we introduced the analytical tools proposed and used them for a qualitative analysis of the system obtaining some results about the dynamics of EJQTDE, 2012 No. 44, p. 15 the system.In this section, we perform a numerical analysis of the model based on the previous results.
In order to numerically experiment the effects of parameter changes in the proposed model, we have run simulations using the standard MATLAB differential equations integrator for the Runge-Kutta method, i.e., routine ODE45 or ODE15S if needed.
From the analysis of section 4, we can see that the infection rate plays a vital role in describing the behavior of the system.Keeping this mind, we have varied the infection rate, b, to observe the dynamics of the system (1.3).We select r = 1, K = 10, d 1 = 0.6, a = 0.2, m = 1, δ = 1.8, d 2 = 0.5, h = 1.5.Using these sets of parametric values and initial values S(0) = 1, I(0) = 1, Q(0) = 0.1, P (0) = 0.1, the following interesting dynamic behavior of the system was observed.Our numerical results show that the infective population will be asymptotically stable for b = 1.2 (Fig. 1).Increasing the infection rate to b = 2.2, we observe that the three species of system (1.3) enter into an oscillatory steady state from a stable situation (Fig. 2).In fact, for b = 1.In most studies of eco-epidemiological models, the predation rate is also an important parameter that plays a vital role in the dynamics of the system.Therefore, we then vary the parameter h keeping the other parameter values.We consider a hypothetical set of parametric values, where r = 1, K = 10, d 1 = 0.6, a = 0.2, m = 1, δ = 1.8, d 2 = 0.5, b = 2.2 and different values of h.We can observe that when h = 1.5 the three species periodically oscillate, i.e., Hopf bifurcation occurs around E 2 (Fig. 2); when h = 4 the three species coexist, i.e, E 2 is asymptotically stable (Fig. 3).
Select r = 1, K = 10, d 1 = 0.6, a = 0.2, m = 1, δ = 1.8, d 2 = 0.5, b = 2.2 and h = 0.4.We can see that the predator-free equilibrium E 1 (0.2727272727, 0.4229249012, 0, 0) is asymptotically stable (Fig. 4).In this paper we proposed and analyzed, both analytically and numerically, an ecoepidemiological model with distributed delay.Firstly, we obtain the positivity and boundedness of solutions.Secondly, by using of Routh-Hurwitz criterion, the asymptotical stability of the non-negative equilibria are obtained.Thirdly, we have obtained conditions for small amplitude periodic solutions bifurcating from the boundary equilibrium E 2 and the positive interior equilibrium E 2 by applying both mathematical and numerical techniques.It is observed that when the infection rate b crosses a critical value, say b * , the system (1.3) enters into Hopf bifurcation that induces periodic solutions near the boundary equilibrium E 2 .And the stability as well as the direction of bifurcation near the positive interior equilibrium E 2 is obtained by applying the algorithm due to Poore.From numerical simulations, we find that the predation rate is also an important parameter that plays a vital role in the dynamics of the system.
In conclusion, this paper emphasizes that the value of the per capita disease contact rate b, the predation rate h and the existence of the predator population have important implications for predator species persistence and the the existence of Hopf bifurcation in a simple eco-epidemiological model where the preys are infected by transmissible disease.

Theorem 3 . 3 .Remark 3 . 3 .Remark 3 . 4 .
If the positive equilibrium E 2 of the system (1.3) exists, then the system around E 2 has a Hopf bifurcation when b = b * .Remark 3.2.Theorem 3.3 shows the importance of b, the infection rate in controlling the system dynamics.The biological implication of Hopf bifurcation for E 2 is that the predator coexists with the susceptible prey and the infected prey, exhibiting oscillatory balance behavior.The period τ 1 of the bifurcating periodic orbits for close to b = b * is given by

(4. 26 )
EJQTDE, 2012 No. 44, p. 11 where the repeated indices within each term imply a sum from 1 to n and all the derivatives of F are evaluated at the equilibrium a 0 .The sign of the real and imaginary parts of the right hand side of expression in Lemma 4.2 are independent of the choice of b(µ ′ ) in Uy(s, µ ′ ) = Φ(s)b(µ ′ ) and the eigenvectors u and v so long as b(0) • b(0) > 0 and uv = 1.So positivity of real part of the above expression in parenthesis really indicates the orbital stability of the periodic solution arising out of Hopf bifurcation.
For positivity of E 1 , we assume bK > d 1 .Positivity of I 2 implies that bmKh + d 2 a > d 1 mh + ah.If h > d 2 , we can get the positivity of Q 2 and P 2 .In order to insure the positivity of S 2 we must have d 1 mh + ah > d 2 a.Hence, for positivity of E 2 , we assume bmKh + d 2 a > d 1 mh + ah and h > d 2 .