# A Predator-Prey Model in the Chemostat with Time Delay

Authors: Guihong Fan; Gail S. K. Wolkowicz
Journal: International Journal of Differential Equations
 (2010)

---

## Abstract

The aim of this paper is to study the dynamics of predator-prey interaction in a chemostat to determine whether including a discrete delay to model the time between the capture of the prey and its conversion to viable biomass can introduce oscillatory dynamics even though there is a globally asymptotically stable equilibrium when the delay is ignored. Hence, Holling type I response functions are chosen so that no oscillatory behavior is possible when there is no delay. It is proven that unlike the analogous model for competition, as the parameter modeling the delay is increased, Hopf bifurcations can occur.


---

## Body

## 1. Introduction

The chemostat, also known as a continuous stir tank reactor (CSTR) in the engineering literature, is a basic piece of laboratory apparatus used for the continuous culture of microorganisms. It has potential applications for such processes as wastewater decomposition and water purification. Some ecologists consider it a lake in a laboratory. It can be thought of as three vessels, the feed bottle that contains fresh medium with all the necessary nutrients, the growth chamber where the microorganisms interact, and the collection vessel. The fresh medium from the feed bottle is continuously added to the growth chamber. The growth chamber is well stirred and its contents are then removed to the collection vessel at a rate that maintains constant volume. For a detailed description of the importance of the chemostat and its application in biology and ecology, one can refer to [1, 2].The following system describes a food chain in the chemostat where a predator population feeds on a prey population of microorganisms that in turn consumes a nonreproducing nutrient that is assumed to be growth limiting at low concentrations(1.1)ṡ(t)=(s0-s(t))D0-x(t)f(s(t))η,ẋ(t)=x(t)(-D+f(s(t)))-y(t)g(x(t))ξ,ẏ(t)=-Δy(t)+y(t)g(x(t)).
Here s(t) represents the concentration of the growth limiting nutrient, x(t) the density of the prey population, and y(t) the density of the predator population. Parameter s0 denotes the concentration of the growth limiting nutrient in the feed vessel, D0 the dilution rate, η(ξ) the growth yield constant, D(Δ) the sum of the dilution rate D0 and the natural species specific death rate of the prey (predator) population, respectively. Here f(s) denotes the functional response of the prey population on the nutrient and g(x) denotes the functional response of the predator on the prey.Butler et al. [3] considered the coexistence of two competing predators feeding on a single prey population growing in the chemostat. As a subsystem of their model, they studied the global stability of system (1.1) with both f(s) and g(x) taking the form of Holling type II. They proved that under certain conditions the interior equilibrium is globally asymptotically stable with respect to the interior of the positive cone. However, they also proved that for certain ranges of the parameters there is at least one nontrivial limit cycle and conjectured that the limit cycle is unique and would be a global attractor with respect to the noncritical orbits in the open positive octant. This conjecture was partially solved by Kuang [4]. He showed that there is a range of parameters for which a unique periodic orbit exists and roughly located the position of the limit cycle.Bulter and Wolkowicz [5] studied predator mediated coexistence in the chemostat assuming D0=D=Δ. Model (1.1) was studied as a submodel. For general monotone response functions, Bulter and Wolkowicz showed that (1.1) is uniformly persistent if the sum of the break even concentrations of substrate and prey is less than the input rate of the nutrient s0. However they showed that it is necessary to specify the form of the response functions in order to discuss the global dynamics of the model. If f(s) is modelled by Holling type I or II and g(x) by Holling type I, Bulter and Wolkowicz proved that (1.1) could have up to three equilibrium points and that there is a transfer of global stability from one equilibrium point to another as different parameters are varied making conditions favorable enough for a new population to survive. In this case, there are no periodic solutions. However, even if f(s) is given by Holling type I, if g(x) is given by Holling type II, they showed that a Hopf bifurcation can occur in (1.1), and numerical simulations indicated that the bifurcating periodic solution was asymptotically stable.We include a time delay in (1.1) to model the time between the capture of the prey and its conversion to viable biomass. Our aim is to show that such a delay can induce nontrivial periodic solutions in a model where there is always a globally asymptotically stable equilibrium when delay is ignored, and hence no such periodic solutions are possible otherwise. For this reason we select the response functions of the simplest form; that is, we choose the Holling type I form for both f(s) and g(x), so that (1.2) always has a globally asymptotically stable equilibrium when the conversion process is assumed to occur instantaneously. It is interesting to note that in the analogous model of competition between two species in the chemostat, delay cannot induce oscillatory behavior for any reasonable monotone response functions (see Wolkowicz and Xia [6]).With delay modelling the time required for the predator to process the prey after it has been captured, the model is given by(1.2)ṡ(t)=(s0-s(t))D0-x(t)f(s(t))η,ẋ(t)=x(t)(-D+f(s(t)))-y(t)g(x(t))ξ,ẏ(t)=-Δy(t)+e-Δτy(t-τ)g(x(t-τ)).
For t∈[-τ,0],(1.3)s(0)=s0∈intℝ+,(x(t),y(t))=(ϕ,ψ)∈ℂ([-τ,0],intℝ+2).
Here variables s(t), x(t), y(t), and parameters s0, D0, η, ξ, D0, D, and Δ have the same interpretation as for model (1.1). Note therefore that D⩾D0, and Δ⩾D0. The additional parameter τ is a nonnegative constant modelling the time required for the conversion process. Hence, e-Δτy(t-τ) represents the concentration of the predator population in the growth chamber at time t that were available at time t-τ to capture prey and were able to avoid death and washout during the τ units of time required to process the captured prey.We analyze the stability of each equilibrium and prove that the coexistence equilibrium can undergo Hopf bifurcations. Numerical simulations appear to show that (1.2) can have a stable periodic solution bifurcating from the coexistence equilibrium as the delay parameter increases from zero. This periodic orbit can then disappear through a secondary Hopf bifurcation as the delay parameter increases further.


## 2. Scaling of the Model and Existence of Solutions

Suppose that functionsf(s) and g(s) are of Holling type I form, that is, f(s)=αs (α>0) and g(x)=kx (k>0). System (1.2) reduces to(2.1)ṡ(t)=(s0-s(t))D0-αx(t)s(t)η,ẋ(t)=x(t)(-D+αs(t))-kx(t)y(t)ξ,ẏ(t)=-Δy(t)+ke-Δτy(t-τ)x(t-τ).t>0,
Introducing the following change of variables gives:(2.2)t̆=D0t,s̆(t̆)=s(t)s0,x̆(t̆)=x(t)s0η,y̆(t̆)=y(t)ξs0η,τ̆=D0τ,D̆=DD0,Δ̆=ΔD0,k̆=ks0ηD0,ᾰ=αs0D0,ds̆(t̆)dt̆=1s0ds(t)dtdtdt̆=1s0D0ds(t)dt=1s0D0((s0-s(t))D0-αx(t)s(t)η)=1-s(t)s0-αs0D0x(t)s0ηs(t)s0=1-s̆(t̆)-ᾰx̆(t̆)s̆(t̆),dx̆(t̆)dt̆=1s0ηdx(t)dtdtdt̆=1s0ηD0dx(t)dt=1s0ηD0(x(t)(-D+αs(t))-kx(t)y(t)ξ)=x(t)s0η(-DD0+αs0D0s(t)s0)-ks0ηD0x(t)s0ηy(t)s0ηξ=x̆(t̆)(-D̆+ᾰs̆(t̆))-k̆x̆(t̆)y̆(t̆),dy̆(t̆)dt̆=1s0ηξdy(t)dtdtdt̆=1s0ηξD0dy(t)dt=1s0ηξD0(-Δy(t)+ke-Δτy(t-τ)x(t-τ))=-Δy(t)s0ηξD0+ke-Δτs0ηξD0y(t-τ)x(t-τ)=-ΔD0y(t)s0ηξ+ks0ηD0e-(Δ/D0)D0τy(t-τ)s0ηξx(t-τ)s0η=-Δ̆y̆(t̆)+k̆e-Δ̆τ̆y̆(t̆-τ̆)x̆(t̆-τ̆).
With this change of variables, omitting the  ̆’s  for convenience, system (2.1) becomes(2.3)ṡ(t)=1-s(t)-αx(t)s(t),ẋ(t)=x(t)(-D+αs(t))-ky(t)x(t),ẏ(t)=-Δy(t)+ke-Δτy(t-τ)x(t-τ),
where Δ⩾1 and D⩾1, with initial data given by (1.3). For biological significance, a point is assumed to be an equilibrium point of (2.3) only if all of its components are nonnegative.Letτ=0. Model (2.3) reduces to a special case of the model considered in [7]. If D>α, the model has only one equilibrium point (1,0,0) and it is globally asymptotically stable. If D<α and 1-D/α-ΔD/k<0, the model has a second equilibrium point (D/α,(α-D)/αD,0) and it is globally asymptotically stable. When 1-D/α-ΔD/k>0, the model has a third equilibrium point (k/(k+αΔ),Δ/k,α/(k+αD)-D/k) and it is the global attractor. Therefore, model (2.3) has no periodic solutions when the time delay is ignored. If g(x) is of Holling type II form, Butler and Wolkowicz [5] proved that a Hopf bifurcation is possible resulting in a periodic solution for a certain range of parameter values. We emphasize again here, that it is for this reason that in this paper we restrict our attention to the simplest case for both response functions, that is, Holling type I, in order to see whether delay can be responsible for periodic solutions in (1.2).Theorem 2.1.
Assuming(s0,ϕ(θ),ψ(θ))∈intℝ+×ℂ([-τ,0],intℝ+2), then there exists a unique solution (s(t),x(t),y(t)) of (2.3) passing through (s0,ϕ(θ),ψ(θ)) with s(t)>0, x(t)>0 and y(t)>0 for t∈[0,∞). The solution is bounded. In particular, given any ϵ0>0, x(t)<1+ϵ0 for all sufficiently large t.Proof.
Fort∈[0,τ], one has t-τ∈[-τ,0],x(t-τ)=ϕ(t-τ), and y(t-τ)=ψ(t-τ). System (2.3) becomes
(2.4)ṡ(t)=1-s(t)-αx(t)s(t),ẋ(t)=x(t)(-D+αs(t))-ky(t)x(t),ẏ(t)=-Δy(t)+ke-Δτϕ(t-τ)ψ(t-τ),
a system of nonautonomous ordinary differential equations with initial conditions s(0)=s0, x(0)=ϕ(0), and y(0)=ψ(0). Since the right-hand side of (2.4) is differentiable in both x and y, by Theorems 2.3, 3.1, and Corollary 4.3 in Miller and Michel [8], there exists a unique solution defined on [0,τ] satisfying (2.4). By using the method of steps in Bellman and Cooke [9], it can be shown that the solution through (s0,ϕ(θ),ψ(θ)) is defined for all t⩾0.
Now we proves(t)>0 for all t>0. From the first equation of (2.3),
(2.5)ṡ(t)=1-s(t)-αx(t)s(t).
Proceed using the method of contradiction. Suppose that there exists a first t⋆ such that s(t⋆)=0 and s(t)>0 for t∈[0,t*). Then ṡ(t⋆)⩽0. But from the first equation of (2.3)
(2.6)ṡ(t⋆)=1-s(t⋆)-αx(t⋆)s(t⋆)=1>0,
a contradiction.
To provex(t)>0 for t∈[0,∞), assume there is a first t¯>0 such that x(t¯)=0, and x(t)>0 for t∈[0,t¯). Divide both sides of the second equation of (2.3) by x(t) and integrate from 0 to t¯, to obtain
(2.7)x(t¯)=ϕ(0)exp(∫0t¯(-D+αs(t)-ky(t))dt)>0,
contradicting x(t¯)=0.
To show thaty(t) is positive on [0,∞), suppose that there exists t⋆>0 such that y(t⋆)=0, and y(t)>0 for t∈[0,t*). Then ẏ(t⋆)⩽0. From the third equation of (2.3), we have
(2.8)ẏ(t⋆)=-Δy(t⋆)+ke-Δτy(t⋆-τ)x(t⋆-τ)=ke-Δτy(t⋆-τ)x(t⋆-τ)>0,
a contradiction.
To prove the boundedness of solutions, define
(2.9)ω(t)=s(t)+x(t)+eΔτy(t+τ)-1,fort⩾0.
It follows that
(2.10)ω̇(t)=1-s(t)-Dx(t)-ΔeΔτy(t+τ)⩽1-s(t)-x(t)-eΔτy(t+τ)⩽-ω(t),
where the first inequality holds since D⩾1, Δ⩾1, x(t)>0 and y(t+τ)>0. It follows that
(2.11)s(t)+x(t)+eΔτy(t+τ)⩽1+(s0+x(0)+eΔτy(τ)-1)e-t→1ast→∞.
Therefore, the solution (s(t),x(t),y(t)) is bounded, and given any ϵ0>0, x(t)<1+ϵ0 for all sufficiently large t.


## 3. Equilibria and Stability

Model (2.3) has three equilibrium points: E1=(1,0,0),E2=(D/α,(α-D)/αD,0), and(3.1)E+=(s+(τ),x+(τ),y+(τ))=(11+(αΔ/k)eΔτ,ΔkeΔτ,αk+αΔeΔτ-Dk).
We call E1 the washout equilibrium, E2 the single species equilibrium, and E+ the coexistence equilibrium. For the sake of biological significance, E+ exists (distinct from E2) if and only if its third coordinate y+(τ)=(αs+(τ)-D)/k>0, that is, s+(τ)>D/α, or equivalently, τ lies between 0 and τc, where(3.2)τc=1Δln(kΔ(1D-1α)).
Note that if (k/Δ)(1/D-1/α)⩽1, the equilibrium E+ does not exist for any τ (⩾0), and if (k/Δ)(1/D-1/α)=1, then E+=E2.The linearization of (2.3) about an equilibrium (s,x,y) is given by(3.3)[ż1(t)ż2(t)ż3(t)]=[-1-αx-αs0αx-D+αs-ky-kx00-Δ][z1(t)z2(t)z3(t)]+[0000000keΔτykeΔτx][z1(t-τ)z2(t-τ)z3(t-τ)].
The associated characteristic equation is given by(3.4)det[-1-αx-λ-αs0αx-D+αs-ky-λ-kx0ke-Δτ-λτy-Δ+ke-Δτ-λτx-λ]=0.
Direct calculation of the left-hand side of (3.4) gives(3.5)(-Δ+ke-(Δ+λ)τx-λ){(-1-αx-λ)(-D+αs-ky-λ)+α2sx}+kxke-(Δ+λ)τy(-1-αx-λ)=(-Δ-λ){(1+αx+λ)(D-αs+ky+λ)+α2sx}+e-(Δ+λ)τkx×{ky(-1-αx-λ)+(1+αx+λ)(D-αs+ky+λ)+α2sx}=(-Δ-λ){(1+αx+λ)(D-αs+ky+λ)+α2sx}+e-(Δ+λ)τkx{(1+αx+λ)(D-αs+λ)+α2sx}=(-Δ-λ){(λ+1)(λ+D+ky)+αx(λ+D+ky)-αs(λ+1)}+e-(Δ+λ)τkx{(λ+1)(λ+D)+αx(λ+D)-αs(λ+1)}.
For convenience, define
(3.6)P(λ):=(-Δ-λ){(λ+1)(λ+D+ky)+αx(λ+D+ky)-αs(λ+1)}+e-(Δ+λ)τkx{(λ+1)(λ+D)+αx(λ+D)-αs(λ+1)}.Theorem 3.1.
EquilibriumE1 is stable if α<D and unstable if α>D.Proof.
Evaluating the characteristic equation atE1 gives
(3.7)P(λ)|E1=-(Δ+λ)(λ+1)(λ+D-α)=0.
The eigenvalues -1 and -Δ are both negative. The third eigenvalue is -D+α. Therefore the equilibrium E1 is stable if α<D and unstable if α>D.Remark 3.2.
Ifα<D, then there is only one equilibrium, E1. If α>D, equilibrium E2 also exists.Lemma 3.3.
Assumeα>D. The characteristic equation evaluated at E2 has two negative eigenvalues, and the remaining eigenvalues are solutions of
(3.8)(λ+Δ)e(λ+Δ)τ=k(1D-1α).
In addition, the characteristic equation evaluated at E2 has zero as an eigenvalue if and only if τ=τc.Proof.
Assumeα>D. Equilibrium E2 exists. Consider the characteristic equation at E2. Since (α-D)/αD=(1-s)/αs at E2,
(3.9)P(λ)|E2={(λ+1)(λ+D)+αx(λ+D)-αs(λ+1)}×(-λ-Δ+e-(Δ+λ)τkx)={(λ+1)(λ+D)+1-ss(λ+D)-D(λ+1)}×(-λ-Δ+e-(Δ+λ)τkα-DαD)={λ(λ+1)-(λ+D)+λ+Ds}(-λ-Δ+kα-DαDe-(Δ+λ)τ)=-(λ2+αDλ+α-D)(λ+Δ-kα-DαDe-(Δ+λ)τ)=-e-(Δ+λ)τ(λ-λ1)(λ-λ2)((λ+Δ)e(Δ+λ)τ-k(1D-1α))=0,
where λ1+λ2=-α/D and λ1λ2=α-D>0. Therefore, λ1 and λ2 have negative real parts. The rest of the eigenvalues are roots of (3.8).
Assuming thatλ=0 is a root of (3.8), we have
(3.10)ΔeΔτ=k(1D-1α).
Solving for τ gives
(3.11)τ=1Δln(kΔ(1D-1α))=τc.Theorem 3.4.
Assume thatD⩾1, Δ⩾1, k>0, α>0, and (k/Δ)(1/D-1/α)⩾1 so that τc⩾0. Equilibrium E2 is locally asymptotically stable if τ>τc and unstable if τ<τc. If D=1, then equilibrium E2 is globally asymptotically stable for τ>(1/Δ)ln(k/Δ).Proof.
Assume thatτ>τc. Assumptions k>0, Δ⩾1, and (k/Δ)(1/D-1/α)⩾1 imply 1/D>1/α, or equivalently α>D. By Lemma 3.3, to prove that equilibrium E2 is locally asymptotically stable, one only needs to show that (3.8) admits no root with nonnegative real part.
Consider the real roots of (3.8) first. Note that 1/D>1/α. Equation (3.8) has no solution for λ⩽-Δ. Otherwise the left-hand side would be less than zero, but the right-hand side would be greater than zero. Assume λ>-Δ. The left-hand side of (3.8) is a monotone increasing function in both λ and τ, takes value 0 at λ=-Δ, and goes to positive infinity as λ→+∞ or τ→+∞. By Lemma 3.3, when τ=τc, then λ=0 is a solution of (3.8). Thus for τ>τc, any real root λ of (3.8) must satisfy -Δ<λ<0.
For anyτ=τ̃<τc, we have (λ+Δ)e(λ+Δ)τ|τ=τ̃,λ=0<k(1/D-1/α) and limλ→+∞(λ+Δ)e(λ+Δ)τ̃=+∞. Therefore there exists at least one λ=λ̃>0 such that (τ̃,λ̃) is a solution of (3.8). Equilibrium E2 is unstable if τ<τc.
In what follows, we prove that ifτ>τc all complex eigenvalues of (3.8) have negative real parts. Suppose that λ+Δ=γ+iβ(β>0) is a solution of (3.8). Using the Euler formula, we have
(3.12)γcos(βτ)-βsin(βτ)+i(γsin(βτ)+βcos(βτ))=k(1D-1α)e-γτ.
Equating the real parts and imaginary parts of the equation, we have
(3.13)γcos(βτ)-βsin(βτ)=k(1D-1α)e-γτγsin(βτ)+βcos(βτ)=0.
Squaring both equations, adding, and taking the square root on both sides give
(3.14)γ2+β2eγτ=k(1D-1α).
The left-hand side of (3.14) is monotonically increasing in γ, β, and τ provided that γ>0. Since (3.14) has solution γ=Δ, β=0 at τ=τc, any roots of (3.14) must satisfy γ<Δ since τ>τc. Hence Re{λ}=γ-Δ<0. Therefore (3.8) has no complex eigenvalue with nonnegative real part and so E2 is locally asymptotically stable for τ>τc.
Assume thatD=1. Now we prove that E2 is globally asymptotically stable when τ>(1/Δ)ln(k/Δ), or equivalently ke-Δτ<Δ. In this case, choose ϵ0>0 small enough such that ke-Δτ(1+ϵ0)<Δ. By Theorem 2.1, for such ϵ0, there exists a T>0 so that 0<x(t)<1+ϵ0 for t>T. Hence, for t>T+τ, ke-Δτx(t-τ)<Δ. In Example 5.1 of Kuang ([10, page 32]), choose ρ(t)=τ, a(t)=Δ, b(t)=ke-Δτx(t-τ), and α=Δ/2. We obtain (ke-Δτ(1+ϵ0))2<Δ2=4(Δ-α)α. Therefore y(t)→0 as t→∞. Let z(t)=s(t)+x(t). Noting D=1, from (2.3), we have ż(t)=1-z(t)-kx(t)y(t). Multiply by the integrating factor et, (z(t)et)′=et(1-kx(t)y(t)). Integrating both sides from 0 to t gives
(3.15)z(t)=e-tz(0)+e-t(et-1)-e-t∫0teskx(s)y(s)ds=1+e-t(z(0)-1)-e-t∫0teskx(s)y(s)ds.
If limt→∞∫0teskx(s)y(s)ds<∞, then limt→∞e-t∫0teskx(s)y(s)ds=0. Therefore limt→∞z(t)=1. If limt→∞∫0teskx(s)y(s)ds=∞, by L'Hôspital's rule,
(3.16)limt→+∞e-t∫0teskx(s)y(s)ds=limt→+∞∫0teskx(s)y(s)dset=limt→+∞etkx(t)y(t)et=limt→+∞kx(t)y(t)=0,
since x(t) is bounded and limt→∞y(t)=0. It again follows that limt→∞z(t)=1. Hence
(3.17)limt→∞s(t)+x(t)=1.
We show thatlimt→∞s(t)=1/α and limt→∞x(t)=(α-1)/α. First assume that the limits exist, that is, limt→∞s(t)=s¯ and limt→∞x(t)=x¯. From (2.3), we know that  ṡ(t) and ẋ(t) are uniformly continuous since s(t), x(t), and y(t) are bounded. By Theorem A.3, it follows that limt→∞ṡ(t)=0 and limt→∞ẋ(t)=0. Note that limt→∞y(t)=0. Letting t→∞ in (2.3) gives
(3.18)(1-s¯)-αx¯s¯=0,x¯(-1+αs¯)=0.
Either (s¯,x¯)=(1,0) or (s¯,x¯)=(1/α,(α-1)/α). Assume that (s¯,x¯)=(1,0), that is, limt→∞s(t)=1 and limt→∞x(t)=0. Note that α>D. There exists ϵ>0 such that α-D-(α+k)ϵ>0. For such ϵ, there exists a sufficiently large t so that s(t)>1-ϵ and 0<y(t)<ϵ. Recalling that x(t)>0, by (2.3)
(3.19)ẋ(t)>x(t)(-D+α(1-ϵ)-kϵ)=x(t)(α-D-αϵ-kϵ)>0,
for all sufficiently large t. Therefore it is impossible for x(t) to approach 0 from above giving a contradiction. Therefore, we must have (s¯,x¯)=(1/α,(α-1)/α).
Now suppose that the limits do not exist. In particular ifx(t) does not converge, then let x¯=limsupt→∞x(t) and x̲=liminft→∞x(t). By Lemma A.2 in the appendix, there exists {tm}↑∞ and {sm}↑∞ such that
(3.20)limm→∞x(tm)=x¯,limm→∞ẋ(tm)=0,limm→∞x(sm)=x̲limm→∞ẋ(sm)=0.
From (2.3),
(3.21)x(tm)(-D+αs(tm)+ky(tm))=0.
Noting that x(tm)>0, we have s(tm)=(1-ky(tm))/α. Since limt→∞y(t)=0, limt→∞s(tm)=1/α. By (3.17), limt→∞x(tm)=limt→∞(x(tm)+s(tm))-s(tm)=1-1/α=(α-1)/α. Therefore x¯=(α-1)/α. Similarly we can show that x̲=(α-1)/α. This implies that limt→∞x(t)=(α-1)/α, a contradiction.
Sinces(t)+x(t) converges and x(t) converges, then s(t) must also converge. Hence limt→∞s(t)=1/α and limt→∞x(t)=(α-1)/α. It follows that E2 is globally asymptotically stable.


## 4. Hopf Bifurcations atE+ Assuming D=Δ=1

Now consider the stability ofE+. The characteristic equation at E+ is(4.1)P(λ)|E+=(-Δ-λ)((1+αx+(τ)+λ)(D-αs+(τ)+ky+(τ)+λ)+α2s+(τ)x+(τ))+e-(Δ+λ)τkx+(τ)((λ+1)(λ+D)+αx+(τ)(λ+D)-αs+(τ)(λ+1))=(-Δ-λ)((1+αx+(τ)+λ)λ+α2s+(τ)x+(τ))+e-(Δ+λ)τkx+(τ)((λ+1)(λ+D)+αx+(τ)(λ+D)-αs+(τ)(λ+1))=(-Δ-λ)((1s+(τ)+λ)λ+α(1-s+(τ)))+e-λτΔ((λ+1)(λ+D)+1-s+(τ)s+(τ)(λ+D)-αs+(τ)(λ+1))=(-Δ-λ)(λ2+λs+(τ)+α(1-s+(τ)))+Δe-λτ((λ+1s+(τ))(λ+D)-αs+(τ)(λ+1))=0.
By assumption Δ=D=1, and so(4.2)P(λ)|E+=-(λ+1)(λ2+λs+(τ)+α(1-s+(τ))+e-λτ(-λ+αs+(τ)-1s+(τ)))=-(λ+1)(λ2+p(τ)λ+β(τ)+e-λτ(qλ+c(τ)))=0,
where(4.3)p(τ)=1s+(τ),β(τ)=α(1-s+(τ)),q=-1,c(τ)=αs+(τ)-1s+(τ).
The characteristic equation at E+ has one eigenvalue equal to -1 and the others are given by solutions of the equation(4.4)λ2+p(τ)λ+β(τ)+e-λτ(qλ+c(τ))=0.Lemma 4.1.
Assumingk>0, α>0, and k(1-1/α)⩾1 so that τc=ln(k(1-1/α))⩾0, then E+ has no zero eigenvalue for τ∈(0,τc).Proof.
Assume thatτ∈(0,τc). By the method of contradiction, suppose that there exists a zero root of (4.4). Therefore
(4.5)β(τ)+c(τ)=α-1s+(τ)=0.
Noting that τc>0 if and only if k(1-1/α)>1, for any 0<τ<τc,
(4.6)α-1s+(τ)=α-1-αkeτ>α-1-α(1-1α)=0,
a contradiction.Lemma 4.2.
Assumek>0, α>0, k(1-1/α)>1. Equilibrium E+ is asymptotically stable when τ=0.Proof.
Forτ=0, (4.4) reduces to
(4.7)λ2+p(0)λ+β(0)+(qλ+c(0))=λ2+(1s+(0)-1)λ+α-1s+(0).
Both coefficients are positive, since
(4.8)1s+(0)-1=αk>0,α-1s+(0)=α-1-αk=α(1-1α-1k)>0,
and k(1-1/α)>1 implies 1-1/α>1/k. Therefore, all the roots of the characteristic equation have negative real parts.Lemma 4.3.
Asτ is increased from 0, a root of (4.4) with positive real part can only appear if a root with negative real part crosses the imaginary axis.Proof.
Takingn=2 and g(λ,τ)=p(τ)λ+(qλ+c(τ))e-λτ+β(τ) in Kuang [10, Theorem 1.4, page 66] gives
(4.9)lim supReλ>0,|λ|→∞|λ-2g(λ,τ)|=0<1.
Therefore, no root of (4.4) with positive real part can enter from infinity as τ increases from 0. Hence roots with positive real part can only appear by crossing the imaginary axis.
Forτ≠0, assuming λ=iω (ω>0) is a root of P(λ)|E+=0,
(4.10)-ω2+ip(τ)ω+β(τ)+e-iωτ(iqω+c(τ))=0.
Substituting eiθ=cosθ+isinθ into (4.10) gives
(4.11)-ω2+β(τ)+qωsin(ωτ)+c(τ)cos(ωτ)+i(p(τ)ω+qωcos(ωτ)-c(τ)sin(ωτ))=0.
Separating the real and imaginary parts, we obtain
(4.12)c(τ)cos(ωτ)+qωsin(ωτ)=ω2-β(τ),c(τ)sin(ωτ)-qωcos(ωτ)=p(τ)ω.
Solving for cos(ωτ) and sin(ωτ) gives
(4.13)sin(ωτ)=c(τ)p(τ)ω+qω(ω2-β(τ))c(τ)2+q2ω2,cos(ωτ)=c(τ)(ω2-β(τ))-qp(τ)ω2c(τ)2+q2ω2.
Noting sin2(ωτ)+cos2(ωτ)=1, squaring both sides of equations (4.13), adding, and rearranging gives
(4.14)ω4+(p2(τ)-q2-2β(τ))ω2+β2(τ)-c2(τ)=0.
Solving for ω, we obtain two roots ω1(τ) and ω2(τ):
(4.15)ω1(τ)=12(q2-p2(τ)+2β(τ)+(q2-p2(τ)+2β(τ))2-4(β2(τ)-c2(τ)))1/2=1s+(τ)2((1-s+(τ))(2αs+2(τ)-s+(τ)-1)+(s+2(τ)-1)2+4αs+2(τ)(s+2(τ)-1)(1-s+(τ))+4s+2(τ)(αs+2(τ)-1)2)1/2ω2(τ)=12(q2-p2(τ)+2β(τ)-(q2-p2(τ)+2β(τ))2-4(β2(τ)-c2(τ)))1/2=1s+(τ)2((1-s+(τ))(2αs+2(τ)-s+(τ)-1)-(s+2(τ)-1)2+4αs+2(τ)(s+2(τ)-1)(1-s+(τ))+4s+2(τ)(αs+2(τ)-1)2)1/2.
Define conditions (H1) and (H2) as follows:
(4.16)(H1)q2-p2(τ)+2β(τ)>0,β2(τ)-c2(τ)>0,(q2-p2(τ)+2β(τ))2-4(β2(τ)-c2(τ))⩾0,(4.17)(H2)β2(τ)-c2(τ)<0,orβ2(τ)-c2(τ)=0andq2-p2(τ)+2β(τ)>0.Lemma 4.4.
If (H1) holds for all τ in some interval I, then (4.14) has two positive roots ω1(τ)⩾ω2(τ) for all τ∈I with ω1(τ)>ω2(τ) when all the inequalities in (H1) are strict. If (H2) holds for all τ in some interval I, then (4.14) has only one positive root, ω1(τ) for all τ∈I. If no interval exists where either (H1) or (H2) holds, then there are no positive real roots of (4.14).Define the interval(4.18)J=[ln(kα(11/4+1/16+1/(2α)-1)),ln(k(α4-1)α)].
When the end points of J are real and J≠∅, define(4.19)I1=[0,τc)∩J.
We prove that (H1) holds for any τ∈I1.FromD=Δ=1,(4.20)τc=1Δln(kΔ(1D-1α))=ln(k(α-1)α).
If α>1, then α>α4. It follows that(4.21)τc>ln(k(α4-1)α).
Therefore,(4.22)I1=[max{0,ln(kα(11/4+1/16+1/(2α)-1))},ln(k(α4-1)α)].Theorem 4.5.
Assumeα>(7+35)/2 and k>α/(α4-1), then I1 is not empty, and for any τ∈I1, but τ≠ln((k/α)(1/(1/4+1/16+1/(2α))-1)), condition (H1) holds and ω1(τ)>ω2(τ)>0. If τ=ln((k/α)(1/(1/4+1/16+1/(2α))-1))∈I1, then ω1(τ)>ω2(τ)=0.Proof.
For anyα>(7+35)/2, we have 1-1/α4>0, and therefore
(4.23)1α4+12-52<1(7+35)/24+12-52=0.
Hence,
(4.24)(1α4-14)2-(116+12α)2=1α-12α4-12α=-12α4((1α4)3+1-2α4)=12α4(1-1α4)((1α4)2+1α4-1)=12α4(1-1α4)((1α4+12)2-54)=12α4(1-1α4)(1α4+12-52)(1α4+12+52)<0.
Therefore, 1/α4-1/4<1/16+1/(2α). Since 1/4+1/16+1/(2α)<1/4+1/16+1/(7+35)<1, it follows that
(4.25)1α4<14+116+12α<1.
Hence,
(4.26)ln(kα(11/4+1/16+1/(2α)-1))<ln(k(α4-1)α).
From k>α/(α4-1), we have ln(k(α4-1)/α)>0. Therefore,
(4.27)max{0,ln(kα(11/4+1/16+1/(2α)-1))}<ln(k(α4-1)α),
and so I1 is not empty. Noting s+(τ)=1/(1+(αΔ/k)eΔτ) and Δ=1, for any τ∈I1, but τ≠ln((k/α)(1/(1/4+1/16+1/(2α))-1)), we have s+(τ)∈[1/α4,1/4+1/16+1/(2α)).
In what follows, we intend to show that for any suchτ, condition (H1) holds. From (4.3),
(4.28)q2-p2(τ)+2β(τ)=(-1)2-1s+2(τ)+2α(1-s+(τ))=(1-s+(τ))2αs+2(τ)(s+(τ)2-s+(τ)2α-12α)=(1-s+(τ))2αs+2(τ)((s+(τ)-14α)2-116α2-12α).
Since s+(τ)<1, to show that the first inequality in (H1) holds, it suffices to show that the factor on the right-hand side of the above expression is positive. Since α>(7+35)/2, 1/α-1/(4α)=(1/α)(1-1/(4α))>0, and
(4.29)(116α2+12α)2-(1α-14α)2=116α2+12α-1α+12αα-116α2=12α(1α-1)<0.
Since 1/α<1/α4 for α>(7+35)/2,
(4.30)14α+116α2+12α<1α<1α4.
For any s+(τ)>1/α4,
(4.31)s+(τ)-14α⩾1α4-14α>116α2+12α.
Hence,
(4.32)(s+(τ)-14α)2⩾116α2+12α.
Next consider the second inequality in (H1). For α>(7+35)/2, since 1/α4>1/α, s+(τ)⩾1/α4>1/α. Therefore, αs+(τ)>1. For s+(τ)∈[1/α4,1/4+1/16+1/(2α))(4.33)β2(τ)-c2(τ)=(β(τ)-c(τ))(β(τ)+c(τ))=(α-2αs+(τ)+1s+(τ))(α-1s+(τ))=-2αs+2(τ)(s+2(τ)-s+(τ)2-12α)(αs+(τ)-1)=-2αs+2(τ)((s+(τ)-14)2-116-12α)(αs+(τ)-1)>0.
Finally,
(4.34)(q2-p2(τ)+2β(τ))2-4(β2(τ)-c2(τ))=(q2-p2(τ))(q2-p2(τ)+4β(τ))+4c2(τ)=(1-1s+2(τ))(1-1s+2(τ)+4α(1-s+(τ)))+4(αs+(τ)-1s+(τ))2=(1-1s+2(τ))2+4α(1-1s+2(τ))(1-s+(τ))+4s+(τ)2α2-8α+4s+2(τ)=4s+(τ)2α2+4α((1-1s+2(τ))(1-s+(τ)))-8α+(1-1s+2(τ))2+4s+2(τ)=4s+(τ)2α2+4α((1-1s+2(τ))(1-s+(τ))-2)+(1+1s+2(τ))2=(α-α1)(α-α2),
where
(4.35)α1=2-(1-1/s+2(τ))(1-s+(τ))+(1/s+2(τ)+2s+(τ)+1)(1/s+(τ)-1)22s+2(τ),α2=2-(1-1/s+2(τ))(1-s+(τ))-(1/s+2(τ)+2s+(τ)+1)(1/s+(τ)-1)22s+2(τ).
Since s+(τ)<1,
(4.36)2-(1-1s+2(τ))(1-s+(τ))=s+(τ)+1+1-s+(τ)s+2(τ)>0,(2-(1-1s+2(τ))(1-s+(τ)))2>(2-(1-1s+2(τ))(1-s+(τ)))2-s+2(τ)(1+1s+2(τ))2=(1s+2(τ)+2s+(τ)+1)(1s+(τ)-1)2>0.
It follows that 0<α2<α1. Again noting that s+(τ)<1,
(4.37)α1<2-(1-1/s+2(τ))(1-s+(τ))+(1/s+2(τ)+2/s+(τ)+1)(1/s+(τ)-1)22s+2(τ)=2-(1-1/s+2(τ))(1-s+(τ))+(1/s+(τ)+1)(1/s+(τ)-1)2s+2(τ)=2-(1-s+(τ)-1/s+2(τ)+1/s+(τ))+1/s+2(τ)-12s+2(τ)=s+(τ)+2/s+2(τ)-1/s+(τ)2s+2(τ)=12(2s+4(τ)-(1s+3(τ)-1s+(τ)))<122s+4(τ)=1s+4(τ).
Hence, for any s+(τ)>1/α4, we have α>1/s+4(τ)>α1>α2. This leads to
(4.38)(q2-p2(τ)+2β(τ))2-4(β2(τ)-c2(τ))=(α-α1)(α-α2)>0.
Therefore, (H1) holds for any τ∈I1. By Lemma 4.4, both ω1(τ)>0 and ω2(τ)>0.
Ifτ=ln((k/α)(1/(1/4+1/16+1/(2α))-1))∈I1, we have s+(τ)=1/4+1/16+1/(2α). Noting (4.25), we obtain
(4.39)q2-p2(τ)+2β(τ)>0,β2(τ)-c2(τ)=0,(q2-p2(τ)+2β(τ))2-4(β2(τ)-c2(τ))>0.
By (4.15), it follows that ω1(τ)>0 and ω2(τ)=0.Now we define intervalI2 and prove that (H2) holds on I2(4.40)I2:=[0,τc)∩(-∞,ln(kα(11/4+1/16+1/(2α)-1))].
In the following theorem, we consider the case that parameters are chosen so that(4.41)I2=[0,ln(kα(11/4+1/16+1/(2α)-1))].Theorem 4.6.
Assumeα>1 and k>α(1/(1/4+1/16+1/(2α))-1)-1. Interval I2 given by (4.41) is not empty. For any τ∈I2, (H2) holds and hence ω1(τ)>0.Proof.
Assumeα>1. Letting
(4.42)G(α)=14+116+12α-1α,
then
(4.43)ddαG(α)=-14α21/16+1/(2α)+1α2=1α2(-11+8/α+1)>0.G(α) is an increasing function of α and G(1)=0. G(α)>G(1) implies that 1/4+1/16+1/(2α)-1/α>0. Therefore
(4.44)1>14+116+12α>1α.
This gives
(4.45)α-1>11/4+1/16+1/(2α)-1>0.
By assumption k>α(1/(1/4+1/16+1/(2α))-1)-1, we obtain
(4.46)k(α-1)α>kα(11/4+1/16+1/(2α)-1)>1.
Noting that D=Δ=1 and recalling the definition of τc given in (3.2),
(4.47)τc=ln(k(α-1)α)>ln(kα(11/4+1/16+1/(2α)-1))>0.
Therefore, I2 given by (4.41) is not empty. For any τ∈I2, noting s+(τ)=1/(1+(αΔ/k)eΔτ) and Δ=1, we have s+(τ)∈[1/4+1/16+1/(2α),1/(1+α)/k)⊂[1/4+1/16+1/(2α),1).
In what follows, we intend to show for anyτ∈I2, or equivalently s+(τ)∈[1/4+1/16+1/(2α),1), (H2) holds. For any s+(τ)∈[1/4+1/16+1/(2α),1), by (4.44), it follows that s+(τ)>1/α and so αs+(τ)>1. Hence,
(4.48)β2(τ)-c2(τ)=-2αs+2(τ)((s+(τ)-14)2-116-12α)(αs+(τ)-1)⩽0.
For any s+(τ)∈[1/4+1/16+1/(2α),1), we have s+(τ)>1/(4α)+1/(16α2)+1/(2α), since 1/4>1/(4α) and 1/16>1/(16α2) imply that
(4.49)14+116+12α>14α+116α2+12α.
Therefore,
(4.50)q2-p2(τ)+2β(τ)=(1-s+(τ))2αs+2(τ)((s+(τ)-14α)2-116α2-12α)>0.
Condition(H2) holds. By Lemma 4.4, ω1(τ)>0.Next, to determine whether (4.2) has a pair of pure imaginary eigenvalues, we consider(4.51)c2(τ)+q2ω2=(αs+(τ)-1s+(τ))2+(-1)2ω2=(αs+(τ)-1s+(τ))2+ω2,c(τ)(p(τ)ω)+qω(ω2-β(τ))=ω(c(τ)p(τ)+q(ω2-β(τ)))=ω(α-1s+2(τ)-ω2+α(1-s+(τ))),c(τ)(ω2-β(τ))+qω(-p(τ)ω)=(αs+(τ)-1s+(τ))(ω2-α(1-s+(τ)))+(-1)ω(-ωs+(τ))=α(s+(τ)ω2-(1-s+(τ))(αs+(τ)-1s+(τ))).
We obtain(4.52)sin(ωτ)=-ωω2-(α-1/s+2(τ))-α(1-s+(τ))(αs+(τ)-1/s+(τ))2+ω2,cos(ωτ)=αs+(τ)ω2-(1-s+(τ))(α-1/s+2(τ))(αs+(τ)-1/s+(τ))2+ω2.If there exists(τ,ω) satisfying (4.52), then (4.2) has a pair of pure imaginary roots ±iω. A necessary condition for (4.52) to have solutions is αs+2(τ)≠1. Otherwise, αs+2(τ)=1, and the second equation of (4.52) becomes cos(ωτ)=αs+(τ). However, for any τ∈(0,τc), we have αs+(τ)>1, since(4.53)αs+(τ)=α1+αeτ/k>α1+(α/k)k(1-1/α)=α1+α(1-1/α)=αα=1.
Hence, the second equation of (4.52) has no solution. Assume αs+2(τ)≠1 for ω⩾0 and τ∈[0,τc] and denote the right-hand sides of (4.52) by(4.54)h1(ω,τ)=c(τ)(p(τ)ω)+qω(ω2-β(τ))c2(τ)+q2ω2=-ωω2-(α-1/s+2(τ))-α(1-s+(τ))(αs+(τ)-1/s+(τ))2+ω2,h2(ω,τ)=c(τ)(ω2-β(τ))+qω(-p(τ)ω)c2(τ)+q2ω2=αs+(τ)ω2-(1-s+(τ))(α-1/s+2(τ))(αs+(τ)-1/s+(τ))2+ω2.
Define functions(4.55)θ1(τ)=arccos(h2(ω1(τ),τ))ifh2(ω1(τ),τ)∈[-1,1],θ2(τ)=arccos(h2(ω2(τ),τ))ifh2(ω2(τ),τ)∈[-1,1].Lemma 4.7.
Assumeα>(7+35)/2 and k>α/(α4-1). For any τ∈I1 given by (4.19), there exists ϵj>0 and θj(τ) with ϵj⩽θj(τ)⩽π (j=1,2) such that
(4.56)sin(θj(τ)+2nπ)=h1(ωj(τ),τ),cos(θj(τ)+2nπ)=h2(ωj(τ),τ).n=0,1,2,…,Proof.
For anyτ∈I1, by Theorem 4.5, ω1(τ)>0 and ω2(τ)⩾0. It is easy to see that h1(0,τ)=0 and limω→+∞h1(ω,τ)=-∞. There are two roots of h1(z,τ)=0, z1=0 and
(4.57)z2(τ)=α-1s+2(τ)+α(1-s+(τ)).
Hence, h1(ω,τ)>0 for 0<ω<z2(τ). For any τ∈I1, as shown in Theorem 4.5, s+(τ)∈[1/α4,1/4+1/16+1/(2α)]. This implies s+(τ)⩾1/α4>1/α>1/α. Therefore, αs+2(τ)>1 and αs+(τ)>1. The function h2(ω,τ) is monotonically increasing for ω⩾0, since
(4.58)∂h2(ω,τ)∂ω=αs+(τ)2ω((αs+(τ)-1/s+(τ))2+ω2)-2ω(ω2-(1-s+(τ))(α-1/s+2(τ)))((αs+(τ)-1/s+(τ))2+ω2)2=2αs+(τ)ω(αs+(τ)-1/s+(τ))2+(1-s+(τ))(α-1/s+2(τ))((αs+(τ)-1/s+(τ))2+ω2)2=2αs+(τ)ω(α-1/s+2(τ))(s+2(τ)(α-1/s+2(τ))+1-s+(τ))((αs+(τ)-1/s+(τ))2+ω2)2=2αs+(τ)ω(α-1/s+2(τ))s+(τ)(αs+(τ)-1)((αs+(τ)-1/s+(τ))2+ω2)2=2αω(αs+2(τ)-1)(αs+(τ)-1)((αs+(τ)-1/s+(τ))2+ω2)2⩾0.
Since s+(τ)<1,
(4.59)h2(0,τ)=-αs+(τ)(1-s+(τ))(α-1/s+2(τ))(αs+(τ)-1/s+(τ))2=-αs+(τ)(1-s+(τ))αs+2(τ)-1<0.
Also, limω→∞h2(ω,τ)=αs+(τ)>1. Therefore, there exists a unique ω=lmax(τ)=α-1/s+2(τ)>0, such that h2(lmax(τ),τ)=1. Solving h2(lmax,τ)=1 for lmax and noting that α-1/s+2(τ)≠0, it is easy to see that lmax(τ)=α-1/s+2(τ)(4.60)lmax2(αs+(τ)-1)=αs+(τ)(1-s+(τ))(α-1s+2(τ))+(αs+(τ)-1s+(τ))2=αs+(τ)(α-1s+2(τ)-αs+(τ)+1s+(τ))+α2s+2(τ)-2α+1s+2(τ)=α2s+(τ)-αs+(τ)-α2s+2(τ)+α+α2s+2(τ)-2α+1s+2(τ)=α2s+(τ)-αs+(τ)-α+1s+2(τ)=(α-1s+2(τ))(αs+(τ)-1).
Then, h2(ω,τ)⩽1 for any ω∈[0,lmax(τ)]. Since s+(τ)<1, lmax(τ)<z2(τ). Therefore, h1(ω,τ)>0 for any ω∈[0,lmax(τ)]. Since ω1(τ) is a positive root of h12(ω,τ)+h22(ω,τ)=1, we have h2(ω1(τ),τ)⩽1, which implies that 0<ω1(τ)⩽lmax(τ)<z2(τ). Therefore, h1(ω1(τ),τ)>0, and so h1(ω1(τ),τ)=1-h2(ω1(τ),τ). In fact,
(4.61)ω1(τ)<lmax(τ),
since
(4.62)h22(lmax(τ),τ)+h12(lmax(τ),τ)=1+h12(lmax(τ),τ)>1.
Thus, θ1(τ) is defined and 0⩽θ1(τ)⩽π. Since cos(θ1(τ)+2nπ)=h2(ω1(τ),τ),
(4.63)sin(θ1(τ)+2nπ)=1-cos2(θ1(τ)+2nπ)=1-h22(ω1(τ),τ)=h1(ω1(τ),τ).
Hence, θ1(τ) satisfies (4.56). From (4.61), h2(ω1(τ),τ)<h2(lmax(τ),τ)=1, and so θ1(τ)>0. Since θ1(τ) is continuous on the interval I1 and I1 is closed, there exists ϵ1>0 such that θ1(τ)⩾ϵ1. Similarly we can prove the existence of θ2(τ).Lemma 4.8.
Assumeα>1 and k>α(1/(1/4+1/16+1/2α)-1)-1. For any τ∈I2 given by (4.41), there exists ϵ>0 and θ1(τ) such that ϵ⩽θ1(τ)<π and θ1(τ) satisfies (4.56) for j=1.Proof.
For anyτ∈I2, by Theorem 4.6, only ω1(τ)>0.
As in Lemma4.7, we have h1(ω,τ)>0 for 0<ω<z2(τ). For any τ∈I2, as shown in Theorem 4.6, s+(τ)∈[1/4+1/16+1/(2α),1). Letting
(4.64)G(α)=14+116+12α-1α,α>1,
we have
(4.65)ddαG(α)=-14α21/16+1/(2α)+12αα=12α(-1α2/4+2α+1α)>0.G(α) is an increasing function and G(1)=0. Since G(α)>G(1), it follows that 1/4+1/16+1/(2α)>1/α. Therefore, for any s+(τ)>1/4+1/16+1/(2α), we obtain s+(τ)>1/α, or equivalently α>1/s+2(τ). Since
(4.66)∂h2(ω,τ)∂ω=2αs+(τ)ω(αs+(τ)-1/s+(τ))2+(1-s+(τ))(α-1/s+2(τ))((αs+(τ)-1/s+(τ))2+ω2)2⩾0,h2(ω,τ) is monotonically increasing for any ω⩾0. For any s+(τ)>1/4+1/16+1/(2α),
(4.67)(s+(τ)-14)2-116-12α=s+2(τ)-s+(τ)2-12α=12α(2αs+2(τ)-αs+(τ)-1)=12α(αs+2(τ)-1-(αs+(τ)-αs+2(τ)))=αs+2(τ)-12α(1-αs+(τ)(1-s+(τ))αs+2(τ)-1)>0,
which implies that (αs+(τ)(1-s+(τ)))/(αs+2(τ)-1)<1. Hence,
(4.68)0>h2(0,τ)=-αs+(τ)(1-s+(τ))(α-1/s+2(τ))(αs+(τ)-1/s+(τ))2=-αs+(τ)(1-s+(τ))αs+2(τ)-1>-1.
For τ∈[0,τc), limω→∞h2(ω,τ)=αs+(τ)>1, since
(4.69)αs+(τ)=α1+αeτ/k>α1+(α/k)k(1-1/α)=α1+α(1-1/α)=αα=1.
As in the proof of Lemma 4.7, there exists a unique lmax(τ)=α-1/s+2(τ)>0 such that h2(lmax(τ),τ)=1. Then lmax(τ)<z2(τ). Therefore, h1(ω,τ)>0 for any ω∈[0,lmax(τ)]. The rest of the proof is similar to the proof of Lemma 4.7. Furthermore, θ1(τ)<π, since h2(ω1(τ),τ)>-1 for any ω1(τ)∈[0,lmax(τ)).Theorem 4.9.
Consider system (2.3) with D=Δ=1. (1)
Supposeα>(7+35)/2, k>α/(α4-1), and τ∈I1 given by (4.22). For τ∈I1 and j=1,2, ωj(τ) is nonnegative and there exists ϵj>0 and θj(τ) such that ϵj⩽θj(τ)⩽π and θj(τ) satisfies (4.56). If there exists n⩾0 such that θj(τ)+2nπ intersects τωj(τ) at some τnj∈I1, then (4.4) has a pair of pure imaginary eigenvalues λ=±iωj(τnj). System (2.3) undergoes a Hopf bifurcation at τnj provided dRe(λ(τ))/dτ|τ=τnj≠0.(2)
Supposeα>1, k>α(1/(1/4+1/16+1/(2α))-1)-1, and τ∈I2 given by (4.41). For τ∈I2, only ω1(τ) is positive and there exists ϵ>0 and θ1(τ) such that ϵ⩽θ1(τ)<π and θ1(τ) satisfies (4.56) for j=1. If there exists n⩾0 such that θ1(τ)+2nπ intersects τω1(τ) at some τn1∈I2, then (4.4) has a pair of pure imaginary eigenvalues λ=±iω1(τn1). System (2.3) undergoes a Hopf bifurcation at τn1 provided dRe(λ(τ))/dτ|τ=τn1≠0.Proof.
AssumeD=Δ=1 in system (2.3).Case 1.
Supposeτ∈I1. By Theorem 4.5, ωj(τ)⩾0 for j=1,2. By Lemma 4.7, there exists ϵj>0 and θj(τ) such that ϵj⩽θj(τ)⩽π and θj(τ) satisfies (4.56). Assume that there exists a positive integer τnj∈I1 such that θj(τnj)+2nπ=τnjωj(τnj) for some integer n⩾0. Then system (4.52) has one solution (τnj,ωj(τnj)). Equation (4.4) has a pair of pure imaginary eigenvalues λ=±iωj(τnj).
In what follows, we show that the conditions required for a Hopf Bifurcation (see TheoremA.1 in the appendix) are satisfied by the linearization (3.3) of (2.3) at E+. In (A.1), choosing τ as the bifurcation parameter and letting
(4.70)D(τ,zt)=[100010001][z1(t)z2(t)z3(t)],L(τ,zt)=[-1-αx+-αs+0αx+-D+αs+-ky+-kx+00-Δ][z1(t)z2(t)z3(t)]+[0000000keΔτy+keΔτx+][z1(t-τ)z2(t-τ)z3(t-τ)],
the linearization (3.3) of (2.3) at E+ is of the form (A.1). Taking a to be any positive real number and b=1/2, hypothesis (S1) in the Hopf Bifurcation Theorem holds, since
(4.71)|det(∑k=0∞Ak(α)e-λrk(α))|=|det[100010001]|=1⩾12,|det(∑k=0∞Ak(α)e-λrk(α)+∫-10A(α,θ)eλθdθ)|=|det[100010001]|=1⩾12
for all τ∈ℝ and |Reλ|<a.
The characteristic equation (4.4) of (3.3) at E+ has a pair of pure imaginary eigenvalues λ=±ω1(τnj) and no other root of (4.4) is an integral multiple of ±ω1(τnj). Hence the hypothesis (S2) in the Hopf Bifurcation Theorem holds. Therefore, (2.3) undergoes a Hopf bifurcation at E+ when τ=τnj provided dRe(λ(τ))/dτ|τ=τnj≠0.Case 2.
Supposeτ∈I2. By Theorem 4.5, only ω1(τ)>0. By Lemma 4.8, there exists ϵ>0 and θ1(τ) such that ϵ⩽θ1(τ)<π and θ1(τ) satisfies (4.56). Assume there exists τn1∈I2 such that θ1(τn1)+2nπ=τn1ω1(τn1) for some integer n⩾0. Then system (4.52) has one solution (τn1,ω1(τn1)). Equation (4.4) has a pair of pure imaginary eigenvalues λ=±iω1(τn1). The rest of the proof is similar to that of Case 1 when j=1.Corollary 4.10.
Consider system (2.3) with D=Δ=1. (1)
Supposeα>(7+35)/2, k>α/(α4-1), and τ∈I1 given by (4.22). For τ∈I1, ωj(τ) is nonnegative and there exists ϵj>0 and θj(τ) such that ϵj⩽θj(τ)⩽π and θj(τ) satisfies (4.56) for j=1,2. If there exists a positive integer nj⩾0 such that minτ∈I1τωj(τ)⩽2njπ and maxτ∈I1τωj(τ)>(2nj+1)π, then θj(τ)+2njπ intersects τωj(τ) at least once at some τnjj∈I1. System (2.3) undergoes a Hopf bifurcation at τnjj provided dRe(λ(τ))/dτ|τ=τnjj≠0.(2)
Supposeα>1, k>α(1/(1/4+1/16+1/(2α))-1)-1, and τ∈I2 defined in (4.41). For τ∈I2, only ω1(τ) is positive. There exists ϵ>0 and θ1(τ) such that ϵ⩽θ1(τ)<π and θ1(τ) satisfies (4.56) for j=1. If there exists a positive integer N⩾0 such that maxτ∈I2τω1(τ)>(2N+1)π, then for any 0⩽n⩽N, θ1(τ)+2nπ intersects τω1(τ) at least once at some τn1∈I2. System (2.3) undergoes a Hopf bifurcation at τn1 provided dRe(λ(τ))/dτ|τ=τn1≠0.Proof.
AssumeD=Δ=1 in system (2.3).Case 1.
Supposeτ∈I1. By Theorem 4.5, ωj(τ)⩾0 for j=1,2. By Lemma 4.7, there exists ϵj>0 and θj(τ) such that ϵj⩽θj(τ)⩽π and θj(τ) satisfies (4.56). Assume that there exists a positive integer nj⩾0 such that minτ∈I1τωj(τ)⩽2njπ and maxτ∈I1τωj(τ)>(2nj+1)π. For such nj,
(4.72)minτ∈I1τωj(τ)<ϵj+2njπ⩽θj(τ)+2njπ⩽(2nj+1)π<maxτ∈I1τωj(τ).
By the Mean Value Theorem, there exists τnjj∈I1 such that θj(τnjj)+2njπ=τnjjωj(τnjj). By Theorem 4.9. Case 1, the conclusion follows.Case 2.
Supposeτ∈I2. By Theorem 4.5, only ω1(τ)>0. By Lemma 4.8, there exists ϵ>0 and θ1(τ) such that ϵ⩽θ1(τ)<π and θ1(τ) satisfies (4.56). Assume that there exists a positive integer N⩾0 such that maxτ∈I2τω1(τ)>(2N+1)π. By (4.41), 0∈I2. Therefore minτ∈I2τω1(τ)=0. For 0⩽n⩽N,
(4.73)minτ∈I2τω1(τ)<ϵ+2nπ⩽θi(τ)+2nπ⩽(2n+1)π<maxτ∈I2τω1(τ).
By the Mean Value Theorem, there exists τn1∈I2 such that θ1(τn1)+2nπ=τn1ω1(τn1). By Theorem 4.9. Case 2, the conclusion follows.Corollary 4.11.
Consider system (2.3) with D=Δ=1.  Assume α>(7+35)/2 and k>α/(α4-1). If 1/4+1/16+1/(2α)>k/(α+k), then I1=[0,ln((k/α)(α4-1))], where I1 was defined in (4.22). For any τ∈I1, ωj(τ) is nonnegative and there exists ϵj>0 and θj(τ) such that ϵj⩽θj(τ)⩽π and θj(τ) satisfies (4.56) for j=1,2. If there exists a positive integer Nj⩾0 (j=1,2) such that maxτ∈I1τωj(τ)>(2Nj+1)π, then for any 0⩽n⩽Nj, θj(τ)+2nπ intersects τωj(τ) at least once at some τnj∈I1. System (2.3) undergoes a Hopf bifurcation at τnj provided dRe(λ(τ))/dτ|τ=τnj≠0.Proof.
Assume1/4+1/16+1/(2α)>k/(α+k). Then  ln((k/α)(1/(1/4+1/16+1/(2α))-1))<0. By (4.22),
(4.74)I1=[max{0,ln(kα(11/4+1/16+1/(2α)-1))},ln(k(α4-1)α)]=[0,ln(kα(α4-1))].
For any τ∈I1, by Theorem 4.5, ωj(τ)⩾0 for j=1,2. By Lemma 4.7, there exists ϵj>0 and θj(τ) such that ϵj⩽θj(τ)⩽π and θj(τ) satisfies (4.56). Noting 0∈I1, minτ∈I1τωj(τ)=0. Assume there exists a positive integer Nj⩾0 (j=1,2) such that maxτ∈I1τωj(τ)>(2Nj+1)π. For any 0⩽n⩽Nj, minτ∈I1τωj(τ)=0⩽2nπ and maxτ∈I1τωj(τ)>(2Nj+1)π⩾(2n+1)π. By Corollary 4.10, the conclusion follows.


## 5. Numerical Results

This section includes bifurcation diagrams involving the interior equilibriumE+ and numerical simulations of periodic solutions of the predator-prey model in the chemostat.


### 5.1. Variation of Eigenvalues

To study the stability switches ofE+, DDEBIFTOOL (see [11, 12]) was chosen to illustrate how the real part of the eigenvalues of (4.2) changes as parameters α and τ vary.First fix parametersD=Δ=1, k=24, and τ=0.5. Taking α as the bifurcation parameter and varying it from 0 to 10, the real part of the eigenvalues with largest real part of (4.2) was plotted in Figure 1. At α≈1.15 and α≈1.5, there is either a zero eigenvalue or a pair of pure imaginary roots. For α∈(1.15,1.5), all eigenvalues have negative real parts. For example, taking α=1.3, Figure 2(a) shows that the eigenvalues of (4.2) with largest real parts (the ones in the circle) have negative real parts. Note that due to the scaling, the eigenvalues in the circle seem to be indistinguishable from zero. But in fact, they are a pair of complex eigenvalues with real parts slightly less than zero. DDEBIFTOOL can keep track of the occurrence of a pair of pure imaginary eigenvalues as α varies in the neighborhood of α=1.5. Figure 2(b) clearly shows that there is a pair of pure imaginary eigenvalues. Hence, Hopf bifurcation is possible. Note that by continuation, the pair of eigenvalues with largest real parts in Figure 2(a) for α=1.3 becomes the pair of pure imaginary eigenvalues in Figure 2(b) for α≈1.5.Figure 1
Variation of the largest real part of the eigenvalues as the bifurcation parameterα is varied. At α≈1.15 and 1.5, the largest real part crosses zero and it seems that there is a zero eigenvalue for α≈1.15 and a pair of pure imaginary eigenvalues for α≈1.5. The largest real part becomes positive as α increases through 1.5. But as α increases further, for α≈17, the largest real part crosses zero again and remains negative thereafter. There is a second Hopf bifurcation at α≈17. This is consistent with what is observed in Figure 3 when τ=0.5 and α varies from 0 to 30. Parameters are D=Δ=1, k=24, and τ=0.5.Eigenvalues with the largest real parts of the characteristic equation (4.2) at E+. Parameters are the same as in Figure 1 except α=1.3 for the TOP and α≈1.5 for the BOTTOM graph. Due to the scaling, the eigenvalues in the circle in the TOP graph seem indistinguishable from zero. In fact, they are a pair of complex eigenvalues with real parts slightly less than zero. As α varies from 1.3 to 1.5, the pair of complex eigenvalues with largest real part becomes a pair of pure imaginary roots in the BOTTOM graph. The eigenvalue with the second largest real part remains equal to -1. This is consistent with our analytical results that showed that (4.2) has a constant eigenvalue -1 when D=Δ=1.Finally fix all parameters as before and vary bothτ and α. In Figure 3, we plot the Hopf bifurcation diagram in α and τ parameter space. The curve at the left upper corner is τ=τc. For any pair (α,τ) below that curve, a coexistence equilibrium E+ exists (i.e., all components are positive). For any pair (α,τ) on the closed curve, there is a Hopf bifurcation. Inside the closed curve, there is a periodic solution surrounding E+. For any (α,τ) outside the closed curve and below τ=τc, the coexistence equilibrium E+ is stable.Figure 3
The two-parameter bifurcation diagram ofE+ in τ and α parameter space. Parameters are the same as in Figure 1 (i.e., D=Δ=1,k=24) except now both τ and α are allowed to vary. For any pair (α,τ) on the closed curve, there is a Hopf bifurcation of E+. Inside the closed curve, there is a periodic solution surrounding E+. For any (α,τ) outside the closed curve and below τ=τc, the coexistence equilibrium E+ is stable.


### 5.2. Simulations Demonstrating Hopf Bifurcations

In this section, we illustrate Theorem4.9 for system (2.3). Take D=Δ=1 and let τ vary. We choose parameters α=100 and k=100 for Case 1 (see Figures 4–12), and α=2 and k=20 for Case 2 (see Figures 13–18).Critical value of delayτ at which a Hopf Bifurcation occurs for D=Δ=1,k=100,α=100.
(a)(b)If the slope ofSnj(τ) is nonzero at τ when Snj(τ)=0 (j=1,2), the transversality condition holds and there is a Hopf bifurcation. The other parameters are the same as for Figure 4.
(a)(b)EquilibriumE+(τ) is stable when τ=0.02<τ01. The other parameters are the same as for Figure 4.
(a)(b)(c)Time series of a periodic solution, forτ=0.03. The other parameters are the same as for Figure 4.
(a)(b)(c)Figure 8
The trajectory in phase space of the periodic solution in Figure7 for τ=0.03.Time series of a periodic solution forτ=0.4. The other parameters are the same as for Figure 4.
(a)(b)(c)Figure 10
The trajectory in phase space of the periodic solution shown in Figure9.The periodic solution disappears at the secondary Hopf bifurcation atτ02≈0.48 and E+ regains stability. In this figure τ=0.5>τ02. The other parameters are the same as for Figure 4.
(a)(b)(c)Figure 12
Bifurcation diagram as the delay varies. The maximum and minimum amplitude of they-component of the solution is plotted on the ordinate axis. Parameters are the same as for Figure 4.Figure 13
Intersections indicate critical values of the delay at which Hopf bifurcations occur. Parameters areD=Δ=1,k=20,α=2.Figure 14
Verification of the transversality condition required for Hopf bifurcation. Parameters are the same as in Figure13.EquilibriumE+(τ) is stable when τ=0.15<τ0,11. The other parameters are the same as in Figure 13.
(a)(b)(c)Time series of a solution with constant initial datas(0)=0.87, x(t)=0.077, and y(t)=0.048 for t∈[-0.3,0], that approaches a stable periodic solution as time increases. In this figure, τ=0.3. The other parameters are the same as in Figure 13.
(a)(b)(c)Figure 17
The attracting periodic solution shown in Figure16 in phase space. Note that τ=0.3>τ0,11.The periodic solution disappears andE+(τ) regains its stability when τ>τ02. In this figure, τ=0.74. The other parameters are the same as in Figure 13.
(a)(b)(c)Case 1.
Note thatI1 is given by (4.22). Since ln((k/α)(1/(1/4+1/16+1/(2α))-1))=-0.03,
(5.1)max{0,ln(kα(11/4+1/16+1/(2α)-1))}=0.
Also, ln(k(α4-1)/α)≈0.77. Therefore I1≈[0,0.77]. By Theorem 4.9, ωj(τ) is positive and θj(τ) satisfies (4.56) for any τ∈I1 and j=1,2.
Figure4 shows that θj(τ) intersects τωj(τ) at some τ0j with τ01≈0.022 and τ02≈0.48. We see that θj(τ)+2nπ has no intersection with τωj(τ) for n⩾2 and j=1,2. By Theorem 4.9, (4.4) has two distinct pairs of pure imaginary eigenvalues λ=±iωj(τ0j). Next we need to check if Re(dλ(τ)/dτ)|τ=τ0j≠0.
As in Beretta and Kuang [13], we can define
(5.2)Snj(τ)=τ-θj(τ)+2nπωj(τ),forj=1,2,n=0,1,2….
Any zero τnj of Snj(τ) is an intersection of θj(τ)+2nπ and τωj(τ) and vice versa. By (4.10) in Beretta and Kuang [13] and noting that (q2-p2(τ)+2α2(τ))2-4(α2(τ)-c2(τ))>0, we have the relation
(5.3)sign{Redλ(τ)dτ|τ=τnj}=±sign{(q2-p2(τ)+2α2(τ))2-4(α2(τ)-c2(τ))}sign{dSnj(τ)dτ|τ=τnj}=±sign{dSnj(τ)dτ|τ=τnj},
where we take + for j=1 and - for j=2.
From Figure5(a), it is observed that Sn1 has only one zero τ01≈0.022 at n=0 with sign{dS01(τ)/dτ|τ=τ01}>0. Hence, sign{Re(dλ(τ)/dτ)|τ=τ01}>0. By Theorem 4.9, system (2.3) undergoes a Hopf bifurcation at τ01. Similarly from Figure 5(b), Sn2 has only one zero τ02≈0.48 for n=0 and (2.3) undergoes a Hopf bifurcation at τ02.
Next we used MATLAB to simulate solutions of model (2.3) for several values of τ. For each fixed delay τ, we chose initial data s(t)=s+(τ)-0.01, x(t)=x+(τ)+0.01, and y(t)=y+(τ)+0.001 for t∈[-τ,0]. From Figure 6, we can see that the equilibrium E+ is stable if τ=0.02<τ01. As delay τ increases past τ01≈0.022, where a Hopf bifurcation occurs, a pair of complex eigenvalues of (4.4) enters the right-half plane. The equilibrium E+ loses its stability and a periodic solution bifurcates from E+ (see Figures 7 and 8). As we increase the delay further to τ=0.4<τ02, the periodic solution still exists and remains stable (see Figures 9 and 10). However, as the delay τ increases further, past τ02, the stable periodic solution disappears in a second Hopf bifurcation, and E+ regains stability (see Figure 11). We provide a bifurcation diagram illustrating the change in dynamics as τ varies (see Figure 12). For any τ∈(τ01,τ02), there is an orbitally asymptotically stable periodic solution.Case 2.
Takek=20 and α=2. For such parameters, I2≈[0,0.85]. By Theorem 4.9, ω1(τ) is positive and θ1(τ) satisfies (4.56) for j=1 and τ∈I2. Figure 13 shows that θ1(τ) intersects τω1(τ) twice. To distinguish these intersections, denote them as τ0,11≈0.2 and τ0,21≈0.7. On the other hand, θ1(τ)+2π has no intersection with τω1(τ).
From Figure14,
(5.4)sign{dS01(τ)dτ|τ=τ0,11}>0,sign{dS11(τ)dτ|τ=τ0,21}<0.
By (5.3),
(5.5)sign{Redλ(τ)dτ|τ=τ0,11}>0,sign{Redλ(τ)dτ|τ=τ0,21}<0.
By Theorem 4.5, system (2.3) undergoes a Hopf bifurcation at τ=τ0,11 and at τ=τ0,21. For τ less than τ0,11, the equilibrium E+ is asymptotically stable (see Figure 15). For τ greater than τ0,11, but less than τ0,21, there is an orbitally asymptotically stable periodic solution surrounding the equilibrium E+ (see Figures 16 and 17). At τ=τ0,21, there is a second Hopf bifurcation, where the periodic solution coalesces with E+. For τ>τ0,21, the periodic orbit no longer exists and E+ regains stability (see Figure 18) until it disappears when τ>τc.


## 5.1. Variation of Eigenvalues

To study the stability switches ofE+, DDEBIFTOOL (see [11, 12]) was chosen to illustrate how the real part of the eigenvalues of (4.2) changes as parameters α and τ vary.First fix parametersD=Δ=1, k=24, and τ=0.5. Taking α as the bifurcation parameter and varying it from 0 to 10, the real part of the eigenvalues with largest real part of (4.2) was plotted in Figure 1. At α≈1.15 and α≈1.5, there is either a zero eigenvalue or a pair of pure imaginary roots. For α∈(1.15,1.5), all eigenvalues have negative real parts. For example, taking α=1.3, Figure 2(a) shows that the eigenvalues of (4.2) with largest real parts (the ones in the circle) have negative real parts. Note that due to the scaling, the eigenvalues in the circle seem to be indistinguishable from zero. But in fact, they are a pair of complex eigenvalues with real parts slightly less than zero. DDEBIFTOOL can keep track of the occurrence of a pair of pure imaginary eigenvalues as α varies in the neighborhood of α=1.5. Figure 2(b) clearly shows that there is a pair of pure imaginary eigenvalues. Hence, Hopf bifurcation is possible. Note that by continuation, the pair of eigenvalues with largest real parts in Figure 2(a) for α=1.3 becomes the pair of pure imaginary eigenvalues in Figure 2(b) for α≈1.5.Figure 1
Variation of the largest real part of the eigenvalues as the bifurcation parameterα is varied. At α≈1.15 and 1.5, the largest real part crosses zero and it seems that there is a zero eigenvalue for α≈1.15 and a pair of pure imaginary eigenvalues for α≈1.5. The largest real part becomes positive as α increases through 1.5. But as α increases further, for α≈17, the largest real part crosses zero again and remains negative thereafter. There is a second Hopf bifurcation at α≈17. This is consistent with what is observed in Figure 3 when τ=0.5 and α varies from 0 to 30. Parameters are D=Δ=1, k=24, and τ=0.5.Eigenvalues with the largest real parts of the characteristic equation (4.2) at E+. Parameters are the same as in Figure 1 except α=1.3 for the TOP and α≈1.5 for the BOTTOM graph. Due to the scaling, the eigenvalues in the circle in the TOP graph seem indistinguishable from zero. In fact, they are a pair of complex eigenvalues with real parts slightly less than zero. As α varies from 1.3 to 1.5, the pair of complex eigenvalues with largest real part becomes a pair of pure imaginary roots in the BOTTOM graph. The eigenvalue with the second largest real part remains equal to -1. This is consistent with our analytical results that showed that (4.2) has a constant eigenvalue -1 when D=Δ=1.Finally fix all parameters as before and vary bothτ and α. In Figure 3, we plot the Hopf bifurcation diagram in α and τ parameter space. The curve at the left upper corner is τ=τc. For any pair (α,τ) below that curve, a coexistence equilibrium E+ exists (i.e., all components are positive). For any pair (α,τ) on the closed curve, there is a Hopf bifurcation. Inside the closed curve, there is a periodic solution surrounding E+. For any (α,τ) outside the closed curve and below τ=τc, the coexistence equilibrium E+ is stable.Figure 3
The two-parameter bifurcation diagram ofE+ in τ and α parameter space. Parameters are the same as in Figure 1 (i.e., D=Δ=1,k=24) except now both τ and α are allowed to vary. For any pair (α,τ) on the closed curve, there is a Hopf bifurcation of E+. Inside the closed curve, there is a periodic solution surrounding E+. For any (α,τ) outside the closed curve and below τ=τc, the coexistence equilibrium E+ is stable.


## 5.2. Simulations Demonstrating Hopf Bifurcations

In this section, we illustrate Theorem4.9 for system (2.3). Take D=Δ=1 and let τ vary. We choose parameters α=100 and k=100 for Case 1 (see Figures 4–12), and α=2 and k=20 for Case 2 (see Figures 13–18).Critical value of delayτ at which a Hopf Bifurcation occurs for D=Δ=1,k=100,α=100.
(a)(b)If the slope ofSnj(τ) is nonzero at τ when Snj(τ)=0 (j=1,2), the transversality condition holds and there is a Hopf bifurcation. The other parameters are the same as for Figure 4.
(a)(b)EquilibriumE+(τ) is stable when τ=0.02<τ01. The other parameters are the same as for Figure 4.
(a)(b)(c)Time series of a periodic solution, forτ=0.03. The other parameters are the same as for Figure 4.
(a)(b)(c)Figure 8
The trajectory in phase space of the periodic solution in Figure7 for τ=0.03.Time series of a periodic solution forτ=0.4. The other parameters are the same as for Figure 4.
(a)(b)(c)Figure 10
The trajectory in phase space of the periodic solution shown in Figure9.The periodic solution disappears at the secondary Hopf bifurcation atτ02≈0.48 and E+ regains stability. In this figure τ=0.5>τ02. The other parameters are the same as for Figure 4.
(a)(b)(c)Figure 12
Bifurcation diagram as the delay varies. The maximum and minimum amplitude of they-component of the solution is plotted on the ordinate axis. Parameters are the same as for Figure 4.Figure 13
Intersections indicate critical values of the delay at which Hopf bifurcations occur. Parameters areD=Δ=1,k=20,α=2.Figure 14
Verification of the transversality condition required for Hopf bifurcation. Parameters are the same as in Figure13.EquilibriumE+(τ) is stable when τ=0.15<τ0,11. The other parameters are the same as in Figure 13.
(a)(b)(c)Time series of a solution with constant initial datas(0)=0.87, x(t)=0.077, and y(t)=0.048 for t∈[-0.3,0], that approaches a stable periodic solution as time increases. In this figure, τ=0.3. The other parameters are the same as in Figure 13.
(a)(b)(c)Figure 17
The attracting periodic solution shown in Figure16 in phase space. Note that τ=0.3>τ0,11.The periodic solution disappears andE+(τ) regains its stability when τ>τ02. In this figure, τ=0.74. The other parameters are the same as in Figure 13.
(a)(b)(c)Case 1.
Note thatI1 is given by (4.22). Since ln((k/α)(1/(1/4+1/16+1/(2α))-1))=-0.03,
(5.1)max{0,ln(kα(11/4+1/16+1/(2α)-1))}=0.
Also, ln(k(α4-1)/α)≈0.77. Therefore I1≈[0,0.77]. By Theorem 4.9, ωj(τ) is positive and θj(τ) satisfies (4.56) for any τ∈I1 and j=1,2.
Figure4 shows that θj(τ) intersects τωj(τ) at some τ0j with τ01≈0.022 and τ02≈0.48. We see that θj(τ)+2nπ has no intersection with τωj(τ) for n⩾2 and j=1,2. By Theorem 4.9, (4.4) has two distinct pairs of pure imaginary eigenvalues λ=±iωj(τ0j). Next we need to check if Re(dλ(τ)/dτ)|τ=τ0j≠0.
As in Beretta and Kuang [13], we can define
(5.2)Snj(τ)=τ-θj(τ)+2nπωj(τ),forj=1,2,n=0,1,2….
Any zero τnj of Snj(τ) is an intersection of θj(τ)+2nπ and τωj(τ) and vice versa. By (4.10) in Beretta and Kuang [13] and noting that (q2-p2(τ)+2α2(τ))2-4(α2(τ)-c2(τ))>0, we have the relation
(5.3)sign{Redλ(τ)dτ|τ=τnj}=±sign{(q2-p2(τ)+2α2(τ))2-4(α2(τ)-c2(τ))}sign{dSnj(τ)dτ|τ=τnj}=±sign{dSnj(τ)dτ|τ=τnj},
where we take + for j=1 and - for j=2.
From Figure5(a), it is observed that Sn1 has only one zero τ01≈0.022 at n=0 with sign{dS01(τ)/dτ|τ=τ01}>0. Hence, sign{Re(dλ(τ)/dτ)|τ=τ01}>0. By Theorem 4.9, system (2.3) undergoes a Hopf bifurcation at τ01. Similarly from Figure 5(b), Sn2 has only one zero τ02≈0.48 for n=0 and (2.3) undergoes a Hopf bifurcation at τ02.
Next we used MATLAB to simulate solutions of model (2.3) for several values of τ. For each fixed delay τ, we chose initial data s(t)=s+(τ)-0.01, x(t)=x+(τ)+0.01, and y(t)=y+(τ)+0.001 for t∈[-τ,0]. From Figure 6, we can see that the equilibrium E+ is stable if τ=0.02<τ01. As delay τ increases past τ01≈0.022, where a Hopf bifurcation occurs, a pair of complex eigenvalues of (4.4) enters the right-half plane. The equilibrium E+ loses its stability and a periodic solution bifurcates from E+ (see Figures 7 and 8). As we increase the delay further to τ=0.4<τ02, the periodic solution still exists and remains stable (see Figures 9 and 10). However, as the delay τ increases further, past τ02, the stable periodic solution disappears in a second Hopf bifurcation, and E+ regains stability (see Figure 11). We provide a bifurcation diagram illustrating the change in dynamics as τ varies (see Figure 12). For any τ∈(τ01,τ02), there is an orbitally asymptotically stable periodic solution.Case 2.
Takek=20 and α=2. For such parameters, I2≈[0,0.85]. By Theorem 4.9, ω1(τ) is positive and θ1(τ) satisfies (4.56) for j=1 and τ∈I2. Figure 13 shows that θ1(τ) intersects τω1(τ) twice. To distinguish these intersections, denote them as τ0,11≈0.2 and τ0,21≈0.7. On the other hand, θ1(τ)+2π has no intersection with τω1(τ).
From Figure14,
(5.4)sign{dS01(τ)dτ|τ=τ0,11}>0,sign{dS11(τ)dτ|τ=τ0,21}<0.
By (5.3),
(5.5)sign{Redλ(τ)dτ|τ=τ0,11}>0,sign{Redλ(τ)dτ|τ=τ0,21}<0.
By Theorem 4.5, system (2.3) undergoes a Hopf bifurcation at τ=τ0,11 and at τ=τ0,21. For τ less than τ0,11, the equilibrium E+ is asymptotically stable (see Figure 15). For τ greater than τ0,11, but less than τ0,21, there is an orbitally asymptotically stable periodic solution surrounding the equilibrium E+ (see Figures 16 and 17). At τ=τ0,21, there is a second Hopf bifurcation, where the periodic solution coalesces with E+. For τ>τ0,21, the periodic orbit no longer exists and E+ regains stability (see Figure 18) until it disappears when τ>τc.



---

*Source: 287969-2010-03-31.xml*