# Mathematical Analysis of a Single-Species Population Model in a Polluted Environment with Discrete Time Delays

Authors: Swarnali Sharma; G. P. Samanta
Journal: Journal of Mathematics
 (2013)

---

## Abstract

We have discussed the dynamical behaviour of a single-species population model in a polluted environment which describes the effect of toxicants on ecological system. Boundedness, positivity, and stability analysis of the model at various equilibrium points is discussed thoroughly. We have also studied the effect of single discrete delay as well as double discrete delays on the population model. Existence conditions of the Hopf bifurcation for single time delay are investigated. The length of delay preserving the stability is also estimated. The direction and the stability criteria of the bifurcating periodic solutions are determined by using the normal form theory and the center manifold theorem. The stability of the model with double time delays is investigated by using the Nyquist criteria. By choosing one of the delays as a bifurcation parameter, the model is found to undergo a Hopf bifurcation. Some numerical simulations for justifying the theoretical results are also illustrated by using MATLAB, which shows the reliability of our model from the practical point of view.


---

## Body

## 1. Introduction

In the world today, the pollution of the environment is a threatening problem due to the rapid development of industrialization. The presence of toxicants in the environment decreases the growth rate of species and its carrying capacity. In recent years, the most serious problem to society is the change in the environment caused by pollution, affecting the long-term survival of species, human life style, and biodiversity of the habitat. A large quantity of the toxicants and contaminants enter into the ecosystem one after another which seriously threaten the survival of the exposed population including humans. Therefore, the study of the effects of toxicants on the population and the assessment of the risk to populations are becoming quite important.The problem of estimating qualitatively the effect of a toxicant on a population by mathematical models is a very effective way. The study of deterministic dynamic population models with toxicant effect was proposed by Hallam and his colleagues in 1980s [1–3]. This model was revisited by many researchers [4, 5]. He and Ma in [6, 7] discussed the survival of a single-species in a polluted environment, considering the organism’s uptake of toxicant from the environment and egestion of toxicant into the environment. In [4] Ma studied a Leslie resource-consumer model and obtained the threshold between the persistence and extinction of the consumer. In [8], Buonomo et al. studied the effect of variation of the population on the toxicant concentration in the organism and environment. Samanta and Maiti [9] studied a dynamical model of a single-species system in a polluted environment under two cases: constant exogenous input of toxicant and rapidly fluctuating random exogenous input of toxicant into the environment by means of ordinary and stochastic differential equations. There are many other works on the effects of a single toxicant on various ecosystems by using mathematical models [10–15].Time delays have been incorporated in ecological models by many authors [16–27]. Introducing such delays in ecological models makes them more realistic and reasonable. Delay may have very complicated impact on the dynamics of a system since delay can cause the change of stability and can bifurcate periodic solutions.In this paper, we have developed a single-species population model in polluted environment. In Section2, we present a brief sketch of the construction of the model. Boundedness and positivity analysis of the solutions is shown which implies that the system is ecologically well behaved. In the next section, we have discussed the existence and stability analysis of various equilibrium points under zero-exogenous input and nonzero constant exogenous input. Next, we obtain necessary conditions for the existence of interior equilibrium E* and local and global stability of the system at E*. The analysis of single delayed model is described in Section 4. This analysis shows that there is a critical value τ* of delay τ below which the system is stable and above which the system becomes unstable at the interior equilibrium E*. The system undergoes a Hopf bifurcation around E* at τ=τ*. In Section 5, we have estimated the length of delay to preserve the stability around E*. The direction and the stability criteria of the bifurcation periodic solutions by using the normal form theory and the center manifold theorem are discussed in Section 6. In the next section, the analysis of double delayed model is discussed. The stability of the model with double delays is investigated by the Nyquist criteria. By choosing one of the delays as a bifurcation parameter, the model is found to undergo a Hopf bifurcation. Some of the important analytic results are numerically verified by using MATLAB in Section 8. Finally, Section 9 contains the general discussions of the paper and ecological implications of our mathematical findings.


## 2. Basic Mathematical Model

In this paper, we analyze a model which describes the effect of toxicant on a single-species population. Here we assume the following:n(t):
concentration of the population biomass,c(t):
concentration of the toxicant in the population,s(t):
concentration of the toxicant in the environment.The model satisfies the following assumptions.(A1)
There is a given toxicant in the environment and the living organisms absorb into their bodies part of this toxicant so that the dynamics of the population is affected by the toxicant.(A2)
For the growth rate of population, we assume that the birth rate isb0-fn(t) and the death rate is d0+αc(t), where b0, f, d0, and α are assumed to be positive constants.We consider the model(1)dndt=n(b0-d0-αc-fn),dcdt=ks-(r+m+b0-fn)c,dsdt=-ksn+(r+d0+αc)cn-hs+u(t),
with initial data n(0)>0, c(0)>0, and s(0)>0. Here,k: depletion rate of toxicant in the environment due to its intake made by the population,r: depletion rate of toxicant in the population due to egestion,m: depletion rate of toxicant in the population due to metabolization process,h: depletion rate of toxicant in the environment,u(t): exogenous toxicant input rate which is assumed to be a smooth bounded nonnegative function of t.We can see that ifb0-d0-αc(t)≤0, then n(t) will be going to extinct. So we suppose that
(2)p=b0-d0>0,b0-d0-αc(t)>0impliesc(t)<b0-d0α,∀t≥0.
The model we have just specified has nine parameters, which make the analysis difficult. To reduce the number of parameters, we make the following transformations: X=fn, Y=αc, and Z=kαs. So, system (1) becomes
(3)dXdt=X(p-Y-X),dYdt=Z-(q-X)Y,dZdt=-aXZ+a(d+Y)XY-hZ+γ(t),(4)withinitialdata:X(0)>0,Y(0)>0,Z(0)>0,
where p=b0-d0, q=r+m+b0, a=k/f, d=r+d0, and γ(t)=kαu(t).Therefore,  from (2), we  have
(5)Y(t)<p,∀t≥0[∵Y(t)=αc(t)].Theorem 1.
Each component of the solution of system (3) subject to initial conditions (4) is positive and bounded.Proof.
From the first equation of the system (3) we get
(6)X(t)=X(0)exp{∫0t‍p-Y(s)-X(s)ds}>0.
From the second equation of (3), we get
(7)dYdt≥-(q-X)Y⟹Y(t)≥Y(0)exp{-∫0t‍(q-X(s))ds}>0.
Similarly, from the third equation of (3), we get
(8)dZdt≥-(aX+h)Z⟹Z(t)≥Z(0)exp{-∫0t‍(aX(s)+h)ds}>0.
Therefore,X(t)>0, Y(t)>0, and Z(t)>0, for all t≥0.
From the first equation of (3), we have
(9)dXdt=X(p-Y-X)≤pX(1-Xp).
By a standard comparison theorem and by (5), we have,
(10)limsupt→∞X(t)≤p,Y(t)<p,∀t≥0.
From the third equation of (3), we have
(11)dZdt≤a(d+p)p2-hZ+γM⟹limsupt→∞Z(t)≤a(d+p)p2+γMh;γM=limsupt→∞γ(t).
So, we conclude that each component of the solution of system (3) subject to (4) is positive and bounded for all t≥0. This completes the proof.


## 3. Stability Behaviour of the Model

### 3.1. Case I: Zero Exogenous Input(γ(t)=0)

Theorem 2.
Ifγ(t)=0, then the model (3) has nonnegative equilibria E0(0,0,0) which is unstable and E1(p,0,0) which is locally asymptotically stable. The interior equilibrium E*(X*,Y*,Z*) is not feasible.Proof.
The variational matrix of system (3) at E0 is
(12)V(E0)=[p000-q100-h].
The eigen values arep, -q, and -h. So, it is obvious that, E0 is unstable (hyperbolic saddle).
Now the varitiational matrix of system (3) at E1 is
(13)V(E1)=[-p-p00-(q-p)10adp-(ap+h)].
The characteristic equation ofV(E1) is  (p+λ)(λ2+Bλ+C)=0, where B=ap+h+q-p=ap+h+d0+r+m>0 and C=(ap+h)(q-p)-apd=ap(q-p-d)+h(q-p)=apm+h(d0+r+m)>0, since p>0.
The eigen values areλ1=-p<0 and λ2,3=(-B±B2-4C)/2.
SinceB>0, and C>0, therefore the signs of the real parts of λ2, λ3 are negative. Hence E1 is locally asymptotically stable.
It is noted here that the other equilibrium point(X*,Y*,andZ*) is not feasible, since
(14)X*=hqh-am,Y*=-h(r+m+d0)-apmh-am
are opposite in sign since p>0, by (2).
This completes the proof.


### 3.2. Case II: Nonzero Exogenous Input(γ(t)=Q>0)

When(γ(t)=Q>0,aconstant), the model (3) has two nonnegative equilibria, E2(0,Q/hq,Q/h) and E*(X*,Y*,Z*). The variational matrix of system (3) at E2 is given by
(15)V(E2)=[p-Qhq00Qhq-q1-aQh+aQhq(d+Qhq)0-h].
The characteristic equation of V(E2) is (p-(Q/hq)-λ)(q+λ)(h+λ)=0.So,E2 is locally asymptotically stable if and only if hpq<Q and if hpq>Q, and then E2 becomes unstable.The interior equilibrium pointE*(X*,Y*,Z*) of system (3) is given by
(16)X*=A-hq+p(am-h)2(am-h),Y*=p(am-h)+hq-A2(am-h),Z*={p(am-h)+hq-A}{(am-h)(2q-p)+hq-A}4(am-h)2,
where
(17)A={hq-p(am-h)}2+4(am-h)(hpq-Q),m=(q-p-d).
It can be provided that the unique interior equilibrium point E*(X*,Y*,Z*) of system (3) exists if and only if the following two conditions are satisfied:
(18)(i)am>h,(ii)hpq>Q.
Summarizing the above analysis we come to the following theorem.Theorem 3.
Ifγ(t)=Q>0, then the equilibrium point E2(0,Q/hq,Q/h) of the system (3) is locally asymptotically stable if and only if hpq<Q and the unique interior equilibrium point E*(X*,Y*,Z*) of system (3) exists if and only if the following two conditions are satisfied:
(19)(i)am>h,(ii)hpq>Q,m=(q-p-d).


### 3.3. Local Stability ofE*(X*,Y*,Z*)

The variational matrix atE* is given by
(20)V(E*)=[m11m120m21m22m23m31m32m33],
where
(21)m11=-X*,m12=-X*,m21=Y*,m22=-q+X*,m23=1,m31=-aZ*+aY*(d+Y*),m32=aX*(d+2Y*),m33=-(aX*+h).
The characteristic equation is
(22)λ3+A1λ2+A2λ+A3=0,
where
(23)A1=-m11-m22-m33=-tr[V(E*)],A2=m11m22+m11m33+m22m33-m23m32-m12m21,A3=-det[V(E*)]A3=m11m23m32+m12m21m33A3=-m11m22m33-m12m23m31.
By the Routh-Hurwitz criterion [28–30], it follows that all eigenvalues of characteristic equation have negative real part if and only if
(24)A1>0,A3>0,A1A2-A3>0.Therefore, we come to the following theorem.Theorem 4.
E
* is locally asymptotically stable if and only if the inequalities (24) are satisfied.Furthermore, the conditions arising through Hopf bifurcation cannot be derived explicitly in terms of system parameters. Let us assume thatμ stands for any parameter involved with the model system (3). If there exists a critical magnitude μ=μ* such that A1(μ*)>0, A3(μ*)>0, (A1A2-A3)μ=μ*=0, and (d/dμ)(A1A2-A3)μ=μ*≠0, then E* undergoes a Hopf bifurcation at μ=μ* (Liu’s criterion) [28–30].


### 3.4.  Global Stability ofE*

E
* is not always globally asymptotically stable. The conditions which guarantee the global stability of E* are stated in the following theorem.Theorem 5.
SinceX(t), Y(t), and Z(t) are bounded, there exist positive constants mi,Mi (i=1,2,3) such that m1≤X(t)≤M1, m2≤Y(t)≤M2, and m3≤Z(t)≤M3. If the following inequalities hold:
(25)(M2-1)2<(q-X*),{aM2(d+M2)-am3}2<(aX*+h),(1+adX*+2aX*M2)2<(q-X*)(h+aX*),
then E* is globally asymptotically stable.Proof.
We consider the following positive definite function which has a strict minimum atE*:
(26)L(X,Y,Z)=(X-X*-X*lnXX*)+12(Y-Y*)2+12(Z-Z*)2.
Differentiating both sides with respect to t along the solution of (3), we get (after some simple calculations)
(27)dLdt=-(X-X*)2+(Y-1)(X-X*)(Y-Y*)dLdt=-(q-X*)(Y-Y*)2dLdt=+(aY(d+Y)-aZ)(X-X*)(Z-Z*)dLdt=-(aX*+h)(Z-Z*)2dLdt=+(1+adX*+aX*(Y+Y*))(Y-Y*)(Z-Z*)dLdt≤-(X-X*)2+(M2-1)(X-X*)(Y-Y*)dLdt=-(q-X*)(Y-Y*)2-(aX*+h)(Z-Z*)2dLdt=+(aM2(d+M2)-am3)(X-X*)(Z-Z*)dLdt=+(1+adX*+2aM2X*)(Y-Y*)(Z-Z*)dLdt=-a11(X-X*)2-a22(Y-Y*)2-a33(Z-Z*)2dLdt=+a12(X-X*)(Y-Y*)+a13(X-X*)(Z-Z*)dLdt=+a23(Y-Y*)(Z-Z*),
where
(28)a11=1,a22=q-X*,a33=aX*+h,a12=M2-1,a13={aM2(d+M2)-am3},a23=1+adX*+2aM2X*.
Now, sufficient conditions for dL/dt to be negative definite are
(29)a122-a11a22<0,a132-a11a33<0,a232-a22a33<0,
that is,
(30)(M2-1)2<(q-X*),{aM2(d+M2)-am3}2<(aX*+h),(1+adX*+2aX*M2)2<(q-X*)(h+aX*).
HenceL is a Lyapunov function and so E* is globally asymptotically stable [28–30]. This competes the proof.


## 3.1. Case I: Zero Exogenous Input(γ(t)=0)

Theorem 2.
Ifγ(t)=0, then the model (3) has nonnegative equilibria E0(0,0,0) which is unstable and E1(p,0,0) which is locally asymptotically stable. The interior equilibrium E*(X*,Y*,Z*) is not feasible.Proof.
The variational matrix of system (3) at E0 is
(12)V(E0)=[p000-q100-h].
The eigen values arep, -q, and -h. So, it is obvious that, E0 is unstable (hyperbolic saddle).
Now the varitiational matrix of system (3) at E1 is
(13)V(E1)=[-p-p00-(q-p)10adp-(ap+h)].
The characteristic equation ofV(E1) is  (p+λ)(λ2+Bλ+C)=0, where B=ap+h+q-p=ap+h+d0+r+m>0 and C=(ap+h)(q-p)-apd=ap(q-p-d)+h(q-p)=apm+h(d0+r+m)>0, since p>0.
The eigen values areλ1=-p<0 and λ2,3=(-B±B2-4C)/2.
SinceB>0, and C>0, therefore the signs of the real parts of λ2, λ3 are negative. Hence E1 is locally asymptotically stable.
It is noted here that the other equilibrium point(X*,Y*,andZ*) is not feasible, since
(14)X*=hqh-am,Y*=-h(r+m+d0)-apmh-am
are opposite in sign since p>0, by (2).
This completes the proof.


## 3.2. Case II: Nonzero Exogenous Input(γ(t)=Q>0)

When(γ(t)=Q>0,aconstant), the model (3) has two nonnegative equilibria, E2(0,Q/hq,Q/h) and E*(X*,Y*,Z*). The variational matrix of system (3) at E2 is given by
(15)V(E2)=[p-Qhq00Qhq-q1-aQh+aQhq(d+Qhq)0-h].
The characteristic equation of V(E2) is (p-(Q/hq)-λ)(q+λ)(h+λ)=0.So,E2 is locally asymptotically stable if and only if hpq<Q and if hpq>Q, and then E2 becomes unstable.The interior equilibrium pointE*(X*,Y*,Z*) of system (3) is given by
(16)X*=A-hq+p(am-h)2(am-h),Y*=p(am-h)+hq-A2(am-h),Z*={p(am-h)+hq-A}{(am-h)(2q-p)+hq-A}4(am-h)2,
where
(17)A={hq-p(am-h)}2+4(am-h)(hpq-Q),m=(q-p-d).
It can be provided that the unique interior equilibrium point E*(X*,Y*,Z*) of system (3) exists if and only if the following two conditions are satisfied:
(18)(i)am>h,(ii)hpq>Q.
Summarizing the above analysis we come to the following theorem.Theorem 3.
Ifγ(t)=Q>0, then the equilibrium point E2(0,Q/hq,Q/h) of the system (3) is locally asymptotically stable if and only if hpq<Q and the unique interior equilibrium point E*(X*,Y*,Z*) of system (3) exists if and only if the following two conditions are satisfied:
(19)(i)am>h,(ii)hpq>Q,m=(q-p-d).


## 3.3. Local Stability ofE*(X*,Y*,Z*)

The variational matrix atE* is given by
(20)V(E*)=[m11m120m21m22m23m31m32m33],
where
(21)m11=-X*,m12=-X*,m21=Y*,m22=-q+X*,m23=1,m31=-aZ*+aY*(d+Y*),m32=aX*(d+2Y*),m33=-(aX*+h).
The characteristic equation is
(22)λ3+A1λ2+A2λ+A3=0,
where
(23)A1=-m11-m22-m33=-tr[V(E*)],A2=m11m22+m11m33+m22m33-m23m32-m12m21,A3=-det[V(E*)]A3=m11m23m32+m12m21m33A3=-m11m22m33-m12m23m31.
By the Routh-Hurwitz criterion [28–30], it follows that all eigenvalues of characteristic equation have negative real part if and only if
(24)A1>0,A3>0,A1A2-A3>0.Therefore, we come to the following theorem.Theorem 4.
E
* is locally asymptotically stable if and only if the inequalities (24) are satisfied.Furthermore, the conditions arising through Hopf bifurcation cannot be derived explicitly in terms of system parameters. Let us assume thatμ stands for any parameter involved with the model system (3). If there exists a critical magnitude μ=μ* such that A1(μ*)>0, A3(μ*)>0, (A1A2-A3)μ=μ*=0, and (d/dμ)(A1A2-A3)μ=μ*≠0, then E* undergoes a Hopf bifurcation at μ=μ* (Liu’s criterion) [28–30].


## 3.4.  Global Stability ofE*

E
* is not always globally asymptotically stable. The conditions which guarantee the global stability of E* are stated in the following theorem.Theorem 5.
SinceX(t), Y(t), and Z(t) are bounded, there exist positive constants mi,Mi (i=1,2,3) such that m1≤X(t)≤M1, m2≤Y(t)≤M2, and m3≤Z(t)≤M3. If the following inequalities hold:
(25)(M2-1)2<(q-X*),{aM2(d+M2)-am3}2<(aX*+h),(1+adX*+2aX*M2)2<(q-X*)(h+aX*),
then E* is globally asymptotically stable.Proof.
We consider the following positive definite function which has a strict minimum atE*:
(26)L(X,Y,Z)=(X-X*-X*lnXX*)+12(Y-Y*)2+12(Z-Z*)2.
Differentiating both sides with respect to t along the solution of (3), we get (after some simple calculations)
(27)dLdt=-(X-X*)2+(Y-1)(X-X*)(Y-Y*)dLdt=-(q-X*)(Y-Y*)2dLdt=+(aY(d+Y)-aZ)(X-X*)(Z-Z*)dLdt=-(aX*+h)(Z-Z*)2dLdt=+(1+adX*+aX*(Y+Y*))(Y-Y*)(Z-Z*)dLdt≤-(X-X*)2+(M2-1)(X-X*)(Y-Y*)dLdt=-(q-X*)(Y-Y*)2-(aX*+h)(Z-Z*)2dLdt=+(aM2(d+M2)-am3)(X-X*)(Z-Z*)dLdt=+(1+adX*+2aM2X*)(Y-Y*)(Z-Z*)dLdt=-a11(X-X*)2-a22(Y-Y*)2-a33(Z-Z*)2dLdt=+a12(X-X*)(Y-Y*)+a13(X-X*)(Z-Z*)dLdt=+a23(Y-Y*)(Z-Z*),
where
(28)a11=1,a22=q-X*,a33=aX*+h,a12=M2-1,a13={aM2(d+M2)-am3},a23=1+adX*+2aM2X*.
Now, sufficient conditions for dL/dt to be negative definite are
(29)a122-a11a22<0,a132-a11a33<0,a232-a22a33<0,
that is,
(30)(M2-1)2<(q-X*),{aM2(d+M2)-am3}2<(aX*+h),(1+adX*+2aX*M2)2<(q-X*)(h+aX*).
HenceL is a Lyapunov function and so E* is globally asymptotically stable [28–30]. This competes the proof.


## 4. 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 environmental toxicant on the population growth will not be instantaneous, but mediated by some discrete time lagτ required for incubation. So, we modify system (3) and (4) as follows:
(31)dXdt=X(p-Y-X)dYdt=Z(t-τ)-(q-X)YdZdt=-aXZ+a(d+Y)XY-hZ+Q
with initial conditions
(32)X(0)>0,Y(0)>0,Z(ϑ)=φ(ϑ),suchthatφ(ϑ)≥0,∀ϑ∈[-τ,0],
where φ(ϑ) is a nonnegative continuous functions on ϑ∈[-τ,0]. For a biological meaning, we further assume that φ(0)>0.All parameters are the same as in system (3) except that the positive constant τ represents the activation period or reaction time of the toxicant in the population and γ(t)=Q(>0) is a constant.The system (31) has the same equilibria as in art.3.2. The main purpose of this section is to study the stability behaviour of E*(X*,Y*,andZ*) in the presence of discrete delay (τ>0). We linearize system (31) by using the following transformations:
(33)X(t)=X*+u(t),Y(t)=Y*+v(t),Z(t)=Z*+w(t).
Then transformed linear system is given by
(34)dUdt=AU(t)+BU(t-τ),
where
(35)U(t)=[uvw]T,A=(aij)3×3,B=(bij)3×3,a11=-X*,a12=-X*,a21=Y*,a22=-q+X*,a31=-aZ*+adY*+a(Y*)2,a32=adX*+2aX*Y*,a33=-(aX*+h)andallotheraij=0;b23=1andallotherbij=0.
Let us take the solution of the system (34) as u(t)=ρ1eλt, v(t)=ρ2eλt, w(t)=ρ3eλt. This leads to the following characteristic equation:
(36)λ3+a1λ2+a2λ+a3+(a4λ+a5)e-λτ=0,
where
(37)a1=-a11-a22-a33,a2=a11a22+a11a33+a22a33-a12a21,a3=a12a21a33-a11a22a33,a4=-a32b23,a5=a11a32b23-a12a31b23.A necessary condition for a stability change ofE* is that the characteristic equation  should have purely imaginary solutions. In order to see the delay induced instability, let us consider τ as bifurcation parameter and assume a purely imaginary solution of (3) in the form λ=iω(ω≠0). Therefore, substituting λ=iω in (37) and separating real and imaginary parts, we get
(38)a1ω2-a3=a5cos(ωτ)+a4ωsin(ωτ),-ω3+a2ω=a5sin(ωτ)-a4ωcos(ωτ).
Eliminating τ by squaring and adding (38), we get the equation for determining ω as
(39)ω6+d1ω4+d2ω2+d3=0,
where
(40)d1=a12-2a2,d2=a22-2a1a3-a42,d3=a32-a52.
Substituting ω2=σ in (40), we get a cubic equation given by
(41)Φ(σ)≡σ3+d1σ2+d2σ+d3=0.
Now, from (38), we get
(42)τk=1ω0arccos(a4ω04+(a1a5-a2a4)ω02-a3a5a52+a4ω02)+2kπω0,k=0,1,2,….
Thus when τ=τk, the characteristic equation (37) has a pair of purely imaginary roots ±iω0. Let λ(τ)=ξ(τ)+iω(τ) be a root of (37) such that the following two conditions hold:
(43)ξ(τk)=0,ω(τk)=ω0.
Differentiating both sides of (37) with respect to τ, we get
(44)(dλdτ)-1=3λ2+2a1λ+a2(a4λ2+a5λ)e-λτ+a4a4λ2+a5λ-τλ.∴Re(dλdτ)λ=iω0-1=1ω0(a42ω02+a52)×[(a2a5-3ω02a5+2a1a4ω02)sin(ω0τ)+(3a4ω03+2a1a5ω0-a2a4ω0)cos(ω0τ)-a42ω0]=[dϑ(σ)/dσ]σ=ω02(a42ω02+a52).
Now, if σ=ω02 is the first positive root of (41) and hence σ=ω0(>0) is the first positive root of (40), and then
(45)Re(dλdτ)λ=iω0-1>0as[dϑ(σ)dσ]σ=ω02>0.
Thus, the transversality condition is satisfied and the steady state becomes unstable when τ>τ*. Moreover, a Hopf bifurcation [28–30] occurs when τ passes through the critical value τ*, where τ* is the smallest positive value of τk given in (42).We notice that, in (41), Φ is continuous everywhere with Φ(0)<0 and Φ(σ)→∞ as σ→∞. Therefore, the cubic equation (41) always has at least one positive root. Consequently the stability criteria of the system for τ=0 will not necessarily ensure the stability of the system for τ≠0. In the following theorem, we have given a criterion for switching the stability behaviour of E*.Theorem 6 (see [28–30]).
IfE* exists with the conditions (24) and σ0=ω02 be a positive root of (41), then there exists a τ=τ* such that(i)
E* is locally asymptotically stable for 0≤τ<τ*,(ii)
E* is unstable for τ>τ*,(iii)
the system (31) undergoes a Hopf-bifurcation around E* at τ=τ*=ming(ω0),
(46)whereg(ω0)=1ω0arccos×(a4ω04+(a1a5-a2a4)ω02-a3a5a52+a4ω02)and the minimum is taken over all positive ω0 such that ω02 is a solution of (41).


## 5. Estimation of the Length of Delay to Preserve Stability

We linearize the system (31) about its interior equilibrium E*(X*,Y*,andZ*) and get
(47)dx1dt=a11x1+a12y1dy1dt=a21x1+a22y1+b23z1(t-τ)dz1dt=a31x1+a32y1+a33z1,
where
(48)X(t)=x1(t)+X*,Y(t)=y1(t)+Y*,Z(t)=z1(t)+Z*.
Taking Laplace transform of the system (47), we get
(49)(s-a11)x1¯(s)=a12y1¯(s)+x1(0)(s-a22)y1¯(s)=a21x1¯(s)+b23e-sτz1¯(s)+b23e-sτk1(s)+y1(0)(s-a33)z1¯(s)=a31x1¯(s)+a32y1¯(s)+z1(0),
where
(50)k1(s)=∫-τ0‍e-stz1(t)dt
and x1¯(s), y1¯(s), and z1¯(s) are the Laplace transform of x1(t), y1(t) and z1(t), respectively.Now, we will use the “Nyquist theorem” [21] which states that if s is the arc length along a curve encircling the right half of the plane, then a curve Γ(s) will encircle the origin a number of times equal to the difference between the number of poles and the number of zeros of Γ(s) in the right half of the plane. Using “Nyquist theorem” [21, 31], it can be shown that the conditions for local asymptotic stability of E*(X*,Y*,Z*) are given by [32]
(51)ImH(iω0)>0,(52)ReH(iω0)=0,
where H(s)=s3+a1s2+a2s+a3+(a4s+a5)e-sτ and ω0 is the smallest positive root of (52).We have already shown thatE*(X*,Y*,Z*) is stable in absence of delay (provided that inequalities given in (24) are satisfied). Hence, by continuity, all eigenvalues will continue to have negative real parts for sufficiently small τ>0 provided that one can guarantee that no eigenvalue with positive real part bifurcates from infinity as τ increases from zero. This can be proved by using Butler’s lemma [32].In this case, conditions (51) and (52) give
(53)a2ω0-ω03>a5sin(ω0τ)-a4ω0cos(ω0τ),(54)a3-a1ω02=-a5cos(ω0τ)-a4ω0sin(ω0τ).Now, (53) and (54), if satisfied simultaneously, are sufficient conditions to guarantee stability. We will utilize them to get an estimate on the length of delay. Our aim is to find an upper bound ω+ on ω0, independent of τ so that (53) holds for all values of ω, 0≤ω≤ω+ and hence in particular at ω=ω0.We rewrite (54) as
(55)a1ω02=a3+a5cos(ω0τ)+a4ω0sin(ω0τ).Maximizinga3+a5cos(ω0τ)+a4ω0sin(ω0τ) subject to |sin(ω0τ)|≤1, |cos(ω0τ)|≤1, we obtain
(56)a1ω02≤a3+|a5|+|a4|ω0.Hence, if(57)ω+=12a1[|a4|+a42+4a1(a3+|a5|)],
then clearly from (56) and (57), we have ω0≤ω+.From the inequality (53), we get
(58)ω02<a2+a4cos(ω0τ)-a5ω0sin(ω0τ).AsE*(X*,Y*,Z*) is locally asymptotically stable for τ=0 (provided that inequalities given in (24) are satisfied), for sufficiently small τ>0, (58) will continue to hold. Substituting (55) into (58) and rearranging, we get
(59)(a1a4-a5){1-cos(ω0τ)}+(a4ω0+a1a5ω0)sin(ω0τ)<a1a2+a1a4-a3-a5.
Using the bounds we obtain
(60)(a1a4-a5){1-cos(ω0τ)}=(a1a4-a5)2sin2(ω0τ2)≤12|a1a4-a5|ω+2τ2,(a4ω0+a1a5ω0)sin(ω0τ)≤(a4ω+2+a1a5)τ.Now, from (59) and (60), we get
(61)l1τ2+l2τ<l3,
where
(62)l1=12|a1a4-a5|ω+2,l2=(a4ω+2+a1a5),l3=a1a2+a1a4-a3-a5.(63)Hence,ifτ+=12l1[-l2+l22+4l1l3],
then stability is preserved for 0≤τ<τ+.Summarizing the previous discussions we come to the following theorem.Theorem 7.
The delayed model (31) and (32) will be locally asymptotically stable at E*(X*,Y*,andZ*) with the conditions (24) if the delay τ lies within the interval (0,τ+) where τ+ is given by (63).


## 6. Direction of Hopf Bifurcation

In this section, we will derive the direction of the Hopf bifurcation [17, 21, 33] and sufficient conditions of the stability of bifurcating periodic solution from the interior equilibrium E* of the system (31) at the critical value τ=τ*. We will employ the approach of the normal form method and center manifold theorem introduced by Hassard et al. [33].Letu1=X-X*, u2=Y-Y*, u3=Z-Z*, τ=τ*+μ, u(t)=(u1(t),u2(t),u3(t))T, and ut(θ)=u(t+θ) for θ∈[-τ*,0]. 𝒞=C([-τ*,0],ℝ3), and we can rewrite the system (31) and (32) as follows:
(64)u˙(t)=Lμ(ut)+F(ut,μ),
where Lμ:ℝ×𝒞→ℝ3, F:ℝ×𝒞→ℝ3 are given, respectively, by
(65)Lμ(ϕ)=B1ϕ(0)+B2ϕ(-τ*),F(ϕ,μ)=[-ϕ12(0)-ϕ1(0)ϕ2(0)ϕ1(0)ϕ2(0)aX*ϕ22(0)+(ad+2aY*)ϕ1(0)ϕ2(0)-aϕ1(0)ϕ3(0)],
where
(66)B1=(a11a120a21a220a31a32a33),B2=(00000b23000),a11=-X*,a12=-X*,a21=Y*,a22=-q+X*,a31=-aZ*+adY*+a(Y*)2,a32=adX*+2aX*Y*,a33=-(aX*+h)andallotheraij=0;b23=1andallotherbij=0.
By the Riesz representation theorem, there exists a matrix function whose components are bounded variation function η(θ,μ) in θ∈[-τ*,0] such that
(67)Lμ(θ)=∫-τ*0‍dη(θ,μ)ϕ(θ),ϕ∈𝒞.
In fact we can choose
(68)η(θ,μ)=B1δ(θ)+B2δ(θ+τ*),
where δ(θ) is the Dirac delta function.Forϕ∈𝒞1([-τ*,0],ℝ3), we define
(69)A(μ)ϕ(θ)={dϕ(θ)dθ,θ∈[-τ*,0)∫-τ*0dη(θ,μ)ϕ(θ),θ=0,(70)R(μ)ϕ(θ)={0,θ∈[-τ*,0)F(μ,ϕ),θ=0.
Then the system (64) is equivalent to
(71)u˙t=A(μ)ut+R(μ)ut.
For ψ∈𝒞1([0,τ*],(ℝ3)*), where (ℝ3)* is the 3-dimensional space of row vectors, the adjoint operator A* of A is defined as
(72)A*(μ)ψ(s)={-dψ(s)ds,s∈(0,τ*]∫-τ*0‍dηT(s,μ)ψ(-s),s=0.
For ϕ∈[-τ*,0) and ψ∈[0,τ*], we define a bilinear inner product taking the first argument to be conjugate linear and the second to be linear of the form
(73)〈ψ,ϕ〉=ψ¯(0)ϕ(0)-∫θ=-τ*0‍∫ξ=0θ‍ψ¯(ξ-θ)dη(θ)ϕ(ξ)dξ,
where η(θ)=η(θ,0). It is easy to verify that A(0) and A* are a pair of adjoint operators.From the discussions in the previous section, we know that±iω0 are the eigenvalues of A(0). Hence, the eigenvalues of A* are ∓iω0.Suppose thatq(θ)=(1,α1,α2)Teiω0θ is eigenvector of A(0) corresponding to iω0, and then we have
(74)A(0)q(θ)=iω0q(θ).
It follows from (67), (68), and the definition of A(0) that
(75)[(B1+B2e-iω0τ*)-iω0I]q(0)=0,
where I is identity matrix of order 3×3, that is,
(76)(a11-iω0a120a21a22-iω0b23e-iω0τ*a31a32a33-iω0)q(0)=(000).
Solving (76), we have q(0)=(1,α1,α2)T, where
(77)α1=iω0-a11a12,α2=-a12a21+(a22-iω0)(iω0-a11)a12b23e-iω0τ*.
Now, suppose that q*(s)=B(1,α1*,α2*)eiω0s is the eigenvector of A* corresponding to -iω0, and then we have
(78)A*(0)q*(s)=-iω0q*(s).
It follows from (67), (68), and the definition of A* that,
(79)[(B1T+B2Teiω0τ*)+iω0I](q*(0))T=0,
where I is identity matrix of order 3×3, that is,
(80)(a11+iω0a21a31a12a22+iω0a320b23eiω0τ*a33+iω0)(q*(0))T=(000).
Solving (80), we have q*(0)=B(1,α1*,α2*), where
(81)α1*=-(a11+iω0)(a33+iω0)a21(a33+iω0)-a31b23eiω0τ*,α2*=(a11+iω0)b23eiω0τ*a21(a33+iω0)-a31b23eiω0τ*.
Now,
(82)〈q*(s),q(θ)〉=B¯{∫θ=-τ*0‍∫ξ=0θ‍(1,α¯1*,α¯2*)(1,α1,α2)T-∫θ=-τ*0‍∫ξ=0θ‍(1,α¯1*,α¯2*)e-i(ξ-θ)ω0dη(θ)×(1,α1,α2)Teiξω0dξ∫θ=-τ*0‍∫ξ=0θ‍}=B¯[∫θ=-τ*0‍1+α1α¯1*+α2α¯2*-∫θ=-τ*0‍(1,α¯1*,α¯2*)θeiω0θdη(θ)(1,α1,α2)T]=B¯[1+α1α¯1*+α2α¯2*+τ*e-iω0τ*(α2α¯1*b23)].
Thus in order to assume 〈q*(s),q(θ)〉=1, we can choose
(83)B=11+α1α¯1*+α2α¯2*+τ*e-iω0τ*(α2α¯1*).
Also since 〈ψ,Aϕ〉=〈A*ψ,ϕ〉, we have
(84)-iω0〈q*,q¯〉=〈q*,Aq¯〉,-iω0〈q*,q¯〉=〈A*q*,q¯〉,-iω0〈q*,q¯〉=〈-iω0q*,q¯〉,-iω0〈q*,q¯〉=iω0〈q*,q¯〉.
Therefore, 〈q*,q¯〉=0.Now, we first compute the coordinates to describe the center manifoldC0 at μ=0. Let ut be the solution of (64) when μ=0. Define
(85)z(t)=〈q*,ut〉,W(t,θ)=ut-z(t)q(θ)-z¯(t)q¯(θ)W(t,θ)=ut(θ)-2Re{z(t)q(θ)}.
On the center manifold C0, we have
(86)W(z,z¯,θ)=W20(θ)z22+W11(θ)zz¯+W02(θ)z¯22+⋯,
where z and z¯ are local coordinates of center manifold C0 in the direction of q* and q¯*. Note that W is real if ut is real. We only consider real solutions.From (85), we get
(87)〈q*,W〉=〈q*,ut-zq-z¯q¯〉〈q*,W〉=〈q*,ut〉-z〈q*,q〉-z¯〈q*,q¯〉=z-z-0=0.
For the solution ut∈C0 of (64) and since μ=0, we have
(88)z˙(t)=〈q*,u˙t〉z˙(t)=〈q*,A(0)ut+R(0)ut〉z˙(t)=〈A*(0)q*,ut〉+q¯*(0)F(0,ut)z˙(t)=〈-iω0q*,ut〉+q¯*(0)f0(z,z¯)z˙(t)=iω0z+q¯*(0)f0(z,z¯)z˙(t)=iω0z+g(z,z¯),
where
(89)g(z,z¯)=q¯*(0)f0(z,z¯)g(z,z¯)=g20z22+g11zz¯+g02z¯22+g21z2z¯2+⋯.
Since ut=W(z,z¯,θ)+zq+z¯q¯, we have
(90)(u1t(θ)u2t(θ)u3t(θ))=(W(1)(θ)W(2)(θ)W(3)(θ))+z(1α1α2)eiω0θ+z¯(1α¯1α¯2)e-iω0θ=(Γ1Γ2Γ3),
where
(91)Γ1=zeiω0θ+z¯e-iω0θ+W20(1)(θ)z22+W11(1)(θ)zz¯Γ1=+W02(1)(θ)z¯22+o(|(z,z¯)|3),Γ2=α1zeiω0θ+α¯1z¯e-iω0θ+W20(2)(θ)z22+W11(2)(θ)zz¯Γ2=+W02(2)(θ)z¯22+o(|z,z¯|3),Γ3=α2zeiω0θ+α¯2z¯e-iω0θ+W20(3)(θ)z22+W11(3)(θ)zz¯Γ3=+W02(3)(θ)z¯22+o(|z,z¯|3).Then,(92)u1t(0)=z+z¯+W20(1)(0)z22+W11(1)(0)zz¯u1t(0)=+W02(1)(0)z¯22+o(|z,z¯|3),u2t(0)=α1z+α¯1z¯+W20(2)(0)z22+W11(2)(0)zz¯u2t(0)=+W02(2)(0)z¯22+o(|z,z¯|3),u3t(0)=α2z+α¯2z¯+W20(3)(0)z22+W11(3)(0)zz¯u3t(0)=+W02(3)(0)z¯22+o(|z,z¯|3),u1t2(0)=z2+2zz¯+z¯2+[W20(1)(0)+2W11(1)(0)]z2z¯u1t2(0)=+h.o.t.,u2t2(0)=α12z2+2|α1|2zz¯+α¯12z¯2u2t2(0)=+[α¯1W20(2)(0)+2α1W11(2)(0)]z2z¯+h.o.t.,u1t(0)u2t(0)=α1z2+(α1+α¯1)zz¯+α¯1z¯2+[W11(1)(0)α1+12W20(2)(0)+12W20(1)(0)α¯1+W11(2)(0)12W20(2)(0)+12W20(1)]z2z¯+h.o.t.,u1t(0)u3t(0)=α2z2+(α2+α¯2)zz¯+α¯1α¯2z¯2+[W11(1)(0)α2+12W20(1)(0)α¯2+12W20(3)(0)+W11(3)(0)+12W20(1)(0)α¯2+12W20(3)]z2z¯+h.o.t.It follows that(93)f0(z,z¯)=[Δ11z2+Δ12zz¯+Δ13z¯2+Δ14z2z¯Δ21z2+Δ22zz¯+Δ23z¯2+Δ24z2z¯Δ31z2+Δ32zz¯+Δ33z¯2+Δ34z2z¯]+⋯,
where
(94)Δ11=-1-α1,Δ12=-2-(α1+α¯1),Δ13=-1-α¯1,Δ14=-[W20(1)(0)+2W11(1)(0)+12W20(1)(0)α¯1+12W20(2)(0)+W11(1)(0)α1+W11(2)(0)],Δ21=α1,Δ22=α1+α¯1,Δ23=α¯1,Δ24=12W20(1)(0)α¯1+12W20(2)(0)+W11(1)(0)α1+W11(2)(0),Δ31=aX*α12+(ad+2aY*)α1-aα2,Δ32=2aX*|α1|2+(ad+2aY*)(α1+α¯1)-a(α2+α¯2),Δ33=aX*α¯12+(ad+2aY*)α¯1-aα¯2,Δ34=aX*[α¯1W20(2)(0)+2α1W11(2)(0)]+(ad+2aY*)Δ34=×[W11(1)(0)α1+12W20(2)(0)+12W20(1)(0)α¯1+W11(2)(0)+12W20(2)(0)+12W20(1)(0)α¯1]Δ34=-a[W11(1)(0)α2+12W20(1)(0)α¯2+12W20(3)(0)+W11(3)(0)+12W20(1)(0)α¯2+12W20(3)(0)].
Since q¯*(0)=B¯(1,α¯1*,α¯2*), then we obtain
(95)g(z,z¯)=q¯*(0)f0(z,z¯)g(z,z¯)=B¯(1,α¯1*,α¯2*)g(z,z¯)=×[Δ11z2+Δ12zz¯+Δ13z¯2+Δ14z2z¯Δ21z2+Δ22zz¯+Δ23z¯2+Δ24z2z¯Δ31z2+Δ32zz¯+Δ33z¯2+Δ34z2z¯]+⋯g(z,z¯)=B¯[(Δ11+α¯1*Δ21+α¯2*Δ31)z2+(Δ12+α¯1*Δ22+α¯2*Δ32)zz¯+(Δ13+α¯1*Δ23+α¯2*Δ33)z¯2+(Δ14+α¯1*Δ24+α¯2*Δ34)z2z¯]+⋯.
The comparison of the coefficients of z2, zz¯, z¯2, z2z¯ from (89) and (95) produces
(96)g20=2B¯(Δ11+α¯1*Δ21+α¯2*Δ31)g11=B¯(Δ12+α¯1*Δ22+α¯2*Δ32)g02=2B¯(Δ13+α¯1*Δ23+α¯2*Δ33)g21=2B¯(Δ14+α¯1*Δ24+α¯2*Δ34).
Clearly, g21 involves the terms W20, W11, and we still need to compute W20, W11 as
(97)W˙=u˙t-z˙q-z¯˙q¯={A(0)W-2Re(q¯*(0)f0(z,z¯)q(θ)),θ∈[-τ*,0)A(0)W-2Re(q¯*(0)f0(z,z¯)q(θ))+f0(z,z¯),θ=0
which is rewritten as
(98)W˙=A(0)W+H(z,z¯,θ),
where
(99)H(z,z¯,θ)=H20(θ)z22+H11(θ)zz¯+H02z¯22+⋯.
On the other hand, on C0(100)W˙=Wzz˙+Wz¯z¯˙.
We derive the following inequalities from (98)–(100):
(101)(A(0)-2iω0)W20(θ)=-H20(θ)A(0)W11(θ)=-H11(θ).
For θ∈[-τ*,0), it follows from (89), (97), and (98) that
(102)H(z,z¯,θ)=-2Re(q¯*(0)f0(z,z¯)q(θ))=-2Re(g(z,z¯)q(θ))=-g(z,z¯)q(θ)-g¯(z,z¯)q¯(θ)=-(g20z22+g11zz¯+g02z¯22+g21z2z¯2+⋯)q(θ)-(g¯20z¯22+g¯11zz¯+g02z22+g¯21z¯2z2+⋯)q¯(θ).
Computing this equation with (99), we get
(103)H20(θ)=-g20q(θ)-g¯02q¯(θ),(104)H11(θ)=-g11q(θ)-g¯11q¯(θ).
For θ∈[-τ*,0), it follows from (69), (101), and (103) that
(105)W˙20(θ)=A(0)W20(θ)W˙20(θ)=2iω0W20(θ)-H20(θ)W˙20(θ)=2iω0W20(θ)+g20q(θ)+g¯02q¯(θ)W˙20(θ)=2iω0W20(θ)+g20q(0)eiω0θ+g¯02q¯(0)e-iω0θ.
Solving for W20(θ), we get
(106)W20(θ)=ig20q(0)ω0eiω0θ+ig¯02q¯(0)3ω0e-iω0θ+E^1e2iω0θ.
And similarly from (69), (101), and (104), it follows that
(107)W11(θ)=-ig11q(0)ω0eiω0θ+ig¯11q¯(0)ω0e-iω0θ+E^2,
where E^1=(E1(1),E1(2),E1(3))T and E^2=(E2(1),E2(2),E2(3))T are both three dimensional constant vectors and can be determined by setting θ=0 in H(z,z¯,θ). This leads to
(108)H(z,z¯,0)=-2Re(q¯*(0)f0(z,z¯)q(0))+f0(z,z¯)=-2Re(g(z,z¯)q(0))+f0(z,z¯)=-g(z,z¯)q(0)-g¯(z,z¯)q¯(0)+f0(z,z¯)=-(g20z22+g11zz¯+g02z¯22+g21z2z¯2+⋯)q(0)-(g¯20z¯22+g¯11zz¯+g02z22+g¯21z¯2z2+⋯)q¯(0)+[Δ11z2+Δ12zz¯+Δ13z¯2+Δ14z2z¯Δ21z2+Δ22zz¯+Δ23z¯2+Δ24z2z¯Δ31z2+Δ32zz¯+Δ33z¯2+Δ34z2z¯]+⋯.
Comparing the corresponding coefficients of the above equation with (99) for θ=0, we get
(109)H20(0)=-g20q(0)-g¯02q¯(0)+2(Δ11Δ21Δ31),(110)H11(0)=-g11q(0)-g¯11q¯(0)+2(Δ12Δ22Δ32).
From the definition of A and (101) for θ=0, we get
(111)∫-τ*0‍dη(θ)W20(θ)=2iω0W20(0)-H20(0),(112)∫-τ*0‍dη(θ)W11(θ)=-H11(0).
Notice that
(113)(iω0I-∫-τ*0‍eiω0θdη(θ))q(0)=0,(-iω0I-∫-τ*0‍e-iω0θdη(θ))q¯(0)=0.
Substituting (106) and (111) into (109), we obtain
(114)(2iω0I-∫-τ*0‍e2iω0θdη(θ))E^1=2(Δ11Δ21Δ31)
which leads to
(115)(2iω0I-B1-B2e-2iω0τ*)E^1=2(Δ11Δ21Δ31),
that is,
(116)(2iω0-a11-a120-a212iω0-a22b23e-2iω0τ*-a31-a322iω0-a33)(E1(1)E1(2)E1(3))=2(Δ11Δ21Δ31).
On simplification, we get
(117)E1(1)=B11C22-B22C11A11B22-A22B11,E1(2)=(2iω0-a11)E1(1)-2Δ11a12,E1(3)=C11A22-A11C22A11B22-A22B11,
where
(118)A11=-a21+(2iω0-a22)(2iω0-a11)a12,B11=-b23e-2iω0τ*,C11=-2Δ11a12(2iω0-a22),A22=-a31-a32a12(2iω0-a11),B22=2iω0-a33,C22=2Δ11a32a12-2Δ31.Similarly substituting (107) and (112) into (110), we get
(119)(-∫-τ*0‍dη(θ))E^2=2(Δ12Δ22Δ32)
which leads to
(120)(-B1-B2)E^2=2(Δ12Δ22Δ32)
that is,
(121)(-a11-a120-a21-a22-b23a31-a32-a33)(E2(1)E2(2)E2(3))=2(Δ12Δ22Δ32).
On simplification, we get
(122)E2(1)=-B11′C22′-B22′C11′A11′B22′-A22′B11′,E2(2)=a11E2(1)-2Δ12a12,E2(3)=C11′A22′-C22′A11′A11′B22′-A22′B11′,
where
(123)A11′=-a21+a11a22a12,B11′=-b23,C11′=2(Δ12a22a12-Δ22),A22′=-a31+a11a32a12,B22′=-a33,C22′=2(Δ12a32a12-Δ32).Thus,W20(θ) and W11(θ) can be determined from (90) and (93). Furthermore, g21 can also be determined. Therefore, based on the previous analysis, gij can be determined in terms of the parameters and delays involved in the system. Thus, we can compute the following values:
(124)C1(0)=i2ω0(g11g20-2|g11|2-|g02|23)+g212,μ2=-Re{C1(0)}Re{λ′(τ0)},β2=2Re{C1(0)},T2=-Im{C1(0)}+μ2Im{λ′(τ0)}ω0,
where sign of μ2 determines the direction of the Hopf bifurcation, sign of β2 determines the stability of the bifurcating periodic solution, and sign of T2 determines the period of the bifurcating periodic solution in the center manifold at the critical value τ*. Therefore, by the results of Hassard et al. [33], we can summarize the properties of the Hopf bifurcation at the critical value τ=τ* in the following theorem.Theorem 8.
In system (31) and (32), the following results hold.(i)
Ifμ2>0(<0), the Hopf bifurcation is supercritical (subcritical).(ii)
Ifβ2<0(>0), the bifurcated periodic solutions are stable (unstable).(iii)
IfT2>0(<0), period of the bifurcating periodic solution increases (decreases).


## 7. Model with Double Delays

In this section we have discussed the dynamical behaviours of the model (3) and (4) with double discrete time delays τ1 and τ2. Here, τ1 represents the activation period or reaction time of the toxicant in the population biomass; τ2 represents the activation period or reaction time of the toxicant in the population from the environment and γ(t)=Q(>0) is a constant. The modified model is as follows:(125a)dXdt=X(p-Y(t-τ1)-X)dYdt=Z(t-τ2)-(q-X)YdZdt=-aXZ+a(d+Y)XY-hZ+Q
with initial conditions
(125b)X(0)>0,Y(ϑ)=φ1(ϑ),Z(ϑ)=φ2(ϑ),suchthatφ1(ϑ)≥0,∀ϑ∈[-τ1,0],φ2(ϑ)≥0,∀ϑ∈[-τ2,0],where φ1(ϑ) is a continuous function on ϑ∈[-τ1,0] and φ2(ϑ) is a continuous function on ϑ∈[-τ2,0]. For a biological meaning, we further assume that φi(0)>0(i=1,2).Now, we will estimate the length of the delay which preserves the stability of the system (125a) and (125b). The corresponding linearized system of (125a) and (125b) about the interior equilibrium E*(X*,Y*,Z*) is given by
(126)dx1dt=a11x1+a12y1(t-τ1)dy1dt=a21x1+a22y1+b23z1(t-τ2)dz1dt=a31x1+a32y1+a33z1,
where
(127)X(t)=x1(t)+X*,Y(t)=y1(t)+Y*,Z(t)=z1(t)+Z*.
Taking Laplace transform of the system (126), we get
(128)(s-a11)x1¯(s)=a12e-sτ1y1¯(s)+a12e-sτ1k1(s)+x1(0)(s-a22)y1¯(s)=a21x1¯(s)+b23e-sτ2z1¯(s)+b23e-sτ2k1(s)+y1(0)(s-a33)z1¯(s)=a31x1¯(s)+a32y1¯(s)+z1(0),
where
(129)k1(s)=∫-τ10‍e-sty1(t)dt,k2(s)=∫-τ20‍e-stz1(t)dt,
and x1¯(s),y1¯(s),z1¯(s) are the Laplace transforms of x1(t), y1(t), and z1(t), respectively.Now, using “Nyquist theorem” [21, 31], it can be shown that the conditions for local asymptotic stability of E*(X*,Y*,andZ*) are given by [32]
(130)ImH(iη0)>0,(131)ReH(iη0)=0,
where
(132)H(s)=s3+D1s2+D2s+D3+D4se-sτ1+D5se-sτ2H(s)=+D6e-sτ1+D7e-sτ2+D8e-s(τ1+τ2),D1=aX*+h+q,D2=aqX*+hq+qX*-(X*)2,D3=aq(X*)2+qhX*-a(X*)3-h(X*)2,D4=X*Y*,D5=-a(d+Y*)X*-aX*Y*,D6=X*Y*(aX*+h),D7=-X*[a(d+Y*)X*-aX*Y*],D8=X*[-aZ*+a(d+Y*)Y*],
and η0 is the smallest positive root of (131).Nyquist criterion [21, 31] implies that if η0 is a solution of (131), then
(133)-D1η02+D3+D4η0sin(η0τ1)+D5η0sin(η0τ2)+D6cos(η0τ1)+D7cos(η0τ2)+D8cos{η0(τ1+τ2)}=0.
Next, we wish to find an upper bound for η (say, η+) independent of τ1, τ2 such that (130) and (131) hold for all η, 0≤η≤η+ and hence in particular for η0 when η0≤η+.Maximizing (133) with |sin(η0τ1)|≤1, |sin(η0τ2)|≤1, |cos(η0τ1)|≤1, |cos(η0τ2)|≤1, |cos{η0(τ1+τ2)}|≤1, we get
(134)D1η02-(|D4|+|D5|)η0-(|D3|+|D6|+|D7|+|D8|)=0.
If η+ is a positive root of (134), then we obtain
(135)η+=12D1[(|D5|)1/2(|D4|+|D5|)η+=12D1+{(|D4|+|D5|)2η+=12D1+4D1(|D3|+|D6|+|D7|+|D8|)(|D5|)2}1/2].
Then, clearly, η+≥η0.Equation (130) indicates that the following inequality should hold:
(136)φ(τ1,τ2,η0)>ψ(τ1,τ2,η0),
where
(137)φ(τ1,τ2,η0)=-η0[η02-D2-D4cos(η0τ1)-D5cos(η0τ2)],ψ(τ1,τ2,η0)=D6sin(η0τ1)+D7sin(η0τ2)+D8sin{η0(τ1+τ2)}.
If we can find φ~(τ1,τ2)>ψ~(τ1,τ2) such that
(138)φ(τ1,τ2,η0)(τ1+τ2)η0≥φ~(τ1,τ2)>ψ~(τ1,τ2)≥ψ(τ1,τ2,η0)(τ1+τ2)η0,
where 0<η0<η+, then the Nyquist criterion holds. By (138), we can estimate the values of τ1 and τ2(139)ψ(τ1,τ2,η0)(τ1+τ2)η0=D6sin(η0τ1)+D7sin(η0τ2)+D8sin{η0(τ1+τ2)}(τ1+τ2)η0≤|D6|+|D7|+|D8|.
We choose
(140)ψ~(τ1,τ2)=|D6|+|D7|+|D8|.
Now,
(141)φ(τ1,τ2,η0)(τ1+τ2)η0=η0{D2+D4cos(η0τ1)+D5cos(η0τ2)-η02}(τ1+τ2)η0φ(τ1,τ2,η0)(τ1+τ2)η0≥{D2+D4cos(η+τ1)+|D5|-η+2}(τ1+τ2).
We can choose
(142)φ~(τ1,τ2)={D2+D4cos(η+τ1)+|D5|-η+2}(τ1+τ2).
From (138), we get
(143)φ~(τ1,τ2)={D2+D4cos(η+τ1)+|D5|-η+2}(τ1+τ2)φ~(τ1,τ2)>|D6|+|D7|+|D8|=ψ~(τ1,τ2).
Thus
(144)τ1+τ2<D2+D4cos(η+τ1)+|D5|-η+2|D6|+|D7|+|D8|τ1+τ2≤|D2|+|D4|+|D5|-η+2|D6|+|D7|+|D8|.
So, we come to the following theorem.Theorem 9.
The delayed model (125a) and (125b) will be locally asymptotically stable at E*(X*,Y*,Z*) with the conditions (24) if
(145)|D2|+|D4|+|D5|>η+2,0<τ1+τ2<|D2|+|D4|+|D5|-η+2|D6|+|D7|+|D8|.Next, we will study the existence of Hopf bifurcation of system (125a) and (125b) by choosing one of the delays as a bifurcation parameter. Let us take τ2 as a bifurcation parameter. The characteristic equation corresponding to system (125a) and (125b) is
(146)λ3+D1λ2+D2λ+D3+D4λe-λτ1+D5λe-λτ2+D6e-λτ1+D7e-λτ2+D8e-λ(τ1+τ2)=0.
Hence, λ=±iη is a simple root of (146) at τ2=τ2c. By separating real and imaginary parts, we obtain
(147)-D1η2+D3+D4ηsin(ητ1)+D5ηsin(ητ2)+D6cos(ητ1)+D7cos(ητ2)+D8cos{η(τ1+τ2)}=0,-η3+D2η+D4ηcos(ητ1)+D5ηcos(ητ2)-D6sin(ητ1)-D7sin(ητ2)+D8sin{η(τ1+τ2)}=0.Differentiating equation (146) implicitly with respect to τ2, we obtain
(148)(dλdτ2)-1=D5λ(D5λ+D7+D8e-λτ1)(dλdτ2)-1=+D8τ1e-λτ1λ(D5λ+D7+D8e-λτ1)(dλdτ2)-1=+3λ2+2D1λ+D2λe-λτ2(D5λ+D7+D8e-λτ1)(dλdτ2)-1=+e-λτ1(D4-D4λτ1-D6τ1)λe-λτ2(D5λ+D7+D8e-λτ1)-τ2λ.Now,(149)Re(dλdτ2)λ=iη0-1=1M2+N2(MP+NR),
where
(150)M=-D5η02+D8η0sin(η0τ1),N=D7η0+D8η0cos(η0τ1),P=D5-D8τ1cos(η0τ1)+(-3η02+D2)cos(η0τ2)P=-2D1η0sin(η0τ2)P=+{(D4-D6τ1)cos(η0τ1)-η0τ1D4sin(η0τ1)}P=×cos(η0τ2)P=+{(D4-D6τ1)sin(η0τ1)+η0τ1D4cos(η0τ1)}P=×sin(η0τ2)R=D8τ1sin(η0τ1)+(-3η02+D2)sin(η0τ2)R=+2D1η0cos(η0τ2)R=+{(D4-D6τ1)cos(η0τ1)-η0τ1D4sin(η0τ1)}R=×sin(η0τ2)R=-{(D4-D6τ1)sin(η0τ1)+η0τ1D4sin(η0τ1)}R=×cos(η0τ2).
Cleary, the sign of M2+N2 in (149) is positive. Because of the coherence of the sign of MP+NR and Re(dλ/dτ2)λ=iη0-1, we get that if MP+NR>0, then Re(dλ/dτ2)λ=iη0-1>0.Therefore, transversality condition holds and Hopf bifurcation occurs atτ2=τ2c, provided that (MP+NR)>0.Summarizing the above analysis, we get the following result.Theorem 10.
If there existsτ2=τ2c such that (147) holds and Re(dλ/dτ2)λ=iη0-1>0, that is, (MP+NR)>0, then a Hopf bifurcation occurs at E*(X*,Y*,Z*) as τ2 passes through τ2c.


## 8. Numerical Simulations

Analytical studies can never be completed without numerical verification of the derived results. In this section, we present computer simulation of some important analytical results discussed in the previous sections. Besides verification of our analytical findings, these numerical solutions are very important from a practical point of view.We take the parameters of the system as given in Table1. Using those values under zero exogenous input, that is, γ=0, it is observed that the equilibrium E0(0,0,0) becomes unstable and interior equilibrium E*(X*,Y*,Z*) is not feasible. Only the axial equilibrium E1(p,0,0)=(2.28,0,0) is locally asymptotically stable, which is shown in Figure 1. Therefore, Theorem 2 is verified. When γ=1.25 (i.e., nonzero exogenous input) using the parameter values given in Table 2, it is observed that the equilibrium E2=(0,4.28082,12.5) exists and is locally asymptotically stable with hpq<γ, but the interior equilibrium E* does not exist, which is shown in Figure 2. On the other hand, when γ=1.25 using the parameter values given in Table 1, it is found that hpq>γ, a(q-p-d)>h, and E2=(0,2.14041,6.25) becomes unstable and E*(X*,Y*,Z*)=(0.43227,1.84773,4.59665) exists, which is shown in Figure 3. Therefore, Theorem 3 is verified.Table 1
Parameter
Values
p
2.28
q
2.92
a
4.14
d
0.54
h
0.2Table 2
Parameter
Values
p
2.28
q
2.92
a
4.14
d
0.54
h
0
.
1Figure 1
Time series plot ofX, Y, and Z for zero exogenous input (γ=0) with X(0)=1, Y(0)=2, and Z(0)=0.2; parameter values are given in Table 1.Figure 2
Time series plot ofX, Y, and Z for nonzero exogenous input (γ=1.25) and hpq<γ with X(0)=3, Y(0)=5, and Z(0)=8; parameter values are given in Table 2.Figure 3
Time series plot ofX, Y, and Z for nonzero exogenous input (γ=1.25) and hpq>γ with X(0)=1, Y(0)=1, and Z(0)=5; parameter values are given in Table 1.Using the parameter values given in Table1 with γ=1.25, it is observed that the conditions of Theorem 4 are satisfied as A1=4.9096>0, A3=0.121507>0, A1A2-A3=0.388821>0. Therefore, by Theorem 4, E* is locally asymptotically stable. Figure 3 shows that X,Y and Z approach to their steady-state values X*, Y*, and Z*, respectively, in finite time. The XY-plane, YZ-plane, XZ-plane, and XYZ-plane projections of the solution are shown in Figures 4, 5, 6, and 7, respectively. Here the conditions of Theorem 5 are also satisfied and consequently E* is globally asymptotically stable. The phase diagram is shown in Figure 8.Figure 4
X
Y-plane projection of the solution for nonzero exogenous input (γ=1.25) with X(0)=1, Y(0)=1, and Z(0)=5; parameter values are given in Table 1.Figure 5
Y
Z-plane projection of the solution for nonzero exogenous input (γ=1.25) with X(0)=1, Y(0)=1, and Z(0)=5; parameter values are given in Table 1.Figure 6
X
Z-plane projection of the solution for nonzero exogenous input (γ=1.25) with X(0)=1, Y(0)=1, and Z(0)=5; parameter values are given in Table 1.Figure 7
X
Y
Z-plane projection of the solution for nonzero exogenous input (γ=1.25) with X(0)=1, Y(0)=1, and Z(0)=5; parameter values are given in Table 1.Figure 8
Global stability ofE*(X*,Y*,Z*) for nonzero exogenous input (γ=1.25) with (X(0),Y(0),Z(0))=(2,1,3),(4,3,5),(6,4,8),(5,3,8),(7,8,4), parameter values are given in Table 1.


## 9. Discussions and Conclusions

In this paper, we have developed a single-species population model under the influence of toxicant in the population and in the environment by means of ordinary differential equations in terms of their concentrations. It is shown that the model system (3) with (4) is bounded and the solutions are positive, which in turn implies that the system is ecologically well behaved. The analysis of the existence of various equilibrium points under zero exogenous toxicant input and nonzero exogenous toxicant input leads us to Theorems 2 and 3, respectively. Theorem 2 states that under zero exogenous input the trivial equilibrium E0 is unstable, the interior equilibrium E* is not feasible, and only the axial equilibrium E1 is locally asymptotically stable. On the other hand, when exogenous input is nonzero, the existence of the interior equilibrium E* depends on some conditions and the equilibrium E2 is locally asymptotically stable under some condition. The stability criteria for E* given in Theorems 4 and 5 are the conditions for locally and globally stable coexistence of the population biomass, population toxicant, and environmental toxicant.It is proved by several researchers that the effect of time delay must be taken into account in the population models to make them ecologically more meaningful and useful. It seems also reasonable that the effect of the environmental toxicant on the organismal activities and the population growth will not be instantaneous. It must be mediated with some discrete time lag required for incubation. From this point of view, we have formulated single and double time delayed model given in (31) and (125a) and (125b), respectively, where the delays may be looked upon as the reaction time of the environment toxicant on the population biomass and the population toxicant.The analysis of the delayed models (31) and (125a) and (125b) shows more complicated behaviour than their nondelayed counterpart. The rigorous analysis of the single delayed model (31) leads us to Theorem 6 which mentions that the stability criteria in absence of delay are no longer enough to guarantee the stability in the presence of delay, but rather there is a critical value τ* of the gestation delay τ such that the system is stable for τ<τ* and becomes unstable for τ>τ* at the interior equilibrium E*(X*,Y*,Z*). From the next part of our analysis, we obtain the estimated length of delay to preserve stability at E* and the direction of Hopf bifurcation by using normal form theory and center manifold theorem. The stability of the double delayed model (125a) and (125b) is investigated by the Nyquist criteria which leads us to Theorem 9. By choosing one of the delays as a bifurcation parameter, the model (125a) and (125b) is found to undergo a Hopf bifurcation under some conditions stated in Theorem 10.Analytical studies can never be completed without numerical verification of the results. So, some numerical simulations for justifying the important theoretical results are also illustrated by using MATLAB, which shows the mathematical and ecological reliability of the proposed model.



---

*Source: 574213-2013-06-16.xml*