A Predator-Prey Model in the Chemostat with Time Delay

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.

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) Here 𝑠(𝑡) represents the concentration of the growth limiting nutrient, 𝑥(𝑡) the density of the prey population, and 𝑦(𝑡) the density of the predator population. Parameter 𝑠⁰ denotes the concentration of the growth limiting nutrient in the feed vessel, 𝐷₀ the dilution rate, 𝜂(𝜉) the growth yield constant, 𝐷(Δ) the sum of the dilution rate 𝐷₀ and the natural species specific death rate of the prey (predator) population, respectively. Here 𝑓(𝑠) denotes the functional response of the prey population on the nutrient and 𝑔(𝑥) 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 𝑓(𝑠) and 𝑔(𝑥) 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 𝐷₀=𝐷=Δ. 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 𝑠⁰. 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 𝑓(𝑠) is modelled by Holling type I or II and 𝑔(𝑥) 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 𝑓(𝑠) is given by Holling type I, if 𝑔(𝑥) 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 𝑓(𝑠) and 𝑔(𝑥), 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) For 𝑡∈[−𝜏,0],

𝑠(0)=𝑠₀∈intℝ₊,(𝑥(𝑡),𝑦(𝑡))=(𝜙,𝜓)∈ℂ[−𝜏,0],intℝ²⁺.(1.3) Here variables 𝑠(𝑡), 𝑥(𝑡), 𝑦(𝑡), and parameters 𝑠⁰, 𝐷₀, 𝜂, 𝜉, 𝐷₀, 𝐷, and Δ have the same interpretation as for model (1.1). Note therefore that 𝐷⩾𝐷₀, and Δ⩾𝐷₀. The additional parameter 𝜏 is a nonnegative constant modelling the time required for the conversion process. Hence, 𝑒^−Δ𝜏𝑦(𝑡−𝜏) represents the concentration of the predator population in the growth chamber at time 𝑡 that were available at time 𝑡−𝜏 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 functions 𝑓(𝑠) and 𝑔(𝑠) are of Holling type I form, that is, 𝑓(𝑠)=𝛼𝑠 (𝛼>0) and 𝑔(𝑥)=𝑘𝑥 (𝑘>0). System (1.2) reduces to

𝑠̇(𝑡)=𝑠⁰−𝑠(𝑡)𝐷₀−𝛼𝑥(𝑡)𝑠(𝑡)𝑥̇(𝑡)=𝑥(𝑡)(−𝐷+𝛼𝑠(𝑡))−𝑘𝜂𝑥(𝑡)𝑦,(𝑡)𝑦̇(𝑡)=−Δ𝑦(𝑡)+𝑘𝑒_−Δ𝜏𝜉,𝑦(𝑡−𝜏)𝑥(𝑡−𝜏).𝑡>0,(2.1) Introducing the following change of variables gives:

𝑡̆=𝐷₀𝑡,𝑠̆𝑡̆=𝑠(𝑡)𝑠⁰,𝑥̆𝑡̆=𝑥(𝑡)𝜏̆=𝐷₀𝑠⁰𝜂,𝑦̆𝑡̆=𝜉𝑦𝑠(⁰𝑡𝜂),𝜏,𝐷̆=𝐷𝐷₀,Δ̆=Δ𝐷₀,𝑘̆=𝑘𝑠⁰𝜂𝐷₀,𝛼̆=𝛼𝑠⁰𝐷₀d𝑠̆𝑡̆,d𝑡̆=𝑠1⁰d𝑠(𝑡)d𝑡dd𝑡𝑡̆=𝑠⁰1𝐷₀d𝑠(𝑡)=1d𝑡𝑠⁰𝐷₀𝑠⁰−𝑠(𝑡)𝐷₀−𝛼𝑥(𝑡)𝑠(𝑡)=1−𝑠(𝑡)𝜂𝑠⁰−𝛼𝑠⁰𝐷₀𝑥(𝑡)𝑠⁰𝜂𝑠𝑠(⁰𝑡)d𝑥̆𝑡̆=1−𝑠̆𝑡̆−𝛼̆𝑥̆𝑡̆𝑠̆𝑡̆,d𝑡̆=𝑠1⁰𝜂d𝑥d𝑡(𝑡)dd𝑡𝑡̆=𝑠⁰𝜂1𝐷₀d𝑥(𝑡)=1d𝑡=𝑠𝑥⁰(𝜂𝑡𝐷)₀𝑥(𝑡)(−𝐷+𝛼𝑠(𝑡))−𝑘𝑥(𝑡𝜉)𝑦(𝑡)𝑠⁰𝜂−𝐷𝐷₀+𝛼𝑠⁰𝐷₀𝑠(𝑡)𝑠⁰−𝑘𝑠⁰𝜂𝐷₀𝑥(𝑡)d𝑦̆𝑡̆=𝑥̆𝑡̆−𝐷̆+𝛼̆𝑠̆𝑡̆−𝑘̆𝑥̆𝑡̆𝑦̆𝑡̆𝑠,⁰𝜂𝑠𝑦⁰(𝜂𝑡𝜉)d𝑡̆=𝑠⁰1𝜂𝜉d𝑦d(𝑡𝑡)dd𝑡𝑡̆=𝑠⁰𝜂1𝜉𝐷₀d𝑦(𝑡)=1d𝑡𝑠⁰𝜂𝜉𝐷₀=−Δ𝑦(𝑡)−Δ𝑦(𝑡)+𝑘𝑒^−Δ𝜏𝑦(𝑡−𝜏)𝑥(𝑡−𝜏)𝑠⁰𝜂𝜉𝐷₀+𝑘𝑒^−Δ𝜏𝑠⁰𝜂𝜉𝐷₀=−Δ𝑦(𝑡−𝜏)𝑥(𝑡−𝜏)𝐷₀𝑦(𝑡)𝑠⁰𝜂𝜉+𝑘𝑠⁰𝜂𝐷₀𝑒^{−(Δ/𝐷0)𝐷0𝜏}𝑦(𝑡−𝜏)=−Δ̆𝑦̆𝑡̆+𝑘̆𝑒_{−̆Δ̆𝜏}𝑠⁰𝜂𝜉𝑥(𝑠𝑡⁰−𝜂𝜏)𝑦̆𝑡̆−𝜏̆𝑥̆𝑡̆−𝜏̆.(2.2) With this change of variables, omitting the ̆‌’s for convenience, system (2.1) becomes

𝑥𝑦𝑠̇̇̇(((𝑡𝑡𝑡)))===𝑥1−Δ(−𝑡)𝑦𝑠((−(𝑡𝑡)𝐷)+−+𝑘𝛼𝑒𝛼𝑥𝑠(𝑡(_−Δ𝜏)𝑡)𝑠)(−𝑡)𝑘,𝑦(𝑡)𝑥(𝑡),𝑦(𝑡−𝜏)𝑥(𝑡−𝜏),(2.3) where Δ⩾1 and 𝐷⩾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 𝐷>𝛼, the model has only one equilibrium point (1,0,0) and it is globally asymptotically stable. If 𝐷<𝛼 and 1−𝐷/𝛼−Δ𝐷/𝑘<0, the model has a second equilibrium point (𝐷/𝛼,(𝛼−𝐷)/𝛼𝐷,0) and it is globally asymptotically stable. When 1−𝐷/𝛼−Δ𝐷/𝑘>0, the model has a third equilibrium point (𝑘/(𝑘+𝛼Δ),Δ/𝑘,𝛼/(𝑘+𝛼𝐷)−𝐷/𝑘) and it is the global attractor. Therefore, model (2.3) has no periodic solutions when the time delay is ignored. If 𝑔(𝑥) 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 (𝑠₀,𝜙(𝜃),𝜓(𝜃))∈intℝ₊×ℂ([−𝜏,0],intℝ²⁺), then there exists a unique solution (𝑠(𝑡),𝑥(𝑡),𝑦(𝑡)) of (2.3) passing through (𝑠₀,𝜙(𝜃),𝜓(𝜃)) with 𝑠(𝑡)>0, 𝑥(𝑡)>0 and 𝑦(𝑡)>0 for 𝑡∈[0,∞). The solution is bounded. In particular, given any 𝜖₀>0, 𝑥(𝑡)<1+𝜖₀ for all sufficiently large 𝑡.

Proof. For 𝑡∈[0,𝜏], one has 𝑡−𝜏∈[−𝜏,0],𝑥(𝑡−𝜏)=𝜙(𝑡−𝜏), and 𝑦(𝑡−𝜏)=𝜓(𝑡−𝜏). System (2.3) becomes 𝑥𝑦𝑠̇̇̇(((𝑡𝑡𝑡)))===𝑥1−Δ(−𝑡)𝑦𝑠((−(𝑡𝑡)𝐷)+−+𝑘𝛼𝑒𝛼𝑥𝑠(𝑡(_−Δ𝜏)𝑡)𝑠)(−𝑡)𝑘,𝑦(𝑡)𝑥(𝑡),𝜙(𝑡−𝜏)𝜓(𝑡−𝜏),(2.4) a system of nonautonomous ordinary differential equations with initial conditions 𝑠(0)=𝑠₀, 𝑥(0)=𝜙(0), and 𝑦(0)=𝜓(0). Since the right-hand side of (2.4) is differentiable in both 𝑥 and 𝑦, 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 (𝑠₀,𝜙(𝜃),𝜓(𝜃)) is defined for all 𝑡⩾0.
Now we prove 𝑠(𝑡)>0 for all 𝑡>0. From the first equation of (2.3),
𝑠̇(𝑡)=1−𝑠(𝑡)−𝛼𝑥(𝑡)𝑠(𝑡).(2.5) Proceed using the method of contradiction. Suppose that there exists a first 𝑡^⋆ such that 𝑠(𝑡^⋆)=0 and 𝑠(𝑡)>0 for 𝑡∈[0,𝑡^∗). Then 𝑠(̇𝑡^⋆)⩽0. But from the first equation of (2.3) 𝑠̇𝑡^⋆=1−𝑠𝑡^⋆−𝛼𝑥𝑡^⋆𝑠𝑡^⋆=1>0,(2.6) a contradiction.
To prove 𝑥(𝑡)>0 for 𝑡∈[0,∞), assume there is a first 𝑡>0 such that 𝑥(𝑡)=0, and 𝑥(𝑡)>0 for 𝑡∈[0,𝑡). Divide both sides of the second equation of (2.3) by 𝑥(𝑡) and integrate from 0 to 𝑡, to obtain
𝑥𝑡=𝜙(0)exp^𝑡0(−𝐷+𝛼𝑠(𝑡)−𝑘𝑦(𝑡))d𝑡>0,(2.7) contradicting 𝑥(𝑡)=0.
To show that 𝑦(𝑡) is positive on [0,∞), suppose that there exists 𝑡^⋆>0 such that 𝑦(𝑡^⋆)=0, and 𝑦(𝑡)>0 for 𝑡∈[0,𝑡^∗). Then 𝑦(̇𝑡^⋆)⩽0. From the third equation of (2.3), we have
𝑦̇𝑡^⋆==𝑘−𝑒Δ𝑦_−Δ𝜏𝑡^⋆+𝑘𝑒^−Δ𝜏𝑦𝑡^⋆−𝜏𝑥𝑡^⋆−𝜏𝑦𝑡^⋆−𝜏𝑥𝑡^⋆−𝜏>0,(2.8) a contradiction.
To prove the boundedness of solutions, define
𝜔(𝑡)=𝑠(𝑡)+𝑥(𝑡)+𝑒^Δ𝜏𝑦(𝑡+𝜏)−1,for𝑡⩾0.(2.9) It follows that 𝜔̇(𝑡)=⩽11−−𝑠𝑠((𝑡𝑡))−−𝐷𝑥(𝑥𝑡)(𝑡−)−𝑒Δ_Δ𝜏𝑒^Δ𝜏𝑦(𝑡+𝜏)⩽−𝜔(𝑡),𝑦(𝑡+𝜏)(2.10) where the first inequality holds since 𝐷⩾1, Δ⩾1, 𝑥(𝑡)>0 and 𝑦(𝑡+𝜏)>0. It follows that 𝑠(𝑡)+𝑥(𝑡)+𝑒^Δ𝜏𝑦(𝑡+𝜏)⩽1+𝑠₀+𝑥(0)+𝑒^Δ𝜏𝑦(𝜏)−1𝑒^−𝑡⟶1as𝑡⟶∞.(2.11) Therefore, the solution (𝑠(𝑡),𝑥(𝑡),𝑦(𝑡)) is bounded, and given any 𝜖₀>0, 𝑥(𝑡)<1+𝜖₀ for all sufficiently large 𝑡.

3. Equilibria and Stability

Model (2.3) has three equilibrium points: 𝐸₁=(1,0,0),𝐸₂=(𝐷/𝛼,(𝛼−𝐷)/𝛼𝐷,0), and

𝐸₊=𝑠₊(𝜏),𝑥₊(𝜏),𝑦₊(𝜏)=11+(𝛼Δ/𝑘)𝑒^Δ𝜏,Δ𝑘𝑒^Δ𝜏,𝛼𝑘+𝛼Δ𝑒^Δ𝜏−𝐷𝑘.(3.1) We call 𝐸₁ the washout equilibrium, 𝐸₂ the single species equilibrium, and 𝐸₊ the coexistence equilibrium. For the sake of biological significance, 𝐸₊ exists (distinct from 𝐸₂) if and only if its third coordinate 𝑦₊(𝜏)=(𝛼𝑠₊(𝜏)−𝐷)/𝑘>0, that is, 𝑠₊(𝜏)>𝐷/𝛼, or equivalently, 𝜏 lies between 0 and 𝜏_𝑐, where

𝜏_𝑐=1ΔlnΔ𝑘𝐷1−𝛼1.(3.2) Note that if (𝑘/Δ)(1/𝐷−1/𝛼)⩽1, the equilibrium 𝐸₊ does not exist for any 𝜏 (⩾0), and if (𝑘/Δ)(1/𝐷−1/𝛼)=1, then 𝐸₊=𝐸₂.

The linearization of (2.3) about an equilibrium (𝑠,𝑥,𝑦) is given by

⎡⎢⎢⎢⎣𝑧̇¹𝑧̇₂(𝑡)𝑧̇₃(𝑡)(𝑡)⎤⎥⎥⎥⎦=⎡⎢⎢⎢⎣−1𝛼−0𝑥𝛼𝑥−𝐷+−𝛼0𝛼𝑠𝑠−𝑘𝑦−−0𝑘Δ𝑥⎤⎥⎥⎥⎦⎡⎢⎢⎢⎣𝑧¹𝑧(𝑡)₂𝑧(𝑡)₃+⎡⎢⎢⎢⎣000𝑘𝑒00_Δ𝜏00(𝑡)⎤⎥⎥⎥⎦𝑦𝑘𝑒^Δ𝜏𝑥⎤⎥⎥⎥⎦⎡⎢⎢⎢⎣𝑧¹𝑧(𝑡−𝜏)₂𝑧(𝑡−𝜏)₃(𝑡−𝜏)⎤⎥⎥⎥⎦.(3.3) The associated characteristic equation is given by

det⎡⎢⎢⎢⎣−1−𝛼0𝛼𝑥𝑥−𝜆−𝐷𝑘+𝑒𝛼−𝑠𝛼−𝑠𝑘_{−Δ𝜏−𝜆𝜏}𝑦−𝜆−0𝑘𝑥𝑦−Δ+𝑘𝑒^{−Δ𝜏−𝜆𝜏}𝑥−𝜆⎤⎥⎥⎥⎦=0.(3.4) Direct calculation of the left-hand side of (3.4) gives

−Δ++𝑘𝑘𝑒𝑥𝑘𝑒^{−(Δ+𝜆)𝜏}𝑥−_{−(Δ+𝜆)𝜏}𝜆(−1−𝛼𝑥−𝜆)(−𝐷+𝛼𝑠−𝑘𝑦−𝜆)+𝛼²𝑠𝑥=(−Δ−𝜆)(1𝑦+(𝛼−𝑥1+−𝜆𝛼)𝑥(𝐷−𝜆−)𝛼𝑠+𝑘𝑦+𝜆)+𝛼₂×𝑘𝑦(−1−𝛼𝑥−𝜆)+(1+𝛼𝑥+𝜆)(𝐷−𝛼𝑠+𝑘𝑠𝑦𝑥++𝜆)𝑒+𝛼^{−(Δ+𝜆)𝜏}₂𝑘𝑥=(−Δ−𝜆)(1+𝛼𝑥+𝜆)(𝐷−𝛼𝑠+𝑘𝑦+𝜆)+𝛼₂𝑠𝑥+𝑒_{−(Δ+𝜆)𝜏}𝑠𝑥=+(−𝑒Δ−𝜆)_{−(Δ+𝜆)𝜏}𝑘{𝑥(𝜆+(11+)(𝛼𝜆𝑥++𝐷𝜆+)(𝐷𝑘𝑦−)+𝛼𝑠𝛼𝑥+(𝜆𝜆)++𝐷𝛼²+𝑠𝑥𝑘𝑦)−𝛼𝑠(𝜆+1)}𝑘𝑥{(𝜆+1)(𝜆+𝐷)+𝛼𝑥(𝜆+𝐷)−𝛼𝑠(𝜆+1)}.(3.5) For convenience, define 𝑃(𝜆)∶=+(𝑒−Δ−𝜆_{−(Δ+𝜆)𝜏}){(𝜆+1)(𝜆+𝐷+𝑘𝑦)+𝛼𝑥(𝜆+𝐷+𝑘𝑦)−𝛼𝑠(𝜆+1)}𝑘𝑥{(𝜆+1)(𝜆+𝐷)+𝛼𝑥(𝜆+𝐷)−𝛼𝑠(𝜆+1)}.(3.6)

Theorem 3.1. Equilibrium 𝐸₁ is stable if 𝛼<𝐷 and unstable if 𝛼>𝐷.

Proof. Evaluating the characteristic equation at 𝐸₁ gives 𝑃(𝜆)||_𝐸1=−(Δ+𝜆)(𝜆+1)(𝜆+𝐷−𝛼)=0.(3.7) The eigenvalues −1 and −Δ are both negative. The third eigenvalue is −𝐷+𝛼. Therefore the equilibrium 𝐸₁ is stable if 𝛼<𝐷 and unstable if 𝛼>𝐷.

Remark 3.2. If 𝛼<𝐷, then there is only one equilibrium, 𝐸₁. If 𝛼>𝐷, equilibrium 𝐸₂ also exists.

Lemma 3.3. Assume 𝛼>𝐷. The characteristic equation evaluated at 𝐸₂ has two negative eigenvalues, and the remaining eigenvalues are solutions of (𝜆+Δ)𝑒^{(𝜆+Δ)𝜏}=𝑘1𝐷−𝛼1.(3.8) In addition, the characteristic equation evaluated at 𝐸₂ has zero as an eigenvalue if and only if 𝜏=𝜏_𝑐.

Proof. Assume 𝛼>𝐷. Equilibrium 𝐸₂ exists. Consider the characteristic equation at 𝐸₂. Since (𝛼−𝐷)/𝛼𝐷=(1−𝑠)/𝛼𝑠 at 𝐸₂, 𝑃(𝜆)||_𝐸2=×{(𝜆−+𝜆1−)(Δ𝜆++𝑒𝐷)+𝛼_{−(Δ+𝜆)𝜏}𝑥(𝜆+𝐷)−𝛼𝑠(𝜆+1)}=(𝜆+1)(𝜆+𝐷)+1𝑘−𝑥𝑠×−𝜆−Δ+𝑒_{−(Δ+𝜆)𝜏}𝑠(𝜆+𝐷)−𝐷(𝜆+1)𝑘𝛼−𝐷=𝜆(𝜆+1)−(𝜆+𝐷)+𝜆𝛼+𝐷𝐷=−𝜆₂𝑠−𝜆−Δ+𝑘𝛼𝛼−𝐷𝐷𝑒^{−(Δ+𝜆)𝜏}+𝛼=−𝑒𝐷_{−(Δ+𝜆)𝜏}𝜆+𝛼−𝐷𝜆+Δ−𝑘𝛼𝛼−𝐷𝐷𝑒^{−(Δ+𝜆)𝜏}𝜆−𝜆₁𝜆−𝜆₂(𝜆+Δ)𝑒^{(Δ+𝜆)𝜏}−𝑘𝐷1−𝛼1=0,(3.9) where 𝜆₁+𝜆₂=−𝛼/𝐷 and 𝜆₁𝜆₂=𝛼−𝐷>0. Therefore, 𝜆₁ and 𝜆₂ 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
Δ𝑒^Δ𝜏=𝑘1𝐷−𝛼1.(3.10) Solving for 𝜏 gives 𝜏=1ΔlnΔ𝑘𝐷1−𝛼1=𝜏^𝑐.(3.11)

Theorem 3.4. Assume that 𝐷⩾1, Δ⩾1, 𝑘>0, 𝛼>0, and (𝑘/Δ)(1/𝐷−1/𝛼)⩾1 so that 𝜏_𝑐⩾0. Equilibrium 𝐸₂ is locally asymptotically stable if 𝜏>𝜏_𝑐 and unstable if 𝜏<𝜏_𝑐. If 𝐷=1, then equilibrium 𝐸₂ is globally asymptotically stable for 𝜏>(1/Δ)ln(𝑘/Δ).

Proof. Assume that 𝜏>𝜏_𝑐. Assumptions 𝑘>0, Δ⩾1, and (𝑘/Δ)(1/𝐷−1/𝛼)⩾1 imply 1/𝐷>1/𝛼, or equivalently 𝛼>𝐷. By Lemma 3.3, to prove that equilibrium 𝐸₂ 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/𝐷>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 𝜏=𝜏_𝑐, then 𝜆=0 is a solution of (3.8). Thus for 𝜏>𝜏_𝑐, any real root 𝜆 of (3.8) must satisfy −Δ<𝜆<0.
For any 𝜏=𝜏̃<𝜏_𝑐, we have (𝜆+Δ)𝑒^{(𝜆+Δ)𝜏}|_𝜏=𝜏̃_,𝜆=0<𝑘(1/𝐷−1/𝛼) and lim_𝜆→+∞(𝜆+Δ)𝑒^(𝜆+Δ)𝜏̃=+∞. Therefore there exists at least one 𝜆=𝜆̃>0 such that (𝜏̃,𝜆̃) is a solution of (3.8). Equilibrium 𝐸₂ is unstable if 𝜏<𝜏_𝑐.
In what follows, we prove that if 𝜏>𝜏_𝑐 all complex eigenvalues of (3.8) have negative real parts. Suppose that 𝜆+Δ=𝛾+𝑖𝛽(𝛽>0) is a solution of (3.8). Using the Euler formula, we have
𝛾cos(𝛽𝜏)−𝛽sin(𝛽𝜏)+𝑖(𝛾sin(𝛽𝜏)+𝛽cos(𝛽𝜏))=𝑘1𝐷−𝛼1𝑒^−𝛾𝜏.(3.12) Equating the real parts and imaginary parts of the equation, we have 𝛾cos(𝛽𝜏)−𝛽sin(𝛽𝜏)=𝑘1𝐷−𝛼1𝑒^−𝛾𝜏𝛾sin(𝛽𝜏)+𝛽cos(𝛽𝜏)=0.(3.13) Squaring both equations, adding, and taking the square root on both sides give 𝛾²+𝛽²𝑒^𝛾𝜏=𝑘1𝐷−𝛼1.(3.14) The left-hand side of (3.14) is monotonically increasing in 𝛾, 𝛽, and 𝜏 provided that 𝛾>0. Since (3.14) has solution 𝛾=Δ, 𝛽=0 at 𝜏=𝜏_𝑐, any roots of (3.14) must satisfy 𝛾<Δ since 𝜏>𝜏_𝑐. Hence Re{𝜆}=𝛾−Δ<0. Therefore (3.8) has no complex eigenvalue with nonnegative real part and so 𝐸₂ is locally asymptotically stable for 𝜏>𝜏_𝑐.
Assume that 𝐷=1. Now we prove that 𝐸₂ is globally asymptotically stable when 𝜏>(1/Δ)ln(𝑘/Δ), or equivalently 𝑘𝑒^−Δ𝜏<Δ. In this case, choose 𝜖₀>0 small enough such that 𝑘𝑒^−Δ𝜏(1+𝜖₀)<Δ. By Theorem 2.1, for such 𝜖₀, there exists a 𝑇>0 so that 0<𝑥(𝑡)<1+𝜖₀ for 𝑡>𝑇. Hence, for 𝑡>𝑇+𝜏, 𝑘𝑒^−Δ𝜏𝑥(𝑡−𝜏)<Δ. In Example 5.1 of Kuang ([10, page 32]), choose 𝜌(𝑡)=𝜏, 𝑎(𝑡)=Δ, 𝑏(𝑡)=𝑘𝑒^−Δ𝜏𝑥(𝑡−𝜏), and 𝛼=Δ/2. We obtain (𝑘𝑒^−Δ𝜏(1+𝜖₀))²<Δ²=4(Δ−𝛼)𝛼. Therefore 𝑦(𝑡)→0 as 𝑡→∞. Let 𝑧(𝑡)=𝑠(𝑡)+𝑥(𝑡). Noting 𝐷=1, from (2.3), we have 𝑧(̇𝑡)=1−𝑧(𝑡)−𝑘𝑥(𝑡)𝑦(𝑡). Multiply by the integrating factor 𝑒^𝑡, (𝑧(𝑡)𝑒^𝑡)^=𝑒^𝑡(1−𝑘𝑥(𝑡)𝑦(𝑡)). Integrating both sides from 0 to 𝑡 gives
𝑧(𝑡)=𝑒^−𝑡𝑧(0)+𝑒^−𝑡𝑒^𝑡−1−𝑒^−𝑡^𝑡0=1+𝑒_−𝑡𝑒^𝑠𝑘𝑥(𝑠)𝑦(𝑠)d𝑠(𝑧(0)−1)−𝑒^−𝑡^𝑡0𝑒^𝑠𝑘𝑥(𝑠)𝑦(𝑠)d𝑠.(3.15) If lim_𝑡→∞∫^𝑡0𝑒^𝑠𝑘𝑥(𝑠)𝑦(𝑠)d𝑠<∞, then lim_𝑡→∞𝑒^−𝑡∫^𝑡0𝑒^𝑠𝑘𝑥(𝑠)𝑦(𝑠)d𝑠=0. Therefore lim_𝑡→∞𝑧(𝑡)=1. If lim_𝑡→∞∫^𝑡0𝑒^𝑠𝑘𝑥(𝑠)𝑦(𝑠)d𝑠=∞, by L'Hôspital's rule, lim_𝑡→+∞𝑒^−𝑡^𝑡0𝑒^𝑠𝑘𝑥(𝑠)𝑦(𝑠)d𝑠=lim_𝑡→+∞∫^𝑡0𝑒^𝑠𝑘𝑥(𝑠)𝑦(𝑠)d𝑠𝑒^𝑡=lim_𝑡→+∞𝑒^𝑡𝑘𝑥(𝑡)𝑦(𝑡)𝑒^𝑡=lim_𝑡→+∞𝑘𝑥(𝑡)𝑦(𝑡)=0,(3.16) since 𝑥(𝑡) is bounded and lim_𝑡→∞𝑦(𝑡)=0. It again follows that lim_t→∞𝑧(𝑡)=1. Hence lim_𝑡→∞𝑠(𝑡)+𝑥(𝑡)=1.(3.17)
We show that lim_𝑡→∞𝑠(𝑡)=1/𝛼 and lim_𝑡→∞𝑥(𝑡)=(𝛼−1)/𝛼. First assume that the limits exist, that is, lim_𝑡→∞𝑠(𝑡)=𝑠 and lim_𝑡→∞𝑥(𝑡)=𝑥. From (2.3), we know that 𝑠(̇𝑡) and 𝑥̇(𝑡) are uniformly continuous since 𝑠(𝑡), 𝑥(𝑡), and 𝑦(𝑡) are bounded. By Theorem A.3, it follows that lim_𝑡→∞𝑠(̇𝑡)=0 and lim_𝑡→∞𝑥̇(𝑡)=0. Note that lim_𝑡→∞𝑦(𝑡)=0. Letting 𝑡→∞ in (2.3) gives
1−𝑠−𝛼𝑥𝑠=0,𝑥−1+𝛼𝑠=0.(3.18) Either (𝑠,𝑥)=(1,0) or (𝑠,𝑥)=(1/𝛼,(𝛼−1)/𝛼). Assume that (𝑠,𝑥)=(1,0), that is, lim_𝑡→∞𝑠(𝑡)=1 and lim_𝑡→∞𝑥(𝑡)=0. Note that 𝛼>𝐷. There exists 𝜖>0 such that 𝛼−𝐷−(𝛼+𝑘)𝜖>0. For such 𝜖, there exists a sufficiently large 𝑡 so that 𝑠(𝑡)>1−𝜖 and 0<𝑦(𝑡)<𝜖. Recalling that 𝑥(𝑡)>0, by (2.3) 𝑥̇(𝑡)>𝑥(𝑡)(−𝐷+𝛼(1−𝜖)−𝑘𝜖)=𝑥(𝑡)(𝛼−𝐷−𝛼𝜖−𝑘𝜖)>0,(3.19) for all sufficiently large 𝑡. Therefore it is impossible for 𝑥(𝑡) to approach 0 from above giving a contradiction. Therefore, we must have (𝑠,𝑥)=(1/𝛼,(𝛼−1)/𝛼).
Now suppose that the limits do not exist. In particular if 𝑥(𝑡) does not converge, then let 𝑥=limsup_𝑡→∞𝑥(𝑡) and 𝑥=liminf_𝑡→∞𝑥(𝑡). By Lemma A.2 in the appendix, there exists {𝑡_𝑚}↑∞ and {𝑠_𝑚}↑∞ such that
lliimm_{𝑚→∞𝑚→∞}𝑥𝑡_𝑚=𝑥,lim_𝑚→∞𝑥̇𝑡_𝑚=0,𝑥𝑠_𝑚=𝑥lim_𝑚→∞𝑥̇𝑠_𝑚=0.(3.20) From (2.3), 𝑥𝑡_𝑚−𝐷+𝛼𝑠𝑡_𝑚+𝑘𝑦𝑡_𝑚=0.(3.21) Noting that 𝑥(𝑡_𝑚)>0, we have 𝑠(𝑡_𝑚)=(1−𝑘𝑦(𝑡_𝑚))/𝛼. Since lim_𝑡→∞𝑦(𝑡)=0, lim_𝑡→∞𝑠(𝑡_𝑚)=1/𝛼. By (3.17), lim_𝑡→∞𝑥(𝑡_𝑚)=lim_𝑡→∞(𝑥(𝑡_𝑚)+𝑠(𝑡_𝑚))−𝑠(𝑡_𝑚)=1−1/𝛼=(𝛼−1)/𝛼. Therefore 𝑥=(𝛼−1)/𝛼. Similarly we can show that 𝑥=(𝛼−1)/𝛼. This implies that lim_𝑡→∞𝑥(𝑡)=(𝛼−1)/𝛼, a contradiction.
Since 𝑠(𝑡)+𝑥(𝑡) converges and 𝑥(𝑡) converges, then 𝑠(𝑡) must also converge. Hence lim_𝑡→∞𝑠(𝑡)=1/𝛼 and lim_𝑡→∞𝑥(𝑡)=(𝛼−1)/𝛼. It follows that 𝐸₂ is globally asymptotically stable.

4. Hopf Bifurcations at 𝐸₊ Assuming 𝐷=Δ=1

Now consider the stability of 𝐸₊. The characteristic equation at 𝐸₊ is

𝑃(𝜆)||_𝐸+=+(−𝑒Δ−𝜆)_{−(Δ+𝜆)𝜏}1+𝛼𝑥₊(𝜏)+𝜆𝐷−𝛼𝑠₊(𝜏)+𝑘𝑦₊(𝜏)+𝜆+𝛼²𝑠₊(𝜏)𝑥₊(𝜏)=(−Δ−𝜆)𝑘𝑥1₊+(𝜏)𝛼𝑥(₊𝜆+1)(𝜆+𝐷)+𝛼𝑥₊(𝜏)(𝜆+𝐷)−𝛼𝑠₊(𝜏)(𝜆+1)+𝑒_{−(Δ+𝜆)𝜏}(𝜏)+𝜆𝜆+𝛼²𝑠₊(𝜏)𝑥₊(𝜏)=(−Δ−𝜆)𝑘𝑥₊(𝜏1)(𝜆+1)(𝜆+𝐷)+𝛼𝑥₊(𝜏)(𝜆+𝐷)−𝛼𝑠₊(𝜏)(𝜆+1)𝑠₊(𝜏)+𝜆𝜆+𝛼1−𝑠⁺+𝑒_−𝜆𝜏(𝜏)Δ(𝜆+1)(𝜆+𝐷)+1−𝑠⁺(𝜏)𝑠₊(𝜏)(𝜆+𝐷)−𝛼𝑠⁺=(−Δ−𝜆)𝜆₂(𝜏)(𝜆+1)+𝜆𝑠₊(𝜏)+𝛼1−𝑠⁺+Δ𝑒_−𝜆𝜏(𝜏)𝜆+𝑠₊1(𝜏)(𝜆+𝐷)−𝛼𝑠⁺(𝜏)(𝜆+1)=0.(4.1) By assumption Δ=𝐷=1, and so

𝑃(𝜆)||_𝐸+=−(𝜆+1)𝜆²+𝜆𝑠₊(𝜏)+𝛼1−𝑠⁺(𝜏)+𝑒^−𝜆𝜏−𝜆+𝛼𝑠⁺(𝜏)−1=−(𝜆+1)𝜆₂𝑠₊(𝜏)+𝑝(𝜏)𝜆+𝛽(𝜏)+𝑒^−𝜆𝜏(𝑞𝜆+𝑐(𝜏))=0,(4.2) where

𝑝(𝜏)=1𝑠₊(𝜏),𝛽(𝜏)=𝛼1−𝑠⁺(𝜏),𝑞=−1,𝑐(𝜏)=𝛼𝑠₊(𝜏)−1𝑠₊(𝜏).(4.3) The characteristic equation at 𝐸₊ has one eigenvalue equal to −1 and the others are given by solutions of the equation

𝜆²+𝑝(𝜏)𝜆+𝛽(𝜏)+𝑒^−𝜆𝜏(𝑞𝜆+𝑐(𝜏))=0.(4.4)

Lemma 4.1. Assuming 𝑘>0, 𝛼>0, and 𝑘(1−1/𝛼)⩾1 so that 𝜏_𝑐=ln(𝑘(1−1/𝛼))⩾0, then 𝐸₊ has no zero eigenvalue for 𝜏∈(0,𝜏_𝑐).

Proof. Assume that 𝜏∈(0,𝜏_𝑐). By the method of contradiction, suppose that there exists a zero root of (4.4). Therefore 𝛽(𝜏)+𝑐(𝜏)=𝛼−1𝑠₊(𝜏)=0.(4.5) Noting that 𝜏_𝑐>0 if and only if 𝑘(1−1/𝛼)>1, for any 0<𝜏<𝜏_𝑐, 𝛼−1𝑠₊(𝜏)=𝛼−1−𝛼𝑘𝑒^𝜏>𝛼−1−𝛼1−1𝛼=0,(4.6) a contradiction.

Lemma 4.2. Assume 𝑘>0, 𝛼>0, 𝑘(1−1/𝛼)>1. Equilibrium 𝐸₊ is asymptotically stable when 𝜏=0.

Proof. For 𝜏=0, (4.4) reduces to 𝜆²+𝑝(0)𝜆+𝛽(0)+(𝑞𝜆+𝑐(0))=𝜆²+1𝑠₊(0)−1𝜆+𝛼−𝑠₊1(0).(4.7) Both coefficients are positive, since 1𝛼−1𝑠₊(0)−1=𝛼𝑘>0,𝑠₊(0)=𝛼−1−𝛼𝑘=𝛼1−𝛼1−𝑘1>0,(4.8) and 𝑘(1−1/𝛼)>1 implies 1−1/𝛼>1/𝑘. 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. Taking 𝑛=2 and 𝑔(𝜆,𝜏)=𝑝(𝜏)𝜆+(𝑞𝜆+𝑐(𝜏))𝑒^−𝜆𝜏+𝛽(𝜏) in Kuang [10, Theorem 1.4, page 66] gives limsup_{Re𝜆>0,|𝜆|→∞}||𝜆⁻²𝑔(𝜆,𝜏)||=0<1.(4.9) 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 𝜆=𝑖𝜔 (𝜔>0) is a root of 𝑃(𝜆)|_𝐸+=0,
−𝜔²+𝑖𝑝(𝜏)𝜔+𝛽(𝜏)+𝑒^{−𝑖𝜔𝜏}(𝑖𝑞𝜔+𝑐(𝜏))=0.(4.10) Substituting 𝑒^𝑖𝜃=cos𝜃+𝑖sin𝜃 into (4.10) gives −𝜔²+𝛽(𝜏)+𝑞𝜔sin(𝜔𝜏)+𝑐(𝜏)cos(𝜔𝜏)+𝑖(𝑝(𝜏)𝜔+𝑞𝜔cos(𝜔𝜏)−𝑐(𝜏)sin(𝜔𝜏))=0.(4.11) Separating the real and imaginary parts, we obtain 𝑐(𝑐𝜏()𝜏c)ossin(𝜔(𝜔𝜏)𝜏+)−𝑞𝜔𝑞𝜔sicno(s𝜔(𝜏𝜔)𝜏=)=𝜔²𝑝−(𝜏𝛽)𝜔(𝜏.),(4.12) Solving for cos(𝜔𝜏) and sin(𝜔𝜏) gives sin(𝜔𝜏)=𝑐(𝜏)𝑝(𝜏)𝜔+𝑞𝜔𝜔²−𝛽(𝜏)𝑐(𝜏)²+𝑞²𝜔²cos(𝜔𝜏)=𝑐(𝜏)𝜔₂,−𝛽(𝜏)−𝑞𝑝(𝜏)𝜔²𝑐(𝜏)²+𝑞²𝜔².(4.13) Noting sin²(𝜔𝜏)+cos²(𝜔𝜏)=1, squaring both sides of equations (4.13), adding, and rearranging gives 𝜔⁴+𝑝²(𝜏)−𝑞²−2𝛽(𝜏)𝜔²+𝛽²(𝜏)−𝑐²(𝜏)=0.(4.14) Solving for 𝜔, we obtain two roots 𝜔₁(𝜏) and 𝜔₂(𝜏): 𝜔₁(𝜏)=1√2𝑞²−𝑝²(𝜏)+2𝛽(𝜏)+𝑞²−𝑝²(𝜏)+2𝛽(𝜏)²−4(𝛽²(𝜏)−𝑐²(𝜏))^1/2=1𝑠₊(𝜏)√21−𝑠⁺+(𝜏)2𝛼𝑠²⁺(𝜏)−𝑠₊(𝜏)−1𝑠²⁺(𝜏)−1²+4𝛼𝑠²⁺(𝜏)𝑠²⁺(𝜏)−11−𝑠₊(𝜏)+4𝑠²⁺(𝜏)𝛼𝑠²⁺(𝜏)−1²^1/2𝜔₂(𝜏)=1√2𝑞²−𝑝²(𝜏)+2𝛽(𝜏)−𝑞²−𝑝²(𝜏)+2𝛽(𝜏)²−4𝛽²(𝜏)−𝑐²(𝜏)^1/2=1𝑠₊(𝜏)√21−𝑠⁺−(𝜏)2𝛼𝑠²⁺(𝜏)−𝑠₊(𝜏)−1𝑠²⁺(𝜏)−1²+4𝛼𝑠²⁺(𝜏)𝑠²⁺(𝜏)−11−𝑠₊(𝜏)+4𝑠²⁺(𝜏)𝛼𝑠²⁺(𝜏)−1²^1/2.(4.15)
Define conditions (𝐻₁) and (𝐻₂) as follows:
𝐻₁𝑞²𝑞−𝑝²(𝜏)+2𝛽(𝜏)>0,𝛽²(𝜏)−𝑐²(𝜏)>0,₂−𝑝²(𝜏)+2𝛽(𝜏)²−4𝛽²(𝜏)−𝑐²(𝜏)⩾0,(4.16)𝐻₂𝛽²(𝜏)−𝑐²(𝜏)<0,or𝛽²(𝜏)−𝑐²(𝜏)=0and𝑞²−𝑝²(𝜏)+2𝛽(𝜏)>0.(4.17)

Lemma 4.4. If (𝐻₁) holds for all 𝜏 in some interval 𝐼, then (4.14) has two positive roots 𝜔₁(𝜏)⩾𝜔₂(𝜏) for all 𝜏∈𝐼 with 𝜔₁(𝜏)>𝜔₂(𝜏) when all the inequalities in (𝐻₁) are strict. If (𝐻₂) holds for all 𝜏 in some interval 𝐼, then (4.14) has only one positive root, 𝜔₁(𝜏) for all 𝜏∈𝐼. If no interval exists where either (𝐻₁) or (𝐻₂) holds, then there are no positive real roots of (4.14).

Define the interval

𝐽=⎡⎢⎢⎣ln𝑘𝛼1/4+√1/116+1/(2𝛼)−1,ln⎛⎜⎜⎝𝑘⁴√𝛼−1𝛼⎞⎟⎟⎠⎤⎥⎥⎦.(4.18) When the end points of 𝐽 are real and 𝐽≠∅, define

𝐼₁=0,𝜏_𝑐∩𝐽.(4.19) We prove that (𝐻₁) holds for any 𝜏∈𝐼₁.

From 𝐷=Δ=1,

𝜏_𝑐=1ΔlnΔ𝑘𝐷1−𝛼1=ln𝑘(𝛼𝛼−1).(4.20) If 𝛼>1, then 𝛼>√⁴𝛼. It follows that

𝜏_𝑐>ln⎛⎜⎜⎝𝑘⁴√𝛼−1𝛼⎞⎟⎟⎠.(4.21) Therefore,

𝐼₁=⎡⎢⎢⎣max0,ln𝑘𝛼1/4+√1/116+1/(2𝛼)−1,ln⎛⎜⎜⎝𝑘⁴√𝛼−1𝛼⎞⎟⎟⎠⎤⎥⎥⎦.(4.22)

Theorem 4.5. Assume 𝛼>(7+3√5)/2 and 𝑘>𝛼/(√⁴𝛼−1), then 𝐼₁ is not empty, and for any 𝜏∈𝐼₁, but 𝜏≠ln((𝑘/𝛼)(1/(1/4+√1/16+1/(2𝛼))−1)), condition (𝐻₁) holds and 𝜔₁(𝜏)>𝜔₂(𝜏)>0. If 𝜏=ln((𝑘/𝛼)(1/(1/4+√1/16+1/(2𝛼))−1))∈𝐼₁, then 𝜔₁(𝜏)>𝜔₂(𝜏)=0.

Proof. For any 𝛼>(7+3√5)/2, we have 1−1/√⁴𝛼>0, and therefore ₄1√𝛼+12−√52<1₄7+3√5/2+12−√52=0.(4.23) Hence, √1⁴𝛼−14²−116+21𝛼²=1=−1√𝛼−2√⁴1𝛼−21𝛼2√⁴𝛼⎛⎜⎜⎝√1⁴𝛼³+1−₄2=1√𝛼⎞⎟⎟⎠2√⁴𝛼1−√⁴1𝛼⎛⎜⎜⎝√1⁴𝛼²+1₄=1√𝛼−1⎞⎟⎟⎠2√⁴𝛼1−√⁴1𝛼⎛⎜⎜⎝√1⁴𝛼+12²−5=14⎞⎟⎟⎠2√⁴𝛼1−√⁴1𝛼√⁴1𝛼+12−√52√⁴1𝛼+12+√52<0.(4.24) Therefore, 1/√⁴𝛼−1/4<√1/16+1/(2𝛼). Since 1/4+√1/16+1/(2𝛼)<1/4+1/16+1/(7+3√5)<1, it follows that ₄1√𝛼<41+116+21𝛼<1.(4.25) Hence, ln𝑘𝛼1/4+√1/116+1/(2𝛼)−1<ln⎛⎜⎜⎝𝑘⁴√𝛼−1𝛼⎞⎟⎟⎠.(4.26) From 𝑘>𝛼/(√⁴𝛼−1), we have ln(𝑘(√⁴𝛼−1)/𝛼)>0. Therefore, max0,ln𝑘𝛼1/4+√1/116+1/(2𝛼)−1<ln⎛⎜⎜⎝𝑘⁴√𝛼−1𝛼⎞⎟⎟⎠,(4.27) and so 𝐼₁ is not empty. Noting 𝑠₊(𝜏)=1/(1+(𝛼Δ/𝑘)𝑒^Δ𝜏) and Δ=1, for any 𝜏∈𝐼₁, but 𝜏≠ln((𝑘/𝛼)(1/(1/4+√1/16+1/(2𝛼))−1)), we have 𝑠₊(𝜏)∈[1/√⁴𝛼,1/4+√1/16+1/(2𝛼)).
In what follows, we intend to show that for any such 𝜏, condition (𝐻₁) holds. From (4.3),
𝑞²−𝑝²(𝜏)+2𝛽(𝜏)=(−1)²−1𝑠²⁺(𝜏)+2𝛼1−𝑠⁺=1−𝑠₊(𝜏)(𝜏)2𝛼𝑠²⁺(𝜏)𝑠⁺(𝜏)²−𝑠⁺(𝜏)=1−𝑠₊2𝛼−21𝛼(𝜏)2𝛼𝑠²⁺(𝜏)𝑠⁺(𝜏)−14𝛼²−116𝛼²−12𝛼.(4.28) Since 𝑠₊(𝜏)<1, to show that the first inequality in (𝐻₁) holds, it suffices to show that the factor on the right-hand side of the above expression is positive. Since 𝛼>(7+3√5)/2, 1/√𝛼−1/(4𝛼)=(1/√𝛼)(1−1/(4√𝛼))>0, and 116𝛼²+12𝛼²−1√𝛼−41𝛼²=116𝛼²+12𝛼−𝛼1+2𝛼1√𝛼−161𝛼²=12𝛼√1𝛼−1<0.(4.29) Since 1/√𝛼<1/√⁴𝛼 for 𝛼>(7+3√5)/2, 14𝛼+116𝛼²+12𝛼<√1𝛼<√1⁴𝛼.(4.30) For any 𝑠₊(𝜏)>1/√⁴𝛼, 𝑠₊(𝜏)−14𝛼⩾√⁴1𝛼−41𝛼>116𝛼²+12𝛼.(4.31) Hence, 𝑠⁺(𝜏)−41𝛼²⩾161𝛼²+21𝛼.(4.32) Next consider the second inequality in (𝐻₁). For 𝛼>(7+3√5)/2, since 1/√⁴𝛼>1/𝛼, 𝑠₊(𝜏)⩾1/√⁴𝛼>1/𝛼. Therefore, 𝛼𝑠₊(𝜏)>1. For 𝑠₊(𝜏)∈[1/√⁴𝛼,1/4+√1/16+1/(2𝛼))𝛽²(𝜏)−𝑐²(𝜏)==(𝛽𝛼(𝜏−)2−𝛼𝑐𝑠₊(𝜏))(𝛽(𝜏)+𝑐(𝜏))(𝜏)+1=−2𝛼𝑠₊(𝜏)𝛼−𝑠₊1(𝜏)𝑠²⁺(𝜏)𝑠²⁺(𝜏)−𝑠⁺(𝜏)2−21𝛼𝛼𝑠⁺=−2𝛼(𝜏)−1𝑠²⁺(𝜏)𝑠⁺(𝜏)−14²−116−21𝛼𝛼𝑠⁺(𝜏)−1>0.(4.33) Finally, 𝑞²−=𝑝²(1𝜏)−+2𝛽1(𝜏)²−4𝛽²(𝜏)−𝑐²(𝜏)=𝑞²−𝑝²(𝜏)𝑞²−𝑝²(𝜏)+4𝛽(𝜏)+4𝑐²(𝜏)𝑠²⁺(𝜏)1−𝑠²⁺1(𝜏)+4𝛼1−𝑠⁺(𝜏)+4𝛼𝑠₊(𝜏)−1𝑠₊(𝜏)²=1−1𝑠²⁺(𝜏)²+4𝛼1−1𝑠²⁺(𝜏)1−𝑠⁺(𝜏)+4𝑠₊(𝜏)²𝛼²−8𝛼+4=4𝑠₊𝑠²⁺(𝜏)(𝜏)²𝛼²+4𝛼1−1𝑠²⁺(𝜏)1−𝑠⁺(𝜏)−8𝛼+1−1𝑠²⁺(𝜏)²+4=4𝑠₊𝑠²⁺(𝜏)(𝜏)²𝛼²+4𝛼1−1𝑠²⁺(𝜏)1−𝑠⁺(𝜏)−2+1+1𝑠²⁺(𝜏)²=𝛼−𝛼₁𝛼−𝛼₂,(4.34) where 𝛼₁=2−1−1/𝑠²⁺(𝜏)1−𝑠₊(𝜏)+1/𝑠²⁺(𝜏)+2𝑠₊(𝜏)+11/𝑠₊(𝜏)−1²𝛼2𝑠²⁺(𝜏)₂,=2−1−1/𝑠²⁺(𝜏)1−𝑠₊(𝜏)−1/𝑠²⁺(𝜏)+2𝑠₊(𝜏)+11/𝑠₊(𝜏)−1²2𝑠²⁺(𝜏).(4.35) Since 𝑠₊(𝜏)<1, 2−1−1𝑠²⁺(𝜏)1−𝑠⁺(𝜏)=𝑠₊(𝜏)+1+1−𝑠⁺(𝜏)2−1−1𝑠²⁺(𝜏)>0,𝑠²⁺(𝜏)1−𝑠⁺(𝜏)²>2−1−1𝑠²⁺(𝜏)1−𝑠⁺(𝜏)²−𝑠²⁺(𝜏)1+1𝑠²⁺(𝜏)²=1𝑠²⁺(𝜏)+2𝑠⁺(𝜏)+11𝑠₊(𝜏)−1²>0.(4.36) It follows that 0<𝛼₂<𝛼₁. Again noting that 𝑠₊(𝜏)<1, 𝛼₁<2−1−1/𝑠²⁺(𝜏)1−𝑠₊(𝜏)+1/𝑠²⁺(𝜏)+2/𝑠₊(𝜏)+11/𝑠₊(𝜏)−1²=2−1−1/𝑠₂₊2𝑠²⁺(𝜏)(𝜏)1−𝑠₊(𝜏)+1/𝑠₊(𝜏)+11/𝑠₊(𝜏)−1=2−1−𝑠₊2𝑠²⁺(𝜏)(𝜏)−1/𝑠²⁺(𝜏)+1/𝑠₊(𝜏)+1/𝑠²⁺(𝜏)−1=𝑠2𝑠²⁺(𝜏)₊(𝜏)+2/𝑠²⁺(𝜏)−1/𝑠₊(𝜏)<12𝑠²⁺(𝜏)=12𝑠⁴⁺2(𝜏)−𝑠³⁺1(𝜏)−𝑠₊1(𝜏)2𝑠⁴⁺2(𝜏)=𝑠⁴⁺1(𝜏).(4.37) Hence, for any 𝑠₊(𝜏)>1/√⁴𝛼, we have 𝛼>1/𝑠⁴⁺(𝜏)>𝛼₁>𝛼₂. This leads to 𝑞²−𝑝²(𝜏)+2𝛽(𝜏)²−4𝛽²(𝜏)−𝑐²(𝜏)=𝛼−𝛼₁𝛼−𝛼₂>0.(4.38) Therefore, (𝐻₁) holds for any 𝜏∈𝐼₁. By Lemma 4.4, both 𝜔₁(𝜏)>0 and 𝜔₂(𝜏)>0.
If 𝜏=ln((𝑘/𝛼)(1/(1/4+√1/16+1/(2𝛼))−1))∈𝐼₁, we have 𝑠₊(𝜏)=1/4+√1/16+1/(2𝛼). Noting (4.25), we obtain
𝑞𝑞²−𝑝²(𝜏)+2𝛽(𝜏)>0,𝛽²(𝜏)−𝑐²(𝜏)=0,₂−𝑝²(𝜏)+2𝛽(𝜏)²−4𝛽²(𝜏)−𝑐²(𝜏)>0.(4.39) By (4.15), it follows that 𝜔₁(𝜏)>0 and 𝜔₂(𝜏)=0.

Now we define interval 𝐼₂ and prove that (𝐻₂) holds on 𝐼₂

𝐼₂∶=0,𝜏_𝑐∩−∞,ln𝑘𝛼1/4+√1/116+1/(2𝛼)−1.(4.40) In the following theorem, we consider the case that parameters are chosen so that

𝐼₂=0,ln𝑘𝛼1/4+√1/116+1/(2𝛼)−1.(4.41)

Theorem 4.6. Assume 𝛼>1 and 𝑘>𝛼(1/(1/4+√1/16+1/(2𝛼))−1)⁻¹. Interval 𝐼₂ given by (4.41) is not empty. For any 𝜏∈𝐼₂, (𝐻₂) holds and hence 𝜔₁(𝜏)>0.

Proof. Assume 𝛼>1. Letting 𝐺(𝛼)=14+116+21𝛼−𝛼1,(4.42) then dd𝛼𝐺(𝛼)=−4𝛼²√1/161+1/(2𝛼)+𝛼1²=1𝛼²√1−+18/𝛼+1>0.(4.43)𝐺(𝛼) is an increasing function of 𝛼 and 𝐺(1)=0. 𝐺(𝛼)>𝐺(1) implies that 1/4+√1/16+1/(2𝛼)−1/𝛼>0. Therefore 1>14+116+21𝛼>𝛼1.(4.44) This gives 𝛼−1>11/4+√1/16+1/(2𝛼)−1>0.(4.45) By assumption 𝑘>𝛼(1/(1/4+√1/16+1/(2𝛼))−1)⁻¹, we obtain 𝑘(𝛼−1)𝛼>𝛼𝑘1/4+√1/116+1/(2𝛼)−1>1.(4.46) Noting that 𝐷=Δ=1 and recalling the definition of 𝜏_𝑐 given in (3.2), 𝜏_𝑐=ln𝑘(𝛼−1)𝛼>ln𝛼𝑘1/4+√1/116+1/(2𝛼)−1>0.(4.47) Therefore, 𝐼₂ given by (4.41) is not empty. For any 𝜏∈𝐼₂, noting 𝑠₊(𝜏)=1/(1+(𝛼Δ/𝑘)𝑒^Δ𝜏) and Δ=1, we have 𝑠₊(𝜏)∈[1/4+√1/16+1/(2𝛼),1/(1+𝛼)/𝑘)⊂[1/4+√1/16+1/(2𝛼),1).
In what follows, we intend to show for any 𝜏∈𝐼₂, or equivalently 𝑠₊(𝜏)∈[1/4+√1/16+1/(2𝛼),1), (𝐻₂) holds. For any 𝑠₊(𝜏)∈[1/4+√1/16+1/(2𝛼),1), by (4.44), it follows that 𝑠₊(𝜏)>1/𝛼 and so 𝛼𝑠₊(𝜏)>1. Hence,
𝛽²(𝜏)−𝑐²(𝜏)=−2𝛼𝑠²⁺(𝜏)𝑠⁺(𝜏)−14²−116−21𝛼𝛼𝑠⁺(𝜏)−1⩽0.(4.48) For any 𝑠₊(𝜏)∈[1/4+√1/16+1/(2𝛼),1), we have 𝑠₊(𝜏)>1/(4𝛼)+√1/(16𝛼²)+1/(2𝛼), since 1/4>1/(4𝛼) and 1/16>1/(16𝛼²) imply that 14+116+21𝛼>41𝛼+116𝛼²+12𝛼.(4.49) Therefore, 𝑞²−𝑝²(𝜏)+2𝛽(𝜏)=1−𝑠⁺(𝜏)2𝛼𝑠²⁺(𝜏)𝑠⁺(𝜏)−14𝛼²−116𝛼²−12𝛼>0.(4.50)
Condition (𝐻₂) holds. By Lemma 4.4, 𝜔₁(𝜏)>0.

Next, to determine whether (4.2) has a pair of pure imaginary eigenvalues, we consider

𝑐²(𝜏)+𝑞²𝜔²=𝛼𝑠₊(𝜏)−1𝑠₊(𝜏)²+(−1)²𝜔²=𝛼𝑠₊(𝜏)−1𝑠₊(𝜏)²𝑐(𝜏)(𝑝(𝜏)𝜔)+𝑞𝜔𝜔₂+𝜔²,=𝜔𝑐(𝜏)𝑝(𝜏)+𝑞−𝜔𝛽₂(𝜏)=𝜔𝛼−1−𝛽(𝜏)𝑠²⁺(𝜏)−𝜔²𝑐(𝜏)𝜔₂+𝛼1−𝑠₊(𝜏),=𝛼−𝑠₊𝛽(𝜏)+𝑞𝜔(−𝑝(𝜏)𝜔)(𝜏)−1𝑠₊(𝜏)𝜔²−𝛼1−𝑠₊(𝜏)+(−1)𝜔−𝜔=𝛼𝑠₊𝑠₊(𝜏)(𝜏)𝜔²−1−𝑠₊(𝜏)𝛼𝑠₊(𝜏)−1𝑠₊(𝜏).(4.51) We obtain

sin(𝜔𝜏)=−𝜔𝜔²−𝛼−1/𝑠²⁺(𝜏)−𝛼1−𝑠₊(𝜏)𝛼𝑠₊(𝜏)−1/𝑠₊(𝜏)²+𝜔²cos(𝜔𝜏)=𝛼𝑠₊,(𝜏)𝜔²−1−𝑠₊(𝜏)𝛼−1/𝑠²⁺(𝜏)𝛼𝑠₊(𝜏)−1/𝑠₊(𝜏)²+𝜔².(4.52)

If there exists (𝜏,𝜔) satisfying (4.52), then (4.2) has a pair of pure imaginary roots ±𝑖𝜔. A necessary condition for (4.52) to have solutions is 𝛼𝑠²⁺(𝜏)≠1. Otherwise, 𝛼𝑠²⁺(𝜏)=1, and the second equation of (4.52) becomes cos(𝜔𝜏)=𝛼𝑠₊(𝜏). However, for any 𝜏∈(0,𝜏_𝑐), we have 𝛼𝑠₊(𝜏)>1, since

𝛼𝑠₊(𝜏)=𝛼1+𝛼𝑒^𝜏/𝑘>1+(𝛼/𝑘)𝛼𝑘(1−1/𝛼)=1+𝛼(1𝛼−1/𝛼)=𝛼𝛼=1.(4.53) Hence, the second equation of (4.52) has no solution. Assume 𝛼𝑠²⁺(𝜏)≠1 for 𝜔⩾0 and 𝜏∈[0,𝜏_𝑐] and denote the right-hand sides of (4.52) by

ℎ₁(𝜔,𝜏)=𝑐(𝜏)(𝑝(𝜏)𝜔)+𝑞𝜔𝜔²−𝛽(𝜏)𝑐²(𝜏)+𝑞²𝜔²=−𝜔𝜔²−𝛼−1/𝑠²⁺(𝜏)−𝛼1−𝑠₊(𝜏)𝛼𝑠₊(𝜏)−1/𝑠₊(𝜏)²+𝜔²ℎ₂,(𝜔,𝜏)=𝑐(𝜏)𝜔²−𝛽(𝜏)+𝑞𝜔(−𝑝(𝜏)𝜔)𝑐²(𝜏)+𝑞²𝜔²=𝛼𝑠₊(𝜏)𝜔²−1−𝑠₊(𝜏)𝛼−1/𝑠²⁺(𝜏)𝛼𝑠₊(𝜏)−1/𝑠₊(𝜏)²+𝜔².(4.54) Define functions

𝜃𝜃₁(𝜏)=arccosℎ₂𝜔₁(𝜏),𝜏ifℎ₂𝜔₁(𝜏),𝜏∈[−1,1₂],(𝜏)=arccosℎ₂𝜔₂(𝜏),𝜏ifℎ₂𝜔₂(𝜏),𝜏∈[−1,1].(4.55)

Lemma 4.7. Assume 𝛼>(7+3√5)/2 and 𝑘>𝛼/(√⁴𝛼−1). For any 𝜏∈𝐼₁ given by (4.19), there exists 𝜖_𝑗>0 and 𝜃_𝑗(𝜏) with 𝜖_𝑗⩽𝜃_𝑗(𝜏)⩽𝜋 (j=1,2) such that csoins𝜃𝜃_𝑗𝑗(𝜏)+2𝑛𝜋=ℎ₁𝜔_𝑗(𝜏),𝜏,(𝜏)+2𝑛𝜋=ℎ₂𝜔_𝑗(𝜏),𝜏.𝑛=0,1,2,…,(4.56)

Proof. For any 𝜏∈𝐼₁, by Theorem 4.5, 𝜔₁(𝜏)>0 and 𝜔₂(𝜏)⩾0. It is easy to see that ℎ₁(0,𝜏)=0 and lim_𝜔→+∞ℎ₁(𝜔,𝜏)=−∞. There are two roots of ℎ₁(𝑧,𝜏)=0, 𝑧₁=0 and 𝑧₂(𝜏)=𝛼−1𝑠²⁺(𝜏)+𝛼1−𝑠⁺(𝜏).(4.57) Hence, ℎ₁(𝜔,𝜏)>0 for 0<𝜔<𝑧₂(𝜏). For any 𝜏∈𝐼₁, as shown in Theorem 4.5, 𝑠₊(𝜏)∈[1/√⁴𝛼,1/4+√1/16+1/(2𝛼)]. This implies 𝑠₊(𝜏)⩾1/√⁴𝛼>1/√𝛼>1/𝛼. Therefore, 𝛼𝑠²⁺(𝜏)>1 and 𝛼𝑠₊(𝜏)>1. The function ℎ₂(𝜔,𝜏) is monotonically increasing for 𝜔⩾0, since 𝜕ℎ₂(𝜔,𝜏)𝜕𝜔=𝛼𝑠⁺(𝜏)2𝜔𝛼𝑠⁺(𝜏)−1/𝑠₊(𝜏)²+𝜔²−2𝜔𝜔²−1−𝑠₊(𝜏)𝛼−1/𝑠²⁺(𝜏)𝛼𝑠⁺(𝜏)−1/𝑠₊(𝜏)²+𝜔²²=2𝛼𝑠₊(𝜏)𝜔𝛼𝑠⁺(𝜏)−1/𝑠₊(𝜏)²+1−𝑠₊(𝜏)𝛼−1/𝑠²⁺(𝜏)𝛼𝑠⁺(𝜏)−1/𝑠₊(𝜏)²+𝜔²²=2𝛼𝑠₊(𝜏)𝜔𝛼−1/𝑠²⁺(𝜏)𝑠²⁺(𝜏)𝛼−1/𝑠²⁺(𝜏)+1−𝑠₊(𝜏)𝛼𝑠⁺(𝜏)−1/𝑠₊(𝜏)²+𝜔²²=2𝛼𝑠₊(𝜏)𝜔𝛼−1/𝑠²⁺(𝜏)𝑠₊(𝜏)𝛼𝑠₊(𝜏)−1𝛼𝑠⁺(𝜏)−1/𝑠₊(𝜏)²+𝜔²²=2𝛼𝜔𝛼𝑠²⁺(𝜏)−1𝛼𝑠₊(𝜏)−1𝛼𝑠⁺(𝜏)−1/𝑠₊(𝜏)²+𝜔²²⩾0.(4.58) Since 𝑠₊(𝜏)<1, ℎ₂(0,𝜏)=−𝛼𝑠₊(𝜏)1−𝑠⁺(𝜏)𝛼−1/𝑠²⁺(𝜏)𝛼𝑠₊(𝜏)−1/𝑠₊(𝜏)²=−𝛼𝑠⁺(𝜏)1−𝑠₊(𝜏)𝛼𝑠²⁺(𝜏)−1<0.(4.59) Also, lim_𝜔→∞ℎ₂(𝜔,𝜏)=𝛼𝑠₊(𝜏)>1. Therefore, there exists a unique 𝜔=𝑙_max(𝜏)=𝛼−1/𝑠²⁺(𝜏)>0, such that ℎ₂(𝑙_max(𝜏),𝜏)=1. Solving ℎ₂(𝑙_max,𝜏)=1 for 𝑙_max and noting that 𝛼−1/𝑠²⁺(𝜏)≠0, it is easy to see that 𝑙_max(𝜏)=𝛼−1/𝑠²⁺(𝜏)𝑙^2max𝛼𝑠₊(𝜏)−1=𝛼𝑠₊(𝜏)1−𝑠₊(𝜏)𝛼−1𝑠²⁺(𝜏)+𝛼𝑠⁺(𝜏)−1𝑠₊(𝜏)²=𝛼𝑠₊(𝜏)𝛼−1𝑠²⁺(𝜏)−𝛼𝑠⁺(𝜏)+1𝑠₊(𝜏)+𝛼²𝑠²⁺(𝜏)−2𝛼+1=𝛼₂𝑠²⁺(𝜏)𝑠₊(𝜏)−𝛼𝑠₊(𝜏)−𝛼²𝑠²⁺(𝜏)+𝛼+𝛼²𝑠²⁺(𝜏)−2𝛼+1=𝛼₂𝑠²⁺(𝜏)𝑠₊(𝜏)−𝛼=𝛼−1𝑠₊(𝜏)−𝛼+𝑠²⁺1(𝜏)𝑠²⁺(𝜏)𝛼𝑠⁺(𝜏)−1.(4.60) Then, ℎ₂(𝜔,𝜏)⩽1 for any 𝜔∈[0,𝑙_max(𝜏)]. Since 𝑠₊(𝜏)<1, 𝑙_max(𝜏)<𝑧₂(𝜏). Therefore, ℎ₁(𝜔,𝜏)>0 for any 𝜔∈[0,𝑙_max(𝜏)]. Since 𝜔₁(𝜏) is a positive root of ℎ²¹(𝜔,𝜏)+ℎ²²(𝜔,𝜏)=1, we have ℎ₂(𝜔₁(𝜏),𝜏)⩽1, which implies that 0<𝜔₁(𝜏)⩽𝑙_max(𝜏)<𝑧₂(𝜏). Therefore, ℎ₁(𝜔₁(𝜏),𝜏)>0, and so ℎ₁(𝜔₁(𝜏),𝜏)=√1−ℎ₂(𝜔₁(𝜏),𝜏). In fact, 𝜔₁(𝜏)<𝑙_max(𝜏),(4.61) since ℎ²²𝑙_max(𝜏),𝜏+ℎ²¹𝑙_max(𝜏),𝜏=1+ℎ²¹𝑙_max(𝜏),𝜏>1.(4.62) Thus, 𝜃₁(𝜏) is defined and 0⩽𝜃₁(𝜏)⩽𝜋. Since cos(𝜃₁(𝜏)+2𝑛𝜋)=ℎ₂(𝜔₁(𝜏),𝜏), sin𝜃₁(𝜏)+2𝑛𝜋==1−cos²𝜃₁(𝜏)+2𝑛𝜋=ℎ₁1−ℎ²²𝜔₁(𝜏),𝜏𝜔₁(𝜏),𝜏.(4.63) Hence, 𝜃₁(𝜏) satisfies (4.56). From (4.61), ℎ₂(𝜔₁(𝜏),𝜏)<ℎ₂(𝑙_max(𝜏),𝜏)=1, and so 𝜃₁(𝜏)>0. Since 𝜃₁(𝜏) is continuous on the interval 𝐼₁ and 𝐼₁ is closed, there exists 𝜖₁>0 such that 𝜃₁(𝜏)⩾𝜖₁. Similarly we can prove the existence of 𝜃₂(𝜏).

Lemma 4.8. Assume 𝛼>1 and 𝑘>𝛼(1/(1/4+√1/16+1/2𝛼)−1)⁻¹. For any 𝜏∈𝐼₂ given by (4.41), there exists 𝜖>0 and 𝜃₁(𝜏) such that 𝜖⩽𝜃₁(𝜏)<𝜋 and 𝜃₁(𝜏) satisfies (4.56) for 𝑗=1.

Proof. For any 𝜏∈𝐼₂, by Theorem 4.6, only 𝜔₁(𝜏)>0.
As in Lemma 4.7, we have ℎ₁(𝜔,𝜏)>0 for 0<𝜔<𝑧₂(𝜏). For any 𝜏∈𝐼₂, as shown in Theorem 4.6, 𝑠₊(𝜏)∈[1/4+√1/16+1/(2𝛼),1). Letting
𝐺(𝛼)=14+116+21𝛼−√1𝛼,𝛼>1,(4.64) we have dd𝛼𝐺(𝛼)=−4𝛼²√1/161+1/(2𝛼)+2𝛼1√𝛼=21𝛼√𝛼²−/41+2𝛼+√1𝛼>0.(4.65)𝐺(𝛼) is an increasing function and 𝐺(1)=0. Since 𝐺(𝛼)>𝐺(1), it follows that 1/4+√1/16+1/(2𝛼)>1/√𝛼. Therefore, for any 𝑠₊(𝜏)>1/4+√1/16+1/(2𝛼), we obtain 𝑠₊(𝜏)>1/√𝛼, or equivalently 𝛼>1/𝑠²⁺(𝜏). Since 𝜕ℎ₂(𝜔,𝜏)𝜕𝜔=2𝛼𝑠⁺(𝜏)𝜔𝛼𝑠⁺(𝜏)−1/𝑠₊(𝜏)²+1−𝑠₊(𝜏)𝛼−1/𝑠²⁺(𝜏)𝛼𝑠⁺(𝜏)−1/𝑠₊(𝜏)²+𝜔²²⩾0,(4.66)ℎ₂(𝜔,𝜏) is monotonically increasing for any 𝜔⩾0. For any 𝑠₊(𝜏)>1/4+√1/16+1/(2𝛼), 𝑠⁺(𝜏)−14²−116−21𝛼=𝑠²⁺(𝜏)−𝑠⁺(𝜏)2−21𝛼=21𝛼2𝛼𝑠²⁺=1(𝜏)−𝛼𝑠₊(𝜏)−12𝛼𝛼𝑠²⁺=𝛼𝑠₂₊(𝜏)−1−𝛼𝑠₊(𝜏)−𝛼𝑠²⁺(𝜏)(𝜏)−12𝛼1−𝛼𝑠⁺(𝜏)1−𝑠₊(𝜏)𝛼𝑠²⁺(𝜏)−1>0,(4.67) which implies that (𝛼𝑠₊(𝜏)(1−𝑠₊(𝜏)))/(𝛼𝑠²⁺(𝜏)−1)<1. Hence, 0>ℎ₂(0,𝜏)=−𝛼𝑠₊(𝜏)1−𝑠⁺(𝜏)𝛼−1/𝑠²⁺(𝜏)𝛼𝑠₊(𝜏)−1/𝑠₊(𝜏)²=−𝛼𝑠⁺(𝜏)1−𝑠₊(𝜏)𝛼𝑠²⁺(𝜏)−1>−1.(4.68) For 𝜏∈[0,𝜏_𝑐), lim_𝜔→∞ℎ₂(𝜔,𝜏)=𝛼𝑠₊(𝜏)>1, since 𝛼𝑠₊(𝜏)=𝛼1+𝛼𝑒^𝜏/𝑘>1+(𝛼/𝑘)𝛼𝑘(1−1/𝛼)=1+𝛼(1𝛼−1/𝛼)=𝛼𝛼=1.(4.69) As in the proof of Lemma 4.7, there exists a unique 𝑙_max(𝜏)=𝛼−1/𝑠²⁺(𝜏)>0 such that ℎ₂(𝑙_max(𝜏),𝜏)=1. Then 𝑙_max(𝜏)<𝑧₂(𝜏). Therefore, ℎ₁(𝜔,𝜏)>0 for any 𝜔∈[0,𝑙_max(𝜏)]. The rest of the proof is similar to the proof of Lemma 4.7. Furthermore, 𝜃₁(𝜏)<𝜋, since ℎ₂(𝜔₁(𝜏),𝜏)>−1 for any 𝜔₁(𝜏)∈[0,𝑙_max(𝜏)).

Theorem 4.9. Consider system (2.3) with 𝐷=Δ=1. (1)Suppose 𝛼>(7+3√5)/2, 𝑘>𝛼/(√⁴𝛼−1), and 𝜏∈𝐼₁ given by (4.22). For 𝜏∈𝐼₁ and 𝑗=1,2, 𝜔_𝑗(𝜏) is nonnegative and there exists 𝜖_𝑗>0 and 𝜃_𝑗(𝜏) such that 𝜖_𝑗⩽𝜃_𝑗(𝜏)⩽𝜋 and 𝜃_𝑗(𝜏) satisfies (4.56). If there exists 𝑛⩾0 such that 𝜃_𝑗(𝜏)+2𝑛𝜋 intersects 𝜏𝜔_𝑗(𝜏) at some 𝜏^𝑗𝑛∈𝐼₁, then (4.4) has a pair of pure imaginary eigenvalues 𝜆=±𝑖𝜔_𝑗(𝜏^𝑗𝑛). System (2.3) undergoes a Hopf bifurcation at 𝜏^𝑗𝑛 provided 𝑑𝑅𝑒(𝜆(𝜏))/𝑑𝜏|_𝜏=𝜏^𝑗𝑛≠0.(2)Suppose 𝛼>1, 𝑘>𝛼(1/(1/4+√1/16+1/(2𝛼))−1)⁻¹, and 𝜏∈𝐼₂ given by (4.41). For 𝜏∈𝐼₂, only 𝜔₁(𝜏) is positive and there exists 𝜖>0 and 𝜃₁(𝜏) such that 𝜖⩽𝜃₁(𝜏)<𝜋 and 𝜃₁(𝜏) satisfies (4.56) for 𝑗=1. If there exists 𝑛⩾0 such that 𝜃₁(𝜏)+2𝑛𝜋 intersects 𝜏𝜔₁(𝜏) at some 𝜏^1𝑛∈𝐼₂, then (4.4) has a pair of pure imaginary eigenvalues 𝜆=±𝑖𝜔₁(𝜏^1𝑛). System (2.3) undergoes a Hopf bifurcation at 𝜏^1𝑛 provided 𝑑𝑅𝑒(𝜆(𝜏))/𝑑𝜏|_𝜏=𝜏^1𝑛≠0.

Proof. Assume 𝐷=Δ=1 in system (2.3).Case 1. Suppose 𝜏∈𝐼₁. By Theorem 4.5, 𝜔_𝑗(𝜏)⩾0 for 𝑗=1,2. By Lemma 4.7, there exists 𝜖_𝑗>0 and 𝜃_𝑗(𝜏) such that 𝜖_𝑗⩽𝜃_𝑗(𝜏)⩽𝜋 and 𝜃_𝑗(𝜏) satisfies (4.56). Assume that there exists a positive integer 𝜏^𝑗𝑛∈𝐼₁ such that 𝜃_𝑗(𝜏^𝑗𝑛)+2𝑛𝜋=𝜏^𝑗𝑛𝜔_𝑗(𝜏^𝑗𝑛) for some integer 𝑛⩾0. Then system (4.52) has one solution (𝜏^𝑗𝑛,𝜔_𝑗(𝜏^𝑗𝑛)). Equation (4.4) has a pair of pure imaginary eigenvalues 𝜆=±𝑖𝜔_𝑗(𝜏^𝑗𝑛).
In what follows, we show that the conditions required for a Hopf Bifurcation (see Theorem A.1 in the appendix) are satisfied by the linearization (3.3) of (2.3) at 𝐸₊. In (A.1), choosing 𝜏 as the bifurcation parameter and letting
𝐷𝜏,𝑧_𝑡=⎡⎢⎢⎢⎣100010001⎤⎥⎥⎥⎦⎡⎢⎢⎢⎣𝑧¹𝑧(𝑡)₂𝑧(𝑡)₃𝐿𝜏,𝑧_𝑡(𝑡)⎤⎥⎥⎥⎦,=⎡⎢⎢⎢⎣−1−𝛼𝑥⁺𝛼𝑥₊−𝛼𝑠₊0−𝐷+𝛼𝑠₊−𝑘𝑦₊−𝑘𝑥₊00−Δ⎤⎥⎥⎥⎦⎡⎢⎢⎢⎣𝑧¹𝑧(𝑡)₂𝑧(𝑡)₃(𝑡)⎤⎥⎥⎥⎦+⎡⎢⎢⎢⎣000𝑘𝑒00^Δ𝜏𝑦₊𝑘𝑒00^Δ𝜏𝑥₊⎤⎥⎥⎥⎦⎡⎢⎢⎢⎣𝑧𝑧₁(𝑡−𝜏)₂𝑧(𝑡−𝜏)₃(𝑡−𝜏)⎤⎥⎥⎥⎦,(4.70) the linearization (3.3) of (2.3) at 𝐸₊ is of the form (A.1). Taking 𝑎 to be any positive real number and 𝑏=1/2, hypothesis (𝑆₁) in the Hopf Bifurcation Theorem holds, since |||||det_∞_𝑘=0𝐴_𝑘(𝛼)𝑒^{−𝜆𝑟𝑘(𝛼)}|||||det_∞|||||=|||||||det⎡⎢⎢⎢⎣100010001⎤⎥⎥⎥⎦|||||||=1⩾21,_𝑘=0𝐴_𝑘(𝛼)𝑒^{−𝜆𝑟𝑘(𝛼)}+⁰⁻¹𝐴(𝛼,𝜃)𝑒^𝜆𝜃d𝜃|||||=|||||||det⎡⎢⎢⎢⎣100010001⎤⎥⎥⎥⎦|||||||=1⩾12(4.71) for all 𝜏∈ℝ and |Re𝜆|<𝑎.
The characteristic equation (4.4) of (3.3) at 𝐸₊ has a pair of pure imaginary eigenvalues 𝜆=±𝜔₁(𝜏^𝑗𝑛) and no other root of (4.4) is an integral multiple of ±𝜔₁(𝜏^𝑗𝑛). Hence the hypothesis (𝑆₂) in the Hopf Bifurcation Theorem holds. Therefore, (2.3) undergoes a Hopf bifurcation at 𝐸₊ when 𝜏=𝜏^𝑗𝑛 provided dRe(𝜆(𝜏))/d𝜏|_𝜏=𝜏^𝑗𝑛≠0.Case 2. Suppose 𝜏∈𝐼₂. By Theorem 4.5, only 𝜔₁(𝜏)>0. By Lemma 4.8, there exists 𝜖>0 and 𝜃₁(𝜏) such that 𝜖⩽𝜃₁(𝜏)<𝜋 and 𝜃₁(𝜏) satisfies (4.56). Assume there exists 𝜏^1𝑛∈𝐼₂ such that 𝜃₁(𝜏^1𝑛)+2𝑛𝜋=𝜏^1𝑛𝜔₁(𝜏^1𝑛) for some integer 𝑛⩾0. Then system (4.52) has one solution (𝜏^1𝑛,𝜔₁(𝜏^1𝑛)). Equation (4.4) has a pair of pure imaginary eigenvalues 𝜆=±𝑖𝜔₁(𝜏^1𝑛). The rest of the proof is similar to that of Case 1 when 𝑗=1.

Corollary 4.10. Consider system (2.3) with 𝐷=Δ=1. (1)Suppose 𝛼>(7+3√5)/2, 𝑘>𝛼/(√⁴𝛼−1), and 𝜏∈𝐼₁ given by (4.22). For 𝜏∈𝐼₁, 𝜔_𝑗(𝜏) is nonnegative and there exists 𝜖_𝑗>0 and 𝜃_𝑗(𝜏) such that 𝜖_𝑗⩽𝜃_𝑗(𝜏)⩽𝜋 and 𝜃_𝑗(𝜏) satisfies (4.56) for 𝑗=1,2. If there exists a positive integer 𝑛_𝑗⩾0 such that min_𝜏∈𝐼1𝜏𝜔_𝑗(𝜏)⩽2𝑛_𝑗𝜋 and max_𝜏∈𝐼1𝜏𝜔_𝑗(𝜏)>(2𝑛_𝑗+1)𝜋, then 𝜃_𝑗(𝜏)+2𝑛_𝑗𝜋 intersects 𝜏𝜔_𝑗(𝜏) at least once at some 𝜏^𝑗𝑛_𝑗∈𝐼₁. System (2.3) undergoes a Hopf bifurcation at 𝜏^𝑗𝑛_𝑗 provided 𝑑𝑅𝑒(𝜆(𝜏))/𝑑𝜏|_𝜏=𝜏^𝑗𝑛𝑗≠0.(2)Suppose 𝛼>1, 𝑘>𝛼(1/(1/4+√1/16+1/(2𝛼))−1)⁻¹, and 𝜏∈𝐼₂ defined in (4.41). For 𝜏∈𝐼₂, only 𝜔₁(𝜏) is positive. There exists 𝜖>0 and 𝜃₁(𝜏) such that 𝜖⩽𝜃₁(𝜏)<𝜋 and 𝜃₁(𝜏) satisfies (4.56) for 𝑗=1. If there exists a positive integer 𝑁⩾0 such that max_𝜏∈𝐼2𝜏𝜔₁(𝜏)>(2𝑁+1)𝜋, then for any 0⩽𝑛⩽𝑁, 𝜃₁(𝜏)+2𝑛𝜋 intersects 𝜏𝜔₁(𝜏) at least once at some 𝜏^1𝑛∈𝐼₂. System (2.3) undergoes a Hopf bifurcation at 𝜏^1𝑛 provided 𝑑𝑅𝑒(𝜆(𝜏))/𝑑𝜏|_𝜏=𝜏^1𝑛≠0.

Proof. Assume 𝐷=Δ=1 in system (2.3).Case 1. Suppose 𝜏∈𝐼₁. By Theorem 4.5, 𝜔_𝑗(𝜏)⩾0 for 𝑗=1,2. By Lemma 4.7, there exists 𝜖_𝑗>0 and 𝜃_𝑗(𝜏) such that 𝜖_𝑗⩽𝜃_𝑗(𝜏)⩽𝜋 and 𝜃_𝑗(𝜏) satisfies (4.56). Assume that there exists a positive integer 𝑛_𝑗⩾0 such that min_𝜏∈𝐼1𝜏𝜔_𝑗(𝜏)⩽2𝑛_𝑗𝜋 and max_𝜏∈𝐼1𝜏𝜔_𝑗(𝜏)>(2𝑛_𝑗+1)𝜋. For such 𝑛_𝑗, min_𝜏∈𝐼1𝜏𝜔_𝑗(𝜏)<𝜖_𝑗+2𝑛_𝑗𝜋⩽𝜃_𝑗(𝜏)+2𝑛_𝑗𝜋⩽2𝑛_𝑗+1𝜋<max_𝜏∈𝐼1𝜏𝜔_𝑗(𝜏).(4.72) By the Mean Value Theorem, there exists 𝜏^𝑗𝑛_𝑗∈𝐼₁ such that 𝜃_𝑗(𝜏^𝑗𝑛_𝑗)+2𝑛_𝑗𝜋=𝜏^𝑗𝑛_𝑗𝜔_𝑗(𝜏^𝑗𝑛_𝑗). By Theorem 4.9. Case 1, the conclusion follows.Case 2. Suppose 𝜏∈𝐼₂. By Theorem 4.5, only 𝜔₁(𝜏)>0. By Lemma 4.8, there exists 𝜖>0 and 𝜃₁(𝜏) such that 𝜖⩽𝜃₁(𝜏)<𝜋 and 𝜃₁(𝜏) satisfies (4.56). Assume that there exists a positive integer 𝑁⩾0 such that max_𝜏∈𝐼2𝜏𝜔₁(𝜏)>(2𝑁+1)𝜋. By (4.41), 0∈𝐼₂. Therefore min_𝜏∈𝐼2𝜏𝜔₁(𝜏)=0. For 0⩽𝑛⩽𝑁, min_𝜏∈𝐼2𝜏𝜔₁(𝜏)<𝜖+2𝑛𝜋⩽𝜃_𝑖(𝜏)+2𝑛𝜋⩽(2𝑛+1)𝜋<max_𝜏∈𝐼2𝜏𝜔₁(𝜏).(4.73) By the Mean Value Theorem, there exists 𝜏^1𝑛∈𝐼₂ such that 𝜃₁(𝜏^1𝑛)+2𝑛𝜋=𝜏^1𝑛𝜔₁(𝜏^1𝑛). By Theorem 4.9. Case 2, the conclusion follows.

Corollary 4.11. Consider system (2.3) with 𝐷=Δ=1. Assume 𝛼>(7+3√5)/2 and 𝑘>𝛼/(√⁴𝛼−1). If 1/4+√1/16+1/(2𝛼)>𝑘/(𝛼+𝑘), then 𝐼₁=[0,ln((𝑘/𝛼)(√⁴𝛼−1))], where 𝐼₁ was defined in (4.22). For any 𝜏∈𝐼₁, 𝜔_𝑗(𝜏) is nonnegative and there exists 𝜖_𝑗>0 and 𝜃_𝑗(𝜏) such that 𝜖_𝑗⩽𝜃_𝑗(𝜏)⩽𝜋 and 𝜃_𝑗(𝜏) satisfies (4.56) for 𝑗=1,2. If there exists a positive integer 𝑁_𝑗⩾0 (𝑗=1,2) such that max_𝜏∈𝐼1𝜏𝜔_𝑗(𝜏)>(2𝑁_𝑗+1)𝜋, then for any 0⩽𝑛⩽𝑁_𝑗, 𝜃_𝑗(𝜏)+2𝑛𝜋 intersects 𝜏𝜔_𝑗(𝜏) at least once at some 𝜏^𝑗𝑛∈𝐼₁. System (2.3) undergoes a Hopf bifurcation at 𝜏^𝑗𝑛 provided 𝑑𝑅𝑒(𝜆(𝜏))/𝑑𝜏|_𝜏=𝜏^𝑗𝑛≠0.

Proof. Assume 1/4+√1/16+1/(2𝛼)>𝑘/(𝛼+𝑘). Then ln((𝑘/𝛼)(1/(1/4+√1/16+1/(2𝛼))−1))<0. By (4.22), 𝐼₁=⎡⎢⎢⎣max0,ln𝑘𝛼1/4+√1/116+1/(2𝛼)−1,ln⎛⎜⎜⎝𝑘⁴√𝛼−1=0,ln𝑘𝛼⎞⎟⎟⎠⎤⎥⎥⎦𝛼⁴√𝛼−1.(4.74) For any 𝜏∈𝐼₁, by Theorem 4.5, 𝜔_𝑗(𝜏)⩾0 for 𝑗=1,2. By Lemma 4.7, there exists 𝜖_𝑗>0 and 𝜃_𝑗(𝜏) such that 𝜖_𝑗⩽𝜃_𝑗(𝜏)⩽𝜋 and 𝜃_𝑗(𝜏) satisfies (4.56). Noting 0∈𝐼₁, min_𝜏∈𝐼1𝜏𝜔_𝑗(𝜏)=0. Assume there exists a positive integer 𝑁_𝑗⩾0 (𝑗=1,2) such that max_𝜏∈𝐼1𝜏𝜔_𝑗(𝜏)>(2𝑁_𝑗+1)𝜋. For any 0⩽𝑛⩽𝑁_𝑗, min_𝜏∈𝐼1𝜏𝜔_𝑗(𝜏)=0⩽2𝑛𝜋 and max_𝜏∈𝐼1𝜏𝜔_𝑗(𝜏)>(2𝑁_𝑗+1)𝜋⩾(2𝑛+1)𝜋. By Corollary 4.10, the conclusion follows.

5. Numerical Results

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

5.1. Variation of Eigenvalues

To study the stability switches of 𝐸₊, 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 parameters 𝐷=Δ=1, 𝑘=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 𝐷=Δ=1, 𝑘=24, and 𝜏=0.5.

Figure 2 

Eigenvalues with the largest real parts of the characteristic equation (4.2) at 𝐸+. 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 𝐷=Δ=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 𝜏=𝜏_𝑐. For any pair (𝛼,𝜏) below that curve, a coexistence equilibrium 𝐸₊ 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 𝐸₊. For any (𝛼,𝜏) outside the closed curve and below 𝜏=𝜏_𝑐, the coexistence equilibrium 𝐸₊ is stable.

Figure 3 

The two-parameter bifurcation diagram of 𝐸+ in 𝜏 and 𝛼 parameter space. Parameters are the same as in Figure 1 (i.e., 𝐷=Δ=1,𝑘=24) except now both 𝜏 and 𝛼 are allowed to vary. For any pair (𝛼,𝜏) on the closed curve, there is a Hopf bifurcation of 𝐸+. Inside the closed curve, there is a periodic solution surrounding 𝐸+. For any (𝛼,𝜏) outside the closed curve and below 𝜏=𝜏𝑐, the coexistence equilibrium 𝐸+ is stable.

5.2. Simulations Demonstrating Hopf Bifurcations

In this section, we illustrate Theorem 4.9 for system (2.3). Take 𝐷=Δ=1 and let 𝜏 vary. We choose parameters 𝛼=100 and 𝑘=100 for Case 1 (see Figures 4–12), and 𝛼=2 and 𝑘=20 for Case 2 (see Figures 13–18).

(a)

(b)

Figure 4 

Critical value of delay 𝜏 at which a Hopf Bifurcation occurs for 𝐷=Δ=1,𝑘=100,𝛼=100.

(a)

(b)

Figure 5 

If the slope of 𝑆𝑗𝑛(𝜏) is nonzero at 𝜏 when 𝑆𝑗𝑛(𝜏)=0 (𝑗=1,2), the transversality condition holds and there is a Hopf bifurcation. The other parameters are the same as for Figure 4.

(a)

(b)

(c)

Figure 6 

Equilibrium 𝐸+(𝜏) is stable when 𝜏=0.02<𝜏10. The other parameters are the same as for Figure 4.

(a)

(b)

(c)

Figure 7 

Time series of a periodic solution, for 𝜏=0.03. The other parameters are the same as for Figure 4.

Figure 8 

The trajectory in phase space of the periodic solution in Figure 7 for 𝜏=0.03.

(a)

(b)

(c)

Figure 9 

Time series of a periodic solution for 𝜏=0.4. The other parameters are the same as for Figure 4.

Figure 10 

The trajectory in phase space of the periodic solution shown in Figure 9.

(a)

(b)

(c)

Figure 11 

The periodic solution disappears at the secondary Hopf bifurcation at 𝜏20≈0.48 and 𝐸+ regains stability. In this figure 𝜏=0.5>𝜏20. The other parameters are the same as for Figure 4.

Figure 12 

Bifurcation diagram as the delay varies. The maximum and minimum amplitude of the 𝑦-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 are 𝐷=Δ=1,𝑘=20,𝛼=2.

Figure 14 

Verification of the transversality condition required for Hopf bifurcation. Parameters are the same as in Figure 13.

(a)

(b)

(c)

Figure 15 

Equilibrium 𝐸+(𝜏) is stable when 𝜏=0.15<𝜏10,1. The other parameters are the same as in Figure 13.

(a)

(b)

(c)

Figure 16 

Time series of a solution with constant initial data 𝑠(0)=0.87, 𝑥(𝑡)=0.077, and 𝑦(𝑡)=0.048 for 𝑡∈[−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.

Figure 17 

The attracting periodic solution shown in Figure 16 in phase space. Note that 𝜏=0.3>𝜏10,1.

(a)

(b)

(c)

Figure 18 

The periodic solution disappears and 𝐸+(𝜏) regains its stability when 𝜏>𝜏20. In this figure, 𝜏=0.74. The other parameters are the same as in Figure 13.

Case 1. Note that 𝐼₁ is given by (4.22). Since ln((𝑘/𝛼)(1/(1/4+√1/16+1/(2𝛼))−1))=−0.03, max0,ln𝑘𝛼1/4+√1/116+1/(2𝛼)−1=0.(5.1) Also, ln(𝑘(√⁴𝛼−1)/𝛼)≈0.77. Therefore 𝐼₁≈[0,0.77]. By Theorem 4.9, 𝜔_𝑗(𝜏) is positive and 𝜃_𝑗(𝜏) satisfies (4.56) for any 𝜏∈𝐼₁ and 𝑗=1,2.
Figure 4 shows that 𝜃_𝑗(𝜏) intersects 𝜏𝜔_𝑗(𝜏) at some 𝜏^𝑗0 with 𝜏¹⁰≈0.022 and 𝜏²⁰≈0.48. We see that 𝜃_𝑗(𝜏)+2𝑛𝜋 has no intersection with 𝜏𝜔_𝑗(𝜏) for 𝑛⩾2 and 𝑗=1,2. By Theorem 4.9, (4.4) has two distinct pairs of pure imaginary eigenvalues 𝜆=±𝑖𝜔_𝑗(𝜏^𝑗0). Next we need to check if Re(d𝜆(𝜏)/d𝜏)|_𝜏=𝜏^𝑗0≠0.
As in Beretta and Kuang [13], we can define
𝑆^𝑗𝑛(𝜏)=𝜏−𝜃^𝑗(𝜏)+2𝑛𝜋𝜔_𝑗(𝜏),for𝑗=1,2,𝑛=0,1,2….(5.2) Any zero 𝜏^𝑗𝑛 of 𝑆^𝑗𝑛(𝜏) is an intersection of 𝜃_𝑗(𝜏)+2𝑛𝜋 and 𝜏𝜔_𝑗(𝜏) and vice versa. By (4.10) in Beretta and Kuang [13] and noting that (𝑞²−𝑝²(𝜏)+2𝛼²(𝜏))²−4(𝛼²(𝜏)−𝑐²(𝜏))>0, we have the relation signRed𝜆(𝜏)=±signd𝜏|||𝑞_𝜏=𝜏^𝑗𝑛₂−𝑝²(𝜏)+2𝛼²(𝜏)²−4𝛼²(𝜏)−𝑐²(𝜏)signd𝑆^𝑗𝑛(𝜏)=±signd𝑆_𝑗𝑛d𝜏||||_{𝜏=𝜏𝑗𝑛}(𝜏)d𝜏||||_{𝜏=𝜏𝑗𝑛},(5.3) where we take + for 𝑗=1 and − for 𝑗=2.
From Figure 5(a), it is observed that 𝑆^1𝑛 has only one zero 𝜏¹⁰≈0.022 at 𝑛=0 with sign{d𝑆¹⁰(𝜏)/d𝜏|_𝜏=𝜏¹⁰}>0. Hence, sign{Re(d𝜆(𝜏)/d𝜏)|_𝜏=𝜏¹⁰}>0. By Theorem 4.9, system (2.3) undergoes a Hopf bifurcation at 𝜏¹⁰. Similarly from Figure 5(b), 𝑆^2𝑛 has only one zero 𝜏²⁰≈0.48 for 𝑛=0 and (2.3) undergoes a Hopf bifurcation at 𝜏²⁰.
Next we used MATLAB to simulate solutions of model (2.3) for several values of 𝜏. For each fixed delay 𝜏, we chose initial data 𝑠(𝑡)=𝑠₊(𝜏)−0.01, 𝑥(𝑡)=𝑥₊(𝜏)+0.01, and 𝑦(𝑡)=𝑦₊(𝜏)+0.001 for 𝑡∈[−𝜏,0]. From Figure 6, we can see that the equilibrium 𝐸₊ is stable if 𝜏=0.02<𝜏¹⁰. As delay 𝜏 increases past 𝜏¹⁰≈0.022, where a Hopf bifurcation occurs, a pair of complex eigenvalues of (4.4) enters the right-half plane. The equilibrium 𝐸₊ loses its stability and a periodic solution bifurcates from 𝐸₊ (see Figures 7 and 8). As we increase the delay further to 𝜏=0.4<𝜏²⁰, the periodic solution still exists and remains stable (see Figures 9 and 10). However, as the delay 𝜏 increases further, past 𝜏²⁰, the stable periodic solution disappears in a second Hopf bifurcation, and 𝐸₊ regains stability (see Figure 11). We provide a bifurcation diagram illustrating the change in dynamics as 𝜏 varies (see Figure 12). For any 𝜏∈(𝜏¹⁰,𝜏²⁰), there is an orbitally asymptotically stable periodic solution.

Case 2. Take 𝑘=20 and 𝛼=2. For such parameters, 𝐼₂≈[0,0.85]. By Theorem 4.9, 𝜔₁(𝜏) is positive and 𝜃₁(𝜏) satisfies (4.56) for 𝑗=1 and 𝜏∈𝐼₂. Figure 13 shows that 𝜃₁(𝜏) intersects 𝜏𝜔₁(𝜏) twice. To distinguish these intersections, denote them as 𝜏^10,1≈0.2 and 𝜏^10,2≈0.7. On the other hand, 𝜃₁(𝜏)+2𝜋 has no intersection with 𝜏𝜔₁(𝜏).
From Figure 14,
sign⎧⎪⎨⎪⎩d𝑆¹⁰(𝜏)d𝜏||||_{𝜏=𝜏10,1}⎫⎪⎬⎪⎭>0,sign⎧⎪⎨⎪⎩d𝑆¹¹(𝜏)d𝜏||||_{𝜏=𝜏10,2}⎫⎪⎬⎪⎭<0.(5.4) By (5.3), signRed𝜆(𝜏)d𝜏|||_𝜏=𝜏^10,1>0,signRed𝜆d(𝜏𝜏)|||_𝜏=𝜏^10,2<0.(5.5) By Theorem 4.5, system (2.3) undergoes a Hopf bifurcation at 𝜏=𝜏^10,1 and at 𝜏=𝜏^10,2. For 𝜏 less than 𝜏^10,1, the equilibrium 𝐸₊ is asymptotically stable (see Figure 15). For 𝜏 greater than 𝜏^10,1, but less than 𝜏^10,2, there is an orbitally asymptotically stable periodic solution surrounding the equilibrium 𝐸₊ (see Figures 16 and 17). At 𝜏=𝜏^10,2, there is a second Hopf bifurcation, where the periodic solution coalesces with 𝐸₊. For 𝜏>𝜏^10,2, the periodic orbit no longer exists and 𝐸₊ regains stability (see Figure 18) until it disappears when 𝜏>𝜏_𝑐.

Appendix

Preliminary Results

To establish the existence of periodic solutions in autonomous delay differential equations, one of the simplest ways is through Hopf Bifurcation. Below is a general Hopf Bifurcation theorem for delay differential equations due to De Oliveira [14]. Before stating the theorem we require some notation.

Consider a one parameter family of neutral delay differential equations:

dd𝑡𝐷𝛼,𝑥^𝑡−𝑔𝛼,𝑥_𝑡=𝐿𝛼,𝑥_𝑡+𝑓𝛼,𝑥_𝑡,𝛼∈ℝ,(A.1) where 𝐷,𝐿,𝑓, and 𝑔 are continuously differentiable in 𝛼 and 𝑥_𝑡∈ℂ([−𝑟,0],ℝ^𝑛) (𝑟 is a constant), 𝑓(𝛼,0)=𝑔(𝛼,0), 𝜕𝑓(𝛼,0)/𝜕𝑥_𝑡=𝜕𝑔(𝛼,0)/𝜕𝑥_𝑡=0, 𝐷(𝛼,𝑥_𝑡) and 𝐿(𝛼,𝑥_𝑡) are linear in 𝑥_𝑡, and

𝐷𝛼,𝑥_𝑡=^∞_𝑘=0𝐿𝛼,𝑥_t𝐴_𝑘(𝛼)𝑥𝑡−𝑟_𝑘(𝛼)+⁰⁻¹𝐴(𝛼,𝜃)𝑥(𝑡+𝜃)d𝜃,=^∞_𝑘=0𝐴_𝑘(𝛼)𝑥𝑡−𝑟_𝑘(𝛼)+⁰⁻¹𝐴(𝛼,𝜃)𝑥(𝑡+𝜃)d𝜃,(A.2) for 𝑥_𝑡∈ℂ([−𝑟,0],ℝ^𝑛). Assume 𝛼∈ℝ, where 𝑟₀(𝛼)=0, 𝑟_𝑘(𝛼)∈(0,1], and 𝐴_𝑘(𝛼), 𝐵_𝑘(𝛼), 𝐴(𝛼,𝜃), and 𝐵(𝛼,𝜃) satisfy

^∞_𝑘=0||𝐴_𝑘(𝛼)||+||𝐵_𝑘(𝛼)||+⁰⁻¹||𝐴(𝛼,𝜃)||+||𝐵(𝛼,𝜃)||d𝜃<∞.(A.3) It is easy to see that the characteristic matrix

Δ(𝛼,𝜆)=𝜆𝐷𝛼,𝑒^𝜆⋅𝐼−𝐿𝛼,𝑒^𝜆⋅𝐼(A.4) is continuously differentiable in 𝛼∈ℝ and Δ(𝛼,𝜆) is an entire function of 𝜆. Making the following assumptions on (A.1).

(𝑆₁) There exist constants 𝑎>0, 𝑏>0 such that, for all complex values 𝜆 such that |𝑅𝑒𝜆|<𝑎 and all 𝛼∈ℝ, the following inequalities hold: |||||det_∞_𝑘=0𝐴_𝑘(𝛼)𝑒^{−𝜆𝑟𝑘(𝛼)}|||||det_∞|||||⩾𝑏,_𝑘=0𝐴_𝑘(𝛼)𝑒^{−𝜆𝑟𝑘(𝛼)}+⁰⁻¹𝐴(𝛼,𝜃)𝑒^𝜆𝜃d𝜃|||||⩾𝑏.(A.5)(𝑆₂) The characteristic equation detΔ(𝛼,𝜆)=0 has, for 𝛼=𝛼₀, a simple purely imaginary root 𝜆₀=𝑖𝑣₀, 𝑣₀>0, and no root of detΔ(𝛼₀,𝜆)=0, other than ±𝑖𝑣₀, is an integral multiple of 𝜆₀.(𝑆₃)𝑅𝑒(𝜕𝜆(𝛼₀)/𝜕𝛼)≠0.

Now we are ready to state the Hopf bifurcation theorem for (A.1).

Theorem A.1 (Hopf bifurcation theorem, see Kuang [10, page 60]). In (A.1), assume that (𝑆₁)-(𝑆₃) hold. Then there is an 𝜖>0 such that, for 𝑎∈ℝ, |𝑎|⩽𝜖, there are functions 𝛼(𝑎)∈ℝ, 𝜔(𝑎)∈ℝ, 𝛼(0)=𝛼₀, 𝜔(0)=2𝜋/𝑣₀, such that (A.1) has an 𝜔(𝛼)-periodic solution 𝑥^∗(𝑎)(𝑡), that is continuously differentiable in 𝑡, and 𝑎 with 𝑥^∗(0)=0. Furthermore, for |𝛼−𝛼_𝑜|<𝜖, |𝜔−(2𝜋/𝑣₀)|<𝜖, every 𝜔-periodic solution 𝑥(𝑡) of (A.1) with |𝑥(𝑡)|<𝜖 must be of this type, except for a translation in phase; that is, there exists 𝑎∈(−𝜖,𝜖) and 𝑏∈ℝ such that 𝑥(𝑡)=𝑥^∗(𝑎)(𝑡+𝑏) for all 𝑡∈ℝ.

The following lemma is usually called the Fluctuation Lemma. For a proof, see Hirsh et al. [15].

Lemma A.2. Let 𝑓∶ℝ⁺→ℝ be a differentiable function. If liminf_𝑡→∞𝑓(𝑡)<limsup_𝑡→∞𝑓(𝑡), then there are sequences 𝑡_𝑚↑∞ and 𝑠_𝑚↑∞ such that for all 𝑚𝑓𝑓̇̇𝑠𝑡_𝑚=0,𝑓𝑡_𝑚⟶limsup_𝑡→∞𝑓(𝑡)as𝑚⟶∞,_𝑚=0,𝑓𝑠_𝑚⟶liminf_𝑡→∞𝑓(𝑡)as𝑚⟶∞,.(A.6)

The proof of the following useful lemma can be found in [16].

Theorem A.3. Let 𝑎∈(−∞,∞) and 𝑓∶[𝑎,∞)→ℝ be a differentiable function. If lim_𝑡→∞𝑓(𝑡) exists (finite) and the derivative function 𝑓̇(𝑡) is uniformly continuous on (𝑎,∞), then lim_𝑡→∞𝑓̇(𝑡)=0.

Acknowledgment

This research is partially supported by NSERC.