Direction and Stability of Bifurcating Periodic Solutions in a Delay-Induced Ecoepidemiological System

A SI-type ecoepidemiological model that incorporates reproduction delay of predator is studied. Considering delay as parameter, we investigate the effect of delay on the stability of the coexisting equilibrium. It is observed that there is stability switches, and Hopf bifurcation occurs when the delay crosses some critical value. By applying the normal form theory and the center manifold theorem, the explicit formulae which determine the stability and direction of the bifurcating periodic solutions are determined. Computer simulations have been carried out to illustrate different analytical findings. Results indicate that the Hopf bifurcation is supercritical and the bifurcating periodic solution is stable for the considered parameter values. It is also observed that the quantitative level of abundance of system populations depends crucially on the delay parameter if the reproduction period of predator exceeds the critical value.


Introduction
Ecoepidemiology is a branch in mathematical biology which considers both the ecological and epidemiological issues simultaneously. After the pioneering work of Anderson  In this model, S, I, and P represent the densities of susceptible prey, infected prey, and the predator populations, respectively. Both susceptible and infected preys contribute to the carrying capacity K , but only susceptible prey can reproduce at the intrinsic growth rate r. Disease spreads horizontally from infected to susceptible prey at a rate λ following the law of mass action. Predator preys on infected prey only and predation process follows Holling Type II 10 response function with search rate m and half-saturation constant a. Here, α is the conversion efficiency of the predator defining the increase in predator's number per unit prey consumption. μ μ 1 μ 2 represents the total death rate of infected prey where μ 1 is the natural death rate and μ 2 is the virulence of the disease. Predators consume both the susceptible and infected preys; however, the predation rate on infected prey may be very high 31 times compare to that on susceptible prey 11 . Based on the experimental observation 11 , it is assumed that predator consumes infected prey only. Predators may have to pay a cost in terms of extra mortality in the tradeoff between the easier predation and the parasitized prey acquisition, but the benefit is assumed to be greater than the cost 12, 13 . So it is assumed that consumption of infected prey contributes positive growth to the predator population. d d 1 d 2 is the total death rate of predator where d 1 is the natural death rate and d 2 is the cost due to parasitized prey acquisition. All parameters are assumed to be positive.
Reproduction of predator after consuming the prey is not instantaneous, but mediated by some time lag. Chattopadhyay and Bairagi 3 did not consider this reproduction delay, defined by the time required for the reproduction of predator after consuming the prey, in their model system. It is well recognized that introduction of reproduction delay makes the model biologically more realistic. If τ >0 is the time required for the reproduction, the model 1.1 can be written as

1.2
We study the delay-induced system 1.2 with the following initial conditions: Hopf bifurcation and its stability in a delay-induced predator-prey system have been studied by many researchers 14-19 . In this paper, we study the effect of reproduction delay on an ecoepidemiological system where predator-prey interaction follows Holling Type II response function, and find the direction and stability of the bifurcating periodic solutions, if any.
The organization of the paper is as follows. Section 2 deals with the linear stability analysis of the model system. In Section 3, direction and stability of Hopf bifurcation are presented. Numerical results to illustrate the analytical findings are presented in Section 4 and, finally, a summary is presented in Section 5.
International Journal of Differential Equations 3

Stability Analysis and Hopf Bifurcation
In epidemiology, the basic reproductive ratio R 0 , the number of new cases acquired directly from a single infected prey when introduced into a population of susceptible, plays a significant role in the spread of the disease. In particular, if R 0 < 1, the disease dies out, but if R 0 > 1, it remains endemic in the host population 20 . For the system 1.2 , the basic reproductive ratio is given by R 0 λK/μ. In ecology, on the other hand, stress is given on the stability of coexisting equilibrium point. We, therefore, concentrate on the study of the stability of the coexisting or endemic equilibrium point of the system 1.2 .
Let x t S t − S * , y t I t − I * , and z t P t − P * be the perturbed variables. Then, the system 1.2 can be expressed in the matrix form after linearization as follows:

2.2
The characteristic equation of the system 2.1 is given by International Journal of Differential Equations

2.6
Here Thus, X > 0 if m > dλK/rα. After some algebraic manipulation, Y can be written as International Journal of Differential Equations 5 So the sufficient condition for Y to be positive is Note that Z m 0 n 0 adrmS * P * /K a I * 2 is always positive. One can write, Since all the terms in the third bracket are positive, so the sufficient condition for the positivity of XY − E is dμ mα Hence, by Routh-Hurwitz criterion and using existence conditions, we state the following theorem for the stability of the interior equilibrium E * of the system 1.2 for τ 0.
where S t , I t , P t is the solution of the system 1.2 which satisfies the condition 1.3 .

6
International Journal of Differential Equations The equilibrium E * is called absolutely stable if it is asymptotically stable for all delays τ ≥ 0 and conditionally stable if it is stable for τ in some finite interval.
Note that the system 1.2 will be stable around the equilibrium E * if all the roots of the corresponding characteristic equation 2.5 have negative real parts. But 2.5 is a transcendental equation and has infinite number of roots. It is difficult to determine the sign of these infinite number of roots. Therefore, we first study the distribution of roots of the cubic exponential polynomial equation 2.5 .

2.15
This two equations give the positive values of τ and ω for which 2.5 can have purely imaginary roots. Squaring and adding, we obtain

2.16
If we assume h ω 2 , then 2. 16  Clearly, if Δ p 2 − 3q ≤ 0, then the function g h is monotonically increasing in h ∈ 0, ∞ . Thus, for s ≥ 0 and Δ ≤ 0, 2.18 has no positive roots for h ∈ 0, ∞ . On the other hand, when s ≥ 0 and Δ < 0, the equation has two real roots Obviously, g h * It follows that h * 1 and h * 2 are the local minimum and the local maximum, respectively. Hence we have the following lemma.
Proof. Noticing that s ≥ 0, h * 1 is the local minimum of g h and lim h → ∞ g h ∞, we immediately know that the sufficiency is true. So we have to prove now the necessity. In contrary, we suppose that either h * Thus, g z cannot have any positive real roots when h * 1 > 0 and g h * 1 > 0. This completes the proof.
Summarizing the above discussions, we obtain the following. ii if s ≥ 0, and Δ p 2 − 3q ≤ 0, then 2.17 has no positive root; iii if s ≥ 0, and Δ p 2 − 3q > 0, then 2.17 has positive roots if and only if h * Suppose that 2.17 has positive roots. Without loss of generality, we assume that it has three positive roots, defined by h 1 , h 2 , and h 3 , respectively. Then, 2.16 has three positive roots Thus, if we denote International Journal of Differential Equations where k 1, 2, 3; j 0, 1, 2, . . ., then ±iω k is a pair of purely imaginary roots of 2.5 . Define We reproduce the following result due to Ruan and Wei 23 to analyze 2.5 .

2.34
Thus, from Lemmas 2.7 and 2.8, we have the following theorem. i When s ≥ 0, and Δ p 2 − 3q ≤ 0, then all roots of 2.5 have negative real parts for all τ ≥ 0 and the equilibrium E * of the system 1.2 is absolutely stable for all τ ≥ 0.
ii If either s < 0 or s ≥ 0, Δ p 2 − 3q > 0, h * 1 −p √ Δ /3 > 0 and g h * 1 ≤ 0 hold, then g h has at least one positive root h k and all roots of 2.5 have negative real parts for all τ ∈ 0, τ 0 k , then the equilibrium E * of the system 1.2 is conditionally stable for iii If all the conditions as stated in (ii) and g h k / 0 hold, then the system 1.2 undergoes a Hopf bifurcation at E * when τ τ j k , j 0, 1, 2, . . . .

Direction and Stability of the Hopf Bifurcation
In the previous section, we obtained some conditions under which system 1.2 undergoes Hopf bifurcation at τ τ j j 0, 1, 2, . . . . In this section, we assume that the system 1.2 undergoes Hopf bifurcation at E * when τ τ j , that is, a family of periodic solutions bifurcate from the positive equilibrium point E * at the critical value τ τ j j 0, 1, 2, . . . . We will use the normal form theory and center manifold presented by Hassard et al. 24 to determine the direction of Hopf bifurcation, that is, to ensure whether the bifurcating branch of periodic solution exists locally for τ > τ j or τ < τ j , and determine the properties of bifurcating periodic solutions, for example, stability on the center manifold and period. Throughout this section, we always assume that system 1.2 undergoes Hopf bifurcation at the positive equilibrium E * S * , I * , P * for τ τ j and then ±iω k is corresponding purely imaginary roots of the characteristic equation.
where τ j is defined by 2.23 and ν ∈ R. Dropping the bars for simplification of notations, system 1.2 can be written as functional differential equation FDE in C C −1, 0 , R 3 aṡ

3.3
By the Riesz representation theorem, there exists a 3 × 3 matrix, η θ, ν −1 ≤ θ ≤ 0 whose elements are bounded variation functions such that International Journal of Differential Equations In fact, we can choose where δ is the Dirac delta function defined by

3.7
Then system 3.1 is equivalent toẋ where η θ η θ, 0 . Then A 0 and A * are adjoint operators. By Theorem 2.9, we know that ±iτ j ω 0 are eigenvalues of A 0 . Thus, they are also eigenvalues of A * . We first need to compute the eigenvalues of A 0 and A * corresponding to iτ j ω 0 and −iτ j ω 0 , respectively.
Suppose that q θ 1, β, γ T e iθω 0 τ j is the eigenvector of A 0 corresponding to iτ j ω 0 . Then A 0 q θ iω 0 τ j q θ . It follows from the definition of A 0 and 3.2 , 3.4 , and 3.5 that Thus, we can easily obtain

3.13
Similarly, let q * s D 1, β * , γ * T e isω 0 τ j be the eigenvector of A * corresponding to −iω 0 τ j . By the definition of A * and 3.2 , 3.3 , and 3.4 , we can compute In order to assure q * s , q θ 1, we need to determine the value of D. From 3.10 , we have

14
International Journal of Differential Equations Thus, we can choose D as

3.16
In the remainder of this section, we use the theory of Hassard et al. 24 to compute the conditions describing center manifold C 0 at ν 0. Let x t be the solution of 3.8 when ν 0. Define On the center manifold C 0 , we have z and z are local coordinates for center manifold C 0 in the direction of q * and q * . Note that W is real if x t is real. We only consider real solutions. For solution x t ∈ C 0 of 3.8 , since ν 0, we havė We rewrite this equation asż t iω 0 τ j z t g z, z , 3.21 where g z, z q * 0 f 0 z, z g 20 z 2 2 g 11 zz g 02 z 2 2 g 21 z 2 z 2 · · · . 3.22 ∈ R 3 is a constant vector. Similarly, from 3.29 and 3.32 , we obtain 3 2 In what follows, we will seek appropriate E 1 and E 2 . From the definition of A and 3.29 , we obtain where η θ η 0, θ . By 3.27 , we have

3.39
Substituting 3.34 and 3.38 into 3.36 and noticing that we obtain

3.41
International Journal of Differential Equations

19
This leads to 3.42 Solving this system for E 1 , we obtain

Numerical Simulations
In this section, we present some numerical simulations to illustrate the analytical results observed in the previous sections. We consider the following set of parameter values: For the above parameter set, the system 1.2 has a unique coexistence equilibrium point E * S * , I * , P * 20.9091, 13.6364, 22.0992 . When τ ≥ 0, the system 1.2 satisfies all conditions of the Theorem 2.9 i . Consequently, the coexistence equilibrium point E * becomes absolutely stable. Figure 1 shows the behavior of the system 1.2 when τ 0, and Figure 2 depicts the same for τ 1 and τ 25. If we change the value of m from 0.45 to 0.72 in the given parameter set, then conditions of the Theorem 2.9 ii are satisfied and the system 1.2 becomes conditionally stable around the coexistence equilibrium point E * for τ ∈ 0, τ 0 see, Figure 3 a and unstable for τ > τ 0 see, Figure 3 b .
For the given parameter set with m 0.72, one can evaluate that τ 0 2.3187 and g h k 0.3312 / 0, so the system 1.2 undergoes a Hopf bifurcation at E * when τ τ 0 following the condition iii of Theorem 2.9. We have constructed a bifurcation diagram see, Figure 4 to observe the dynamics of the system when τ varies. For this, we have run the system 1.2 for 500 time-steps and have plotted the successive maxima and minima of the prey and predator populations with τ as a variable parameter. This figure shows that the coexisting equilibrium E * is stable if τ is less than its critical value τ 0 2.3187 and unstable if τ > τ 0 and a Hopf bifurcation occurs at τ τ 0 . Using Theorem 3.1, one can determine the values of v 2 , β 2 and T 2 . For the given parameter set with m 0.72, one can evaluate that v 2 58.4107 >0 , β 2 −1.3298 <0 , and T 2 1.9364 >0 . Since v 2 > 0 and β 2 < 0, the Hopf bifurcation is supercritical and the bifurcating periodic solutions exist when τ crosses τ 0 from left to right. Also, the bifurcating periodic solution is stable as β 2 < 0 and its period increases with τ as T 2 > 0 . From the bifurcation diagram Figure 4 , it is clear that when the delay, τ, exceeds the critical value τ 0 2.3187 days approximately , the system 2.4 bifurcates from stable focus to stable limit cycle. One can also notice that the amplitude of the oscillations increases with increasing τ.

Summary
In this paper, we have studied the effects of reproduction delay on an ecoepidemiological system where predator-prey interaction follows Holling Type II response function. We have obtained sufficient conditions on the parameters for which the delay-induced system is asymptotically stable around the positive equilibrium for all values of the delay parameter and if the conditions are not satisfied, then there exists a critical value of the delay parameter below which the system is stable and above which the system is unstable. By applying the normal form theory and the center manifold theorem, the explicit formulae which determine the stability and direction of the bifurcating periodic solutions have been determined. Our analytical and simulation results show that when τ passes through the critical value τ 0 , the coexisting equilibrium E * losses its stability and a Hopf bifurcation occurs, that is, a family of periodic solutions bifurcate from E * . Also, the amplitude of oscillations increases with increasing τ. For the considered parameter values, it is observed that the Hopf bifurcation is supercritical and the bifurcating periodic solution is stable. The quantitative level of abundance of system populations depends crucially on the delay parameter if the reproduction period of predator exceeds the critical value τ 0 .