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

Authors: N. Bairagi
Journal: International Journal of Differential Equations
 (2011)

---

## Abstract

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.


---

## Body

## 1. Introduction

Ecoepidemiology is a branch in mathematical biology which considers both the ecological and epidemiological issues simultaneously. After the pioneering work of Anderson and May [1], literature in the field of ecoepidemiology has grown enormously [2–9]. Chattopadhyay and Bairagi [3] studied the following ecoepidemiological model with mα=θ: (1.1)dSdt=rS(1-S+IK)-λIS,dIdt=λIS-mIPa+I-μI,dPdt=mαIPa+I-dP.
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(=d1+d2) is the total death rate of predator where d1 is the natural death rate and d2 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)dSdt=rS(1-S+IK)-λIS,dIdt=λIS-mIPa+I-μI,dPdt=mαI(t-τ)P(t-τ)a+I(t-τ)-dP.
We study the delay-induced system (1.2) with the following initial conditions: (1.3)S(θ)=ψ1(θ)≥0,I(θ)=ψ2(θ)≥0,P(θ)=ψ3(θ)≥0,θ∈(-τ,0].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. Section2 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.


## 2. Stability Analysis and Hopf Bifurcation

In epidemiology, the basic reproductive ratioR0, 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 R0<1, the disease dies out, but if R0>1, it remains endemic in the host population [20]. For the system (1.2), the basic reproductive ratio is given by R0=λ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). The ecoepidemiological system (1.2) has a unique interior equilibrium point E*(S*,I*,P*), where S*=K-(ad(r+λK)/r(mα-d)),I*=ad/(mα-d), and P*=(1/m)(a+I*)(λS*-μ). Note that I* exists if m>d/α, S* exists if m>d/α+(ad(r+λK)/rKα) and P* exists if m>(d/α)+(adλ(r+λK)/rα(λK-μ))=(d/α)+(ad(r+λK)/rKα-(rαμ/λ)) with λ>μ/K. Thus, the conditions for coexisting equilibrium point E* are(i)
λ>μ/K, that is, R0>1,(ii)
m>(d/α)+(adλ(r+λK)/rα(λK-μ)).Letx(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.1)ddt(x(t)y(t)z(t))=A1′(x(t)y(t)z(t))+A2′(x(t-τ)y(t-τ)z(t-τ)),
where   (2.2)A1′=(-rS*K-(λ+rK)S*0λI*mI*P*(a+I*)2-dα00-d),A2′=(0000000αamP*(a+I*)2d).
The characteristic equation of the system (2.1) is given by (2.3)|A1′+A2′e-ξτ-ξI|=0,
that is, (2.4)Φ(ξ,τ)=ξ3+[A+Be-ξτ]ξ2+[C+De-ξτ]ξ+E+Fe-ξτ=0,
where (2.4)′A=rKS*+d-mI*P*(a+I*)2,B=-d,C=-mdI*P*(a+I*)2+drS*K-rmS*I*P*K(a+I*)2+(rK+λ)λS*I*,D=mdI*P*(a+I*)2-drS*K+admP*(a+I*)2,E=-drmS*I*P*K(a+I*)2+(rK+λ)λdS*I*,F=drmS*I*P*K(a+I*)2+adrmS*P*K(a+I*)2-(rK+λ)λdS*I*.
Equation (2.4) can be written as (2.5)Φ(ξ,τ)=ξ3+m2ξ2+m1ξ+m0+(n2ξ2+n1ξ+n0)e-ξτ=0,
where (2.5)′m2=A,n2=B,m1=C,n1=D,m0=E,n0=F,Σni2≠0,(i=0,1,2).
For τ=0, (2.5) becomes (2.6)Φ(ξ,0)=ξ3+(m2+n2)ξ2+(m1+n1)ξ+(m0+n0)=ξ3+Xξ2+Yξ+Z=0.
Here (2.7)X=m2+n2=rS*K-mP*I*(a+I*)2=(rK-dλmα)S*+dμmα.
Thus, X>0 if m>dλK/rα. After some algebraic manipulation, Y can be written as (2.8)Y=m1+n1=rmaS*P*K(a+I*)2+madP*(a+I*)2+S*{rμK+λI*(λ+rK)-rλS*K}.
So the sufficient condition for Y to be positive is (2.9)rμK+λI*(λ+rK)>rλS*Korm<dα+2λad(r+λK)rα(λK-μ).
Note that Z=m0+n0=adrmS*P*/K(a+I*)2 is always positive. One can write, (2.10)XY-Z=[(rK-dλmα)S*][{rmaS*P*K(a+I*)2+madP*(a+I*)2}+{rμS*K+λS*I*(λ+rK)-rλS*2K}]+dμmα[madP*(a+I*)2]+dμmα[rμS*K+λS*I*(λ+rK)-λrS*2K]+dμmαrmaS*P*K(a+I*)2-rmadS*P*K(a+I*)2.
Since all the terms in the third bracket are positive, so the sufficient condition for the positivity of (XY-E) is (2.11)dμmαrmaS*P*K(a+I*)2>rmadS*P*K(a+I*)2orm<μα.
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.Theorem 2.1.
If(i)
R0>1orλ>μ/K,(ii)
m̲<m<m̅,
wherem̲=max[dλK/rα,(d/α)+(adλ(r+λK)/rα(λK-μ))] and m̅=min[μ/α,(d/α)+(2λad(r+λK)/rα(λK-μ))], then the system (1.2) is locally asymptotically stable without delay around the positive interior equilibrium E*.We now reproduce some definitions given by [21, 22].Definition 2.2.
The equilibriumE* is called asymptotically stable if there exists a δ>0 such that
(2.12)sup⁡-τ≤θ≤0⁡[|ψ1(θ)-S*|,|ψ2(θ)-I*|,|ψ3(θ)-P*|]<δ
implies that
(2.13)lim⁡t→∞⁡(S(t),I(t),P(t))=(S*,I*,P*),
where (S(t),I(t),P(t)) is the solution of the system (1.2) which satisfies the condition (1.3).Definition 2.3.
The equilibriumE* 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).We know thatiω(ω>0) is a root of (2.5) if and only if ω satisfies (2.14)-iω3-ω2m2+m1iω+m0+(-n2ω2+n1iω+n0)(cos⁡⁡ωτ-isin⁡⁡ωτ)=0.
Separating real and imaginary parts, we get(2.15)m2ω2-m0=-n2ω2cos⁡⁡ωτ+n1ωsin⁡⁡ωτ+n0cos⁡⁡ωτ,ω3-m1ω=n2ω2sin⁡⁡ωτ+n1ωcos⁡⁡ωτ-n0sin⁡⁡ωτ.
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)ω6+pω4+qω2+s=0,
where (2.16)′p=m22-2m1-n22,q=m12-2m0m2+2n0n2-n12,s=(m02-n02).
If we assume h=ω2, then (2.16) reduces to (2.17)h3+ph2+qh+s=0.
Denote (2.18)g(h)=h3+ph2+qh+s.
Note that g(0)=s and limh→+∞g(h)=+∞. Thus, if s<0, then (2.18) has at least one positive root.From (2.18), we have (2.19)g′(h)=3h2+2ph+q.
Clearly, if Δ=p2-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 (2.20)3h2+2ph+q=0
has two real roots (2.21)h1*=-p+Δ3,h2*=-p-Δ3.
Obviously, g″(h1*)=2Δ>0 and g″(h2*)=-2Δ<0. It follows that h1* and h2* are the local minimum and the local maximum, respectively. Hence we have the following lemma.Lemma 2.4.
Suppose thats≥0 and Δ>0.  Then (2.17) has positive roots if and only if h1*>0,g(h1*)≤0.Proof.
Noticing thats≥0,  h1* is the local minimum of g(h) and limh→+∞g(h)=+∞, we immediately know that the sufficiency is true. So we have to prove now the necessity. In contrary, we suppose that either h1*≤0 or h1*>0 and g(h1*)>0. Since g(h) is increasing for h≥h1* and g(0)=s≥0, it follows that g(h) has no positive real roots for h1*≤0 and g(h1*)>0. If h1*>0 and g(h1*)>0, since h2* is the local maximum value, it follows that g(h1*)<g(h2*). Thus, g(z) cannot have any positive real roots when h1*>0 and g(h1*)>0. This completes the proof.Summarizing the above discussions, we obtain the following.Lemma 2.5.
One has the following results on the distribution of roots of (2.17). (i)
Ifs<0, then (2.17) has at least one positive root;(ii)
ifs≥0, and Δ=p2-3q≤0, then (2.17) has no positive root;(iii)
ifs≥0, and Δ=p2-3q>0, then (2.17) has positive roots if and only if h1*=(-p+Δ)/3>0 and g(h1*)≤0, where g(z)=h3+ph2+qh+s.Suppose that (2.17) has positive roots. Without loss of generality, we assume that it has three positive roots, defined by h1, h2, and h3, respectively. Then, (2.16) has three positive roots ω1=h1,ω2=h2, and ω3=h3.From (2.15), we have (2.22)cos⁡⁡ωkτ=(n1-n2m2)ω4+(n2m0+n0m2-n1m1)ω2-n0m0(n2ω2-n0)2+n12w2,k=1,2,3.
Thus, if we denote (2.23)τkj=1ωk{arccos⁡⁡[(n1-n2m2)ωk4+(n2m0+n0m2-n1m1)ωk2-n0m0(n2ωk2-n0)2+n12wk2]+2jπ},
where k=1,2,3;  j=0,1,2,…, then ±iωk is a pair of purely imaginary roots of (2.5). Define (2.24)τ0=τk0(0)=min⁡k∈{1,2,3}⁡{τk(0)},ω0=ωk0.
We reproduce the following result due to Ruan and Wei [23] to analyze (2.5).Lemma 2.6.
Consider the exponential polynomial(2.25)P(ξ,e-ξτ1,…,e-ξτm)=ξn+p1(0)ξn-1+⋯+pn-1(0)ξ+pn(0)+[p1(1)ξn-1+⋯+pn-1(1)ξ+pn(1)]e-ξτ1+⋯+[p1(m)ξn-1+⋯+pn-1(m)ξ+pn(m)]e-ξτm,
where τi≥0(i=1,2,…,m) and pj(i), (i=0,1,2,…,m;j=1,2,…,n) are constants. As (τ1,τ2,…,τm) vary, the sum of the order of zeros of P(ξ,e-ξτ1,…,e-ξτm) on the open right half hand can change only if a zero appears on or crosses the imaginary axis.Using Lemmas2.5 and 2.6, we can easily obtain the following results on the distribution of roots of the transcendental (2.5).Lemma 2.7.
For the third degree exponential polynomial equation (2.5), one has(i)
ifs≥0, and Δ=p2-3q≤0, then all roots with positive real parts of (2.5) have the same sum as those of the polynomial equation (2.6) for all τ≥0,(ii)
if eithers<0 or s≥0, Δ=p2-3q>0,h1*=(-p+Δ)/3>0 and g(h1*)≤0, then all roots with positive real parts of (2.5) have the same sum as those of the polynomial equation (2.6) for all τ∈[0,τ0).Let(2.26)ξ(τ)=η(τ)+iω(τ),
where η and ω are real, be the roots of (2.5) near τ=τk(j) satisfying (2.27)η(τk(j))=0,ω(τk(j))=ωk.
Then the following transversality condition holds.Lemma 2.8.
Suppose thathk=ωk2 and g′(hk)≠0, where g(h) is defined by (2.18). Then,
(2.28)ddτ[Re{ξ(τk(j))}]≠0,
and the sign of (d/dτ)[Re{ξ(τk(j))}] is consistent with that of g′(hk).Proof.
Differentiating (2.5) with respect to τ, we obtain
(2.29)[(3ξ2+2m2ξ+m1)+e-ξτ{2n2ξ+n1-τ(n2ξ2+n1ξ+n0)}]dξdτ=ξ(n2ξ2+n1ξ+n0)e-ξτ.
This gives
(2.30)(dξdτ)-1=(3ξ2+2m2ξ+m1)eξτξ(n2ξ2+n1ξ+n0)+2n2ξ+n1ξ(n2ξ2+n1ξ+n0)-τξ.
It follows from (2.15) that
(2.31)[ξ(n2ξ2+n1ξ+n0)]τ=τk(j)=-n1ωk2+i(n0ωk-n2ωk3),[(3ξ2+2m2ξ+m1)eξτ]τ=τk(j)=[(m1-3ωk2)cos⁡⁡ωkτk(j)-2m2ωksin⁡⁡ωkτk(j)]+i[(m1-3ωk2)sin⁡⁡ωkτk(j)+2m2ωkcos⁡⁡ωkτk(j)],[2n2ξ+n1]τ=τk(j)=n1+i2n2ωk.
Using (2.31) in (2.30), we get
(2.32)ddτ[Re{ξ(τ)}]τ=τk(j)-1=Re{(3ξ2+2m2ξ+m1)eξτξ(n2ξ2+n1ξ+n0)}τ=τk(j)+Re{2n2ξ+n1ξ(n2ξ2+n1ξ+n0)}τ=τk(j)-Re{τξ}τ=τk(j)=1Λ[ωk(m1-3ωk2){-n1ωkcos⁡⁡ωkτk(j)+(n0-n2ωk2)sin⁡⁡ωkτk(j)}+2m2ωk2{n1ωksin⁡⁡ωkτk(j)+(n0-n2ωk2)cos⁡⁡ωkτk(j)}-n12ωk2+2n2ωk2(n0-n2ωk2){-n1ωkcos⁡⁡ωkτk(j)+(n0-n2ωk2)sin⁡⁡ωkτk(j)}]=1Λ[ωk(m1-3ωk2)(m1ωk-ωk3)+2m2ωk2(m2ωk2-m0)-n12ωk2+2n2ωk2(n0-n2ωk2)]=1Λ[3ωk6+2(m22-2m1-n22)ωk4+(m12-2m0m2-n12+2n0n2)ωk2]=1Λ[3ωk6+2pωk4+qωk2]=1Λhkg′(hk),
where Λ=n12ωk4+(n0-n2ωk2)2(>0). Thus, we have
(2.33)sign⁡{ddτReξ(τ)}τ=τk(j)=sign⁡{ddτReξ(τ)}τ=τk(j)-1=sign⁡{hkΛg′(hk)}≠0.
Since Λ,hk are positive, we conclude that the sign of {(d/dτ)Reξ(τ)}τ=τk(j) is determined by that of g′(hk). This proves the lemma.From(2.5)′ and (2.16)′, we have (2.34)p=m22-2m1-n22=A2-2C2-B2,q=m12-2m0m2+2n0n2-n12=C2-2AE+2BF-D2,s=m02-n02=E2-F2.
Thus, from Lemmas 2.7 and 2.8, we have the following theorem.Theorem 2.9.
Letmi,ni(i=0,1,2);p,q,s and τj are defined by (2.5)′, (2.34), and (2.23), respectively. Suppose that conditions of Theorem (2.1) hold. Then the following results hold.(i)
Whens≥0, and Δ=p2-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 eithers<0 or s≥0,Δ=p2-3q>0,h1*=(-p+Δ)/3>0 and g(h1*)≤0 hold, then g(h) has at least one positive root hk and all roots of (2.5) have negative real parts for all τ∈[0,τk(0)), then the equilibrium E* of the system (1.2) is conditionally stable for τ∈[0,τk(0)).(iii)
If all the conditions as stated in (ii) andg′(hk)≠0 hold, then the system (1.2) undergoes a Hopf bifurcation at E* when τ=τk(j), (j=0,1,2,…).


## 3. 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.Letx1=S-S*, x2=I-I*,  x3=P-P*, x̅i(t)=xi(τt), τ=τ(j)+ν, 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],R3) as (3.1)ẋ(t)=Lν(xt)+f(ν,xt),
where x(t)=(x1,x2,x3)T∈R3, and Lν:C→R,  f:R×C→R are given, respectively, by (3.2)Lν(ϕ)=(τ(j)+ν)(-rS*K-(rK+λ)S*0λI*mI*P*(a+I*)2-dα00-d)(ϕ1(0)ϕ2(0)ϕ3(0))+(τ(j)+ν)(0000000αamP*(a+I*)2d)(ϕ1(-1)ϕ2(-1)ϕ3(-1)),(3.3)f(ν,ϕ)=(τ(j)+ν)(-rK(ϕ12(0)+ϕ1(0)ϕ2(0))-λϕ1(0)ϕ2(0)λϕ1(0)ϕ2(0)-mϕ2(0)ϕ3(0)a+ϕ2(0)αmϕ2(-1)ϕ3(-1)a+ϕ2(-1)).
By the Riesz representation theorem, there exists a (3×3) matrix, η(θ,ν)(-1≤θ≤0) whose elements are bounded variation functions such that(3.4)Lνϕ=∫-10dη(θ,ν)ϕ(θ),forϕ∈C.
In fact, we can choose(3.5)η(θ,ν)=(τ(j)+ν)(-rS*K-(rK+λ)S*0λI*mI*P*(a+I*)2-dα00-d)δ(θ)-(τ(j)+ν)(0000000αamP*(a+I*)2d)δ(θ+1),
where δ is the Dirac delta function defined by (3.6)δ(θ)={0,θ≠0,1,θ=0.Forϕ∈C1([-1,0],R3), define the operator A(ν) as(3.7)A(ν)ϕ(θ)={dϕ(θ)dθ,θ∈[-1,0),∫-10dη(ν,s)ϕ(s),θ=0,R(ν)ϕ(θ)={0,θ∈[-1,0),f(ν,ϕ),θ=0.
Then system (3.1) is equivalent to (3.8)ẋ(t)=A(ν)xt+R(ν)xt,
where xt(θ)=x(t+θ)forθ∈[-1,0].Forψ∈C1([0,1],(R3)*), define(3.9)A*ψ(s)={-dψ(s)ds,s∈(0,1]∫-10dηT(t,0)ψ(-t),s=0
and a bilinear inner product(3.10)〈ψ(s),ϕ(θ)〉=ψ̅(0)ϕ(0)-∫-10∫ξ=0θψ̅(ξ-θ)dη(θ)ϕ(ξ)dξ,
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 thatq(θ)=(1,β,γ)Teiθω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(3.11)τ(j)(iω0+rS*K(rK+λ)S*0-λI*iω0-mI*P*(a+I*)2dα0-αamP*(a+I*)2e-iω0τ(j)iω0+d-de-iω0τ(j))q(0)=(000).
Thus, we can easily obtain(3.12)q(0)=(1,β,γ)T,
where(3.13)β=-iω0K+rS*(r+λK)S*,γ=-αamP*(iω0K+rS*)e-iω0τ(j)(r+λK)S*(iω0+d-de-iω0τ(j)).
Similarly, let q*(s)=D(1,β*,γ*)Teisω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(3.14)q*(s)=D(1,β*,γ*)eisω0τ(j)=D(1,-iω0K+rS*λKI*,d(iω0K-rS*)α(-iω0+d-de-iω0τ(j)))eisω0τ(j).
In order to assure 〈q*(s),q(θ)〉=1, we need to determine the value of D. From (3.10), we have(3.15)〈q*(s),q(θ)〉=D̅(1,β̅*,γ̅*)(1,β,γ)T-∫-10∫ξ=0θD̅(1,β̅*,γ̅*)e-iω0τ(j)(ξ-θ)dη(θ)(1,β,γ)Teiω0ξτ(j)dξ=D̅{1+ββ̅*+γγ̅*-∫-10(1,β̅*,γ̅*)θeiω0θτ(j)dη(θ)(1,β,γ)T}=D̅{1+ββ̅*+γγ̅*+τ(j)γ̅*αβamP*(a+I*)2e-iω0τ(j)}.
Thus, we can choose D̅ as (3.16)D̅=11+ββ̅*+γγ̅*+τ(j)γ̅*(αβamP*/(a+I*)2)e-iω0τ(j),∴D=11+β̅β*+γ̅γ*+τ(j)γ*(αβamP*/(a+I*)2)eiω0τ(j).
In the remainder of this section, we use the theory of Hassard et al. [24] to compute the conditions describing center manifold C0 at ν=0. Let xt be the solution of (3.8) when ν=0. Define (3.17)z(t)=〈q*,xt〉,W(t,θ)=xt(θ)-2Re{z(t)q(θ)}.
On the center manifold C0, we have (3.18)W(t,θ)=W(z(t),z̅(t),θ),
where (3.19)W(z,z̅,θ)=W20(θ)z22+W11(θ)zz̅+W02(θ)z̅22+W30(θ)z36+⋯,z and z̅ are local coordinates for center manifold C0 in the direction of q* and q̅*. Note that W is real if xt is real. We only consider real solutions. For solution xt∈C0 of (3.8), since ν=0, we have (3.20)ż(t)=iω0τ(j)z+q̅*(0)f(0,W(z,z̅,0)+2Re{zq(θ)})≝iω0τ(j)z+q̅*(0)f0(z,z̅).
We rewrite this equation as (3.21)ż(t)=iω0τ(j)z(t)+g(z,z̅),
where (3.22)g(z,z̅)=q̅*(0)f0(z,z̅)=g20z22+g11zz̅+g02z̅22+g21z2z̅2+⋯.
We have xt(θ)=(x1t(θ),x2t(θ),x3t(θ)) and q(θ)=(1,β,γ)Teiθω0τ(j), so from (3.17) and (3.19) it follows that (3.23)xt(θ)=W(t,θ)+2Re{z(t)q(t)}=W20(θ)z22+W11(θ)zz̅+W02(θ)z̅22+(1,β,γ)Teiω0τ(j)θz+(1,β̅,γ̅)Te-iω0τ(j)θz̅+⋯
and then we have(3.24)x1t(0)=z+z̅+W201(0)z22+W11(1)(0)zz̅+W02(1)(0)z̅22+⋯,x2t(0)=βz+β̅z̅+W20(2)(0)z22+W11(2)(0)zz̅+W02(2)(0)z̅22+⋯,x3t(0)=γz+γ̅z̅+W20(3)(0)z22+W11(3)(0)zz̅+W02(3)(0)z̅22+⋯,x1t(-1)=ze-iω0τ(j)+z̅eiω0τ(j)+W20(1)(-1)z22+W11(1)(-1)zz̅+W02(1)(-1)z̅22+⋯,x2t(-1)=βze-iω0τ(j)+β̅z̅eiω0τ(j)+W20(2)(-1)z22+W11(2)(-1)zz̅+W02(2)(-1)z̅22+⋯,x3t(-1)=γze-iω0τ(j)+γ̅z̅eiω0τ(j)+W20(3)(-1)z22+W11(3)(-1)zz̅+W02(3)(-1)z̅22+⋯.It follows together with (3.3) that (3.25)g(z,z̅)=q̅*(0)f0(z,z̅)=q̅*(0)f(0,xt)=τ(j)D̅(1,β̅*,γ̅*)(-rK(x1t2(0)+x1tx2t(0))-λx1t(0)x2t(0)λx1t(0)x2t(0)-mx2t(0)x3t(0)a+x2t(0)αmx2t(-1)x3t(-1)a+x2t(-1))=z22{2τ(j)D[-(rK+β(r+λK)K)+β̅*(βλ-mβγa)+γ̅*mαβγe-2iω0τ(j)a]}+zz̅{2τ(j)[D-(rK+(r+λK)KRe{β})+β̅*(λRe{β}-maRe{βγ̅})+γ̅*mαaRe{βγ̅}(rK+(r+λK)KRe{β})](rK+(r+λK)KRe{β})}+z̅22{2τ(j)D[-(rK+β̅(r+λK)K)+β̅*(β̅λ-mβ̅γ̅a)+γ̅*mαβ̅γ̅e2iω0τ(j)a(r+λK)K]}+z2z̅2{-τ(j)DrK{4W11(1)(0)+2W20(1)(0)}-τ(j)D(r+λK-β̅*λK)K{2βW11(1)(0)+2W11(2)(0)+β̅W20(1)(0)+W20(2)(0)}-τ(j)Dβ̅*ma{2a2βW11(3)(0)+2γW11(2)(0)+β̅W20(3)(0)+γ̅W20(2)(0)-2a(β2γ̅+2ββ̅γ)}+τ(j)Dγ̅*mαa[2e-iω0τ(j){βW11(3)(-1)+γW11(2)(-1)}+eiω0τ(j){β̅W20(3)(-1)+γ̅W20(2)(-1)}-2β(βγ̅+2β̅γ)e-iω0τ(j)]τ(j)DrK}⋯.Comparing the coefficients with (3.22), we have (3.26)g20=2τ(j)D[-(rK+β(r+λK)K)+β̅*(βλ-mβγa)+γ̅*mαβγe-2iω0τ(j)a],g11=2τ(j)D[-(rK+(r+λK)KRe{β})+β̅*(λRe{β}-maRe{βγ̅})+γ̅*mαaRe{βγ̅}],g02=2τ(j)D[-(rK+β̅(r+λK)K)+β̅*(β̅λ-mβ̅γ̅a)+γ̅*mαβ̅γ̅e2iω0τ(j)a],g21=-τ(j)DrK{4W11(1)(0)+2W20(1)(0)}-τ(j)D(r+λK-β̅*λK)K{2βW11(1)(0)+2W11(2)(0)+β̅W20(1)(0)+W20(2)(0)}-τ(j)Dβ̅*ma{2βW11(3)(0)+2γW11(2)(0)+β̅W20(3)(0)+γ̅W20(2)(0)-2a(β2γ̅+2ββ̅γ)}+τ(j)Dγ̅*mαa[2e-iω0τ(j){βW11(3)(-1)+γW11(2)(-1)}+eiω0τ(j){β̅W20(3)(-1)+γ̅W20(2)(-1)}-2β(βγ̅+2β̅γ)e-iω0τ(j)].
Since there are W20(θ) and W11(θ) in g21, we still need to compute them.From (3.8) and (3.17), we have(3.27)Ẇ=ẋt-żq-z̅̇q̅={AW-2Re{q̅*(0)f0q(θ)},θ∈[-1,0),AW-2Re{q̅*(0)f0q(θ)}+f0,θ=0,≝AW+H(z,z̅,θ),
where (3.28)H(z,z̅,θ)=H20(θ)z22+H11(θ)zz̅+H02(θ)z̅22+⋯.
Substituting the corresponding series into (3.27) and comparing the coefficients, we obtain (3.29)(A-2iω0τj)W20(θ)=-H20(θ),AW11(θ)=-H11(θ).
From (3.27), we know that for θ∈[-1,0),(3.30)H(z,z̅,θ)=-q̅*(0)f0q(θ)-q*(0)f̅0q̅(θ)=-g(z,z̅)q(θ)-g̅(z,z̅)q̅(θ).
Comparing the coefficients with (3.28), we get (3.31)H20(θ)=-g20q(θ)-g̅02q̅(θ),(3.32)H11(θ)=-g11q(θ)-g̅11q̅(θ).
From (3.29) and (3.31) and the definition of A, it follows that (3.33)Ẇ20(θ)=2iω0τ(j)W20(θ)+g20q(θ)+g̅02q̅(θ).
Notice that q(θ)=(1,β,γ)Teiω0τ(j)θ, hence (3.34)W20(θ)=g20ω0τ(j)q(0)eiω0τ(j)θ+ig̅023ω0τ(j)q̅(0)e-iω0τ(j)θ+E1e2iω0τ(j)θ,
where E1=(E1(1),E1(2),E1(3))∈R3 is a constant vector. Similarly, from (3.29) and (3.32), we obtain (3.35)W11(θ)=-ig11ω0τ(j)q(0)eiω0τ(j)θ+ig̅11ω0τ(j)q̅(0)e-iω0τ(j)θ+E2,
where E2=(E2(1),E2(2),E2(3))∈R3 is also a constant vector.In what follows, we will seek appropriateE1 and E2. From the definition of A and (3.29), we obtain (3.36)∫-10dη(θ)W20(θ)=2iω0τ(j)W20(0)-H20(0),(3.37)∫-10dη(θ)W11(θ)=-H11(0),
where η(θ)=η(0,θ).By (3.27), we have (3.38)H20(0)=-g20q(0)-g̅02q̅(0)+2τ(j)(-rK-β(rK+λ)βλ-mβγamαβγae-2iω0τ(j)),(3.39)H11(0)=-g11q(0)-g̅11q̅(0)+2τ(j)(-rK-(rK+λ)Re{β}λRe{β}-maRe{βγ̅}mαaRe{βγ̅}).
Substituting (3.34) and (3.38) into (3.36) and noticing that (3.40)(iω0τ(j)I-∫-10eiω0τ(j)θdη(θ))q(0)=0,(-iω0τ(j)I-∫-10e-iω0τ(j)θdη(θ))q̅(0)=0,
we obtain (3.41)(2iω0τ(j)I-∫-10e2iω0τ(j)θdη(θ))E1=2τ(j)(-rK-β(rK+λ)βλ-mβγamαβγae-2iω0τ(j)).
This leads to(3.42)(2iω0+rKS*(rK+λ)S*0-λI*2iω0-mI*P*(a+I*)2dα0-amαP*(a+I*)2e-2iω0τ(j)2iω0+d-de-2iω0τ(j))E1=2(-rK-β(rK+λ)βλ-mβγamαβγae-2iω0τ(j)).
Solving this system for E1, we obtain (3.43)E1(1)=2A|-rK-β(rK+λ)(rK+λ)S*0βλ-mβγa2iω0-mI*P*(a+I*)2dαmαβγae-2iω0τ(j)-amαP*(a+I*)2e-2iω0τ(j)2iω0+d-de-2iω0τ(j)|,E1(2)=2A|2iω0+rKS*-rK-β(rK+λ)0-λI*βλ-mβγadα0mαβγae-2iω0τ(j)2iω0+d-de-2iω0τ(j)|,E1(3)=2A|2iω0+rKS*(rK+λ)S*-rK-β(rK+λ)-λI*2iω0-mI*P*(a+I*)2βλ-mβγa0-amαP*(a+I*)2e-2iω0τ(j)mαβγae-2iω0τ(j)|,
where(3.44)A=|2iω0+rKS*(rK+λ)S*0-λI*2iω0-mI*P*(a+I*)2dα0-amαP*(a+I*)2e-2iω0τ(j)2iω0+d-de-2iω0τ(j)|.
Similarly, substituting (3.35) and (3.39) into (3.37), we get(3.45)(rKS*(rK+λ)S*0-λI*-mI*P*(a+I*)2dα0-amαP*(a+I*)20)E2=2(-rK-(rK+λ)Re{β}λRe{β}-maRe{βγ̅}mαaRe{βγ̅}),
and hence(3.46)E2(1)=2B|-rK-(rK+λ)Re{β}(rK+λ)S*0λRe{β}-maRe{βγ̅}-mI*P*(a+I*)2dαmαaRe{βγ̅}-amαP*(a+I*)20|,E2(2)=2B|rKS*-rK-(rK+λ)Re{β}0-λI*λRe{β}-maRe{βγ̅}dα0mαaRe{βγ̅}0|,E2(3)=2B|rKS*(rK+λ)S*-rK-(rK+λ)Re{β}-λI*-mI*P*(a+I*)2λRe{β}-maRe{βγ̅}0-amαP*(a+I*)2mαaRe{βγ̅}|,
where(3.47)B=|rKS*(rK+λ)S*0-λI*-mI*P*(a+I*)2dα0-amαP*(a+I*)20|.Thus, we can determineW20(θ) and W11(θ) from (3.34) and (3.35). Furthermore, g21 in (3.26) can be expressed by the parameters and delay. Thus, we can compute the following values:
(3.48)c1(0)=i2ω0τ(j)(g20g11-2|g11|2-|g02|23)+g212,ν2=-Re{c1(0)}Re{ξ′(τ(j))},β2=2Re{c1(0)},T2=-Im⁡{c1(0)}+ν2Im⁡{ξ′(τ(j))}ω0τ(j),
which determine the qualities of bifurcating periodic solution in the center manifold at the critical value τ(j).Theorem 3.1.
ν2  determines the direction of the Hopf bifurcation. If ν2>0, then the Hopf bifurcation is supercritical and the bifurcating periodic solutions exist for τ>τ(j). If ν2<0, then the Hopf bifurcation is subcritical and the bifurcating periodic solutions exist for τ<τ(j). β2 determines the stability of the bifurcating periodic solutions: the bifurcating periodic solutions are stable if β2<0  and unstable if β2>0. T2  determines the period of the bifurcating periodic solutions: the period increase if T2>0 and decrease if T2<0.


## 4. 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:(4.1)r=3,K=40,λ=0.03,m=0.45,a=15,μ=0.28,α=0.42,d=0.09.
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)).Figure 1
Stable coexistence of prey and predator forτ=0. Parameter values are given in the text.Time series solutions of the prey and predator populations of the system (1.2): (a) τ=1, (b) τ=25. Parameter values are given in the text. This figure shows that the coexisting equilibrium E* is absolutely stable for all delay.
(a)(b)Behavior of the system (1.2) for different τ: (a) τ=1, (b) τ=3. All parameters are as in Figure 2 except m=0.72. This figure represents the conditional stability of the coexisting equilibrium E*.
(a)(b)For the given parameter set withm=0.72, one can evaluate that τ0=2.3187 and g′(hk)=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.Figure 4
Bifurcation diagram of the susceptible prey, infected prey and predator populations with respect to the delayτ. Parameters are as in Figure 3(b). This figure shows that the coexisting equilibrium E* is stable if τ<τ0=2.3187 and unstable if τ>τ0.Using Theorem3.1, one can determine the values of v2,β2 and T2. For the given parameter set with m=0.72, one can evaluate that v2=58.4107 (>0), β2=-1.3298 (<0), and T2=1.9364 (>0). Since v2>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 T2>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 τ.


## 5. 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.



---

*Source: 978387-2011-11-02.xml*