# Stability and Bifurcation Analysis of a Delayed Leslie-Gower Predator-Prey System with Nonmonotonic Functional Response

Authors: Jiao Jiang; Yongli Song
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

A delayed Leslie-Gower predator-prey model with nonmonotonic functional response is studied. The existence and local stability of the positive equilibrium of the system with or without delay are completely determined in the parameter plane. Using the method of upper and lower solutions and monotone iterative scheme, a sufficient condition independent of delay for the global stability of the positive equilibrium is obtained. Hopf bifurcations induced by the ratio of the intrinsic growth rates of the predator and prey and by delay, respectively, are found. Employing the normal form theory, the direction and stability of Hopf bifurcations can be explicitly determined by the parameters of the system. Some numerical simulations are given to support and extend our theoretical results. Two limit cycles enclosing an equilibrium, one limit cycle enclosing three equilibria and different types of heteroclinic orbits such as connecting two equilibria and connecting a limit cycle and an equilibrium are also found by using analytic and numerical methods.


---

## Body

## 1. Introduction

Predation is one of the most fundamental interactions between different species and is universal in nature. Recent studies have shown that predator-prey relationships (not competitive or mutualistic relationships) provide the necessary stability for almost infinite numbers of species to exist in ecosystems and make the rich biodiversity of complex ecosystems possible [1]. The mathematical model and the related dynamics are very useful for qualitatively understanding the interaction between predators and preys. Thus, the predator-prey dynamics has been extensively studied and continue to be of interest to both applied mathematicians and ecologists.The pioneering predator-prey model was initially proposed by Lotka [2] and Volterra [3] (known as Lotka-Volterra predator-prey model) and successfully captured the oscillations in populations of a predator and its prey. Holling extended this model in [4–6] where various types of functional responses were proposed to better understand the components of predator-prey interactions. The Lotka-Volterra predator-prey model with Holling functional response has been a cornerstone in the field of mathematical biology, and much literature has been devoted to study variants of these equations.Another important type of well-known predator-prey models was initially proposed and studied by Leslie and Gower [7–9]. Compared with Lotka-Volterra model, the novel feature of Leslie-Gower model is that both interacting species are assumed to grow according to the logistic law, the environmental carrying capacity of the predator is not a constant but a function of the available prey quantity, that there are upper limits to the rates of increase of both prey and predator, which are not recognized in the Lotka-Volterra model. Later, Holling functional response was introduced to Leslie-Gower model by May [10] and the corresponding model is described by the following ordinary differential equations:
(1)x˙=rx(1-xK)-p(x)y,y˙=s(1-ynx)y,
where x(t) and y(t) denote the prey and predator densities, respectively, p(x) represents the Holling type functional response [6], the positive constants r and s are the intrinsic growth rates of the prey and predators, respectively, the carrying capacity for the prey is a positive constant K, n is a measure of the food quality that the prey provides for conversion into predator births, and then nx is considered as the carrying capacity for the predator. The model (1) is also known as the Holling-Tanner predator-prey model [11–13]. When the functional response p(x) is Holling types I, II, and III functions respectively, it has been proven that the dynamics of model (1) are quite interesting and the global stability of the positive equilibrium, the existence of limit cycles, and so forth have been investigated by many researchers (see, e.g., [12–16] and references therein). Holling type I, II, and III functional responses are similar in that each per capita rate of predation increases with prey density to a maximum [4] and it is shown that the model (1) with the type I, II, and III functional responses exhibits qualitatively similar bifurcation and stability behavior [17]. Holling type IV functional response described by Holling in [18] incorporates prey interference with predation. The main characteristic of Holling type IV functional response is that the per capita predation rate increases with prey density to a maximum at a critical prey density beyond which it decreases, while for other Holling functional responses, the per capita predation rate always increases with prey density. Collings [17] has also numerically shown that the dynamics of the model (1) with Holling type IV functional response is very different from that of the model (1) with other Holling functional responses at higher levels of prey interference. The qualitative analysis for the model (1) with a simplified Holling type IV functional response, proposed by Sokol and Howell [19] as follows:
(2)p(x)=mxa+x2,
has been studied by Li and Xiao [20]. They mainly studied the case when the model has two equilibria: one is an unstable multiple focus with multiplicity one and the other is a cusp of codimension 2 and the dynamics of the model near the small neighborhoods of these two equilibria.On the other hand, the introduction of time delay into the population model is more realistic to model the interaction between the predator and prey populations and the population models with time delay are of current research interest in mathematical biology [21, 22]. There is extensive literature about the effects of delay on the dynamics of predator-prey models. We refer to a recent survey paper by Ruan [22] and the references cited therein for studies on delayed Lotka-Volterra predator-prey models with Holling functional response. The effects of delay on the dynamics of the Leslie-Gower predator-prey model with and its modified version have been recently studied by many authors (see, e.g., [23–26] and references therein). However, most of them considered the case of the monotonic functional response, under which the corresponding system has only one positive equilibrium. When the nonmonotonic functional response (Holling type IV) is introduced to the Leslie-Gower predator-prey model, there is a possibility of the existence of one, two or three positive equilibria and the dynamics is much more complex [17, 20]. In [27], the authors considered the following delayed Leslie-Gower model with nonmonotonic functional response:
(3)x′(t)=rx(t)(1-x(t)k1)-mx(t)a+x2(t)y(t),y′(t)=y(t)[s(1-hy(t-τ)x(t-τ)+k2)],
where τ is incorporated because of the negative feedback of the predator density. Under the assumption that the positive equilibrium exists and it is stable for τ=0, the local stability and delay-induced Hopf bifurcations have been investigated. The supercritical stable Hopf bifurcations were also found by normal form theory and numerical simulation.In this paper, following the work of Li and Xiao [20] and Lian and Xu [27], we study how the hunting delay affects the dynamics of the Leslie-Gower predator-prey model nonmonotonic functional response. Introducing the hunting delay into the Leslie-Gower predator-prey model with nonmonotonic functional response, we have the following delay differential equations:
(4)x′(t)=rx(t)(1-x(t)K)-mx(t)y(t-τ)x2(t)+a,y′(t)=s(1-y(t)nx(t))y(t),
where the biologic meaning of the parameters r,  s,  n is the same as in model (1), m is the maximal predator per capita consumption rate, that is, the maximum number of prey that can be eaten by a predator in each time unit, and a is the number of prey necessary to achieve one-half of the maximum rate m. To the best of our knowledge, although the dynamics of the delayed Lotka-Volterra predator-prey model with nonmonotonic functional response has been studied in [28, 29], there are a few research papers on the dynamics of the delayed Leslie-Gower predator-prey model (4). Our goal is to obtain the explicit condition for the existence of the positive equilibrium, to investigate the local and global stability, and to study how the hunting delay affects the stability of the positive equilibrium and induces the Hopf bifurcations and their properties.The paper is organized as follows. In Section2, we study the explicit condition for the existence of the positive equilibria and their stability and the Hopf bifurcation induced by the intrinsic growth rate of the predator for the system without delay. In Section 3, we investigate the effect of the delay on the dynamics of the system. The local stability is studied by analyzing the related characteristic equation and the global stability is studied by the method of upper and lower solutions and monotone iterative scheme, and the direction and stability of the delay-induced Hopf bifurcation are determined by using the normal form theory. Numerical simulations are performed to illustrate and extend the obtained results in Section 4. A brief discussion is given in Section 5 to conclude the paper.


## 2. Linear Stability and Hopf Bifurcation Analysis for the System without Delay

In order to reduce the number of parameters and simplify the calculus in system (4), introducing new variables
(5)t¯=rt,x¯(t¯)=x(t)K,y¯(t¯)=my(t)rK2,τ¯=rτ
and then dropping the bars, we derive
(6)x˙(t)=x(t)(1-x(t))-x(t)y(t-τ)x2(t)+α,y˙(t)=γ(1-y(t)βx(t))y(t),
where α=a/K2,  β=mn/Kr, and γ=s/r are positive parameters. When τ=0, system (6) becomes the following ordinary differential equations:
(7)x˙(t)=x(t)(1-x(t))-x(t)y(t)x2(t)+α,y˙(t)=γ(1-y(t)βx(t))y(t).
Systems (6) and (7) have the same equilibria, but the stability of system (6) is affected by delay. We first consider the existence of equilibria and their stability for system (7).


### 2.1. Existence of Equilibria and Their Stability

From the biological viewpoint, we are only interested in the dynamics of system (7) in the closed first quadrant R+2 in the (x,y) plane. Letting E*(x*,y*) denote the equilibrium of system (7), then (x*,y*) are the positive real roots of the following equations:
(8)1-x=yx2+α,y=βx.We first have the following results on the existence of the positive equilibrium of system (7).Lemma 1.
(i) Ifα≥1/27, then system (7) has a unique positive equilibrium.
(ii) If0<α<1/27, then system (7) has three positive equilibria for β1(α)<β<β2(α), two positive equilibria for β=β1(α) or β=β2(α), and a unique positive equilibrium for β<β1(α) or β>β2(α), where βj(α)=-3Tj2+2Tj-α,  j=1,2, with
(9)T1=1+cos(arccos(-1+54α)/3)6-3sin(arccos(-1+54α)/3)6,T2=1+cos(arccos(-1+54α)/3)6+3sin(arccos(-1+54α)/3)6.Proof.
(i) From (8), we have (1-x)(x2+α)=βx. Let y1(x)=(1-x)(x2+α), and y2(x)=βx. Notice that y1(0)=α>0 and y1(x) is decreasing as the x-value increases for α≥1/3 since y1′(x)=-3x2+2x-α≤0 for α≥1/3. So, when α≥1/3, the two curves y1(x) and y2(x) have only one interaction point for x>0 and β>0.
When0<α<1/3, the curve y1(x) has one minimum at x=m* and one maximum at x=m* and decreases in the interval x∈(0,m*)∪(m*,1) and increases in the interval (m*,m*), where
(10)m*=1-1-3α3,m*=1+1-3α3.
The tangent point of two curves y1(x) and y2(x) is the key point for determining how many intersection points of these two curves. The tangent point of two curves y1(x) and y2(x) can be determined by the roots of the following equations
(11)(1-x)(α+x2)=βx,-3x2+2x-α=β,
which yields
(12)2x3-x2+α=0.
It is easy to verify that the function p(x)=2x3-x2+α has a minimum at x=1/3 for x>0. Noticing that p(1/3)=0 is equivalent to α=1/27, we can deduce that (12) has no positive real roots for 1/27<α<1/3, one positive real root for α=1/27, and two different positive real roots for 0<α<1/27. Thus, when 1/27<α<1/3, the straight line y2(x) is not tangent to the curve y1(x) for any x>0 and β>0. When α=1/27, the straight line y2(x) is tangent to the curve y1(x) only at β=8/27 and x=1/3, as shown in Figure 1(a). Consequently, when 1/27≤α<1/3, the two curves y1(x) and y2(x) have only one interaction point for any x>0 and β>0. This completes the proof of conclusion (i).
(ii) When0<α<1/27, it follows from Shengjin’s theorem [30] that the two positive real roots of (12) are T1 and T2 defined by (9). This, together with (11), implies that two curves y1(x) and y2(x) are tangent at points βj=-3Tj2+2Tj-α,  j=1,2, as shown in Figure 1(b) for α=1/40. This completes the proof of conclusion (ii).Diagram of positive equilibrium. (a) Whenα≥1/27, there exists only one positive equilibrium for any β>0. Here, α=1/27, and β=8/27. (b) When 0<α<1/27, there are three cases: (i) only one positive equilibrium for β<β1 or β>β2, (ii) two positive equilibria for β=β1 or β=β2, and (iii) three positive equilibria for β1<β<β2. Here, α=1/40,  β1=0.2599, and β2=0.2782.
(a)
(b)According to Lemma1, the distribution of the positive equilibria of system (7) in the α-β plane can be illustrated in Figure 2.Bifurcation diagram of positive equilibrium in theα-β plane. System (7) has three positive equilibria in the inner region enclosed by the curves β1 and β2, a unique positive equilibrium outside this region and on the point P(1/27,8/27) and two positive equilibria on the two curves on the boundary curves β1 and β2 except the point P(1/27,8/27).
(a)
(b)The Jacobian matrix of system (7) at the positive equilibrium point E* is
(13)J=(-r1-r2βγ-γ),
where
(14)r1=-x*y1′(x*)(x*)2+α,r2=x*(x*)2+α>0.
Then the characteristic equation at E* is
(15)λ2+(r1+γ)λ+γ(r1+r2β)=0.Lettingβ1* and β2* denote the slopes of the straight line y2(x) crossing the maximum and minimum points of y1(x), respectively, we have
(16)β1*=(1-m*)(α+(m*)2)m*,β2*=(1-m*)(α+m*2)m*,
where m* and m* are defined by (10). When α<1/3, we have β1*<β2*, as shown in Figure 2. The curves β1,  β2,  β1*,  β2* divide the α-β plane into four regions Di,  i=1,2,3,4, as shown in Figure 2. From Lemma 1, we know that system (7) has only one positive equilibrium in the regions D1 and D2, two positive equilibria in the curves β1(α) and β2(α) with α≤1/27, three positive equilibria in the regions D3 and D4. Notice that
(17)r1=-x*y1′(x*)(x*)2+α,r1+r2β=x*(x*)2+α(β-y1′(x*)),
and then by a direct but tedious analysis, we have the following results on the signs of r1 and r1+r2β, which determine the stability of the positive equilibrium.In the regionsD1 and D2, we have
(18)r1{>0,for(α,β)∈D1,=0,for(α,β)lyingonthecurvesβ1*(α)orβ2*(α),<0,for(α,β)∈D2,r1+r2β>0.In the regionsD3 and D4, denote the three positive equilibria by Ei*(xi*,yi*),  i=1,2,3 with x1*<x2*<x3*. Notice that x1*<T1,  T1<x2*<T2,  x3*>T2. Then we have
(19)r1{<0,for(α,β)∈D3,allthreeequilibriaEi*,i=1,2,3,<0,for(α,β)∈D4,theequilibriaE1*andE2*,>0,for(α,β)∈D4,theequilibriaE3*,(20)r1+r2β{<0,forE2*,>0,forE1*andE3*.
In addition, notice that if r1≥0, then r1+γ>0 for any γ>0. However, when r1<0, we can choose γ as a parameter such that r1+γ≥0 for γ≥-r1 and r1+γ<0 for 0<γ<-r1. Thus, by (18), (19), and (20), we have the following results on the distribution of roots of (15).Lemma 2.
Assume thatT1,  T2 are defined by (9) and γ0=x*y1′(x*)/((x*)2+α).(i)
In the regionD1, for the unique equilibrium E*, all roots of (15) have negative real parts for any γ>0.(ii)
In the regionD2, for the unique equilibrium E*, all roots of (15) have negative real parts for γ>γ0, a pair of purely imaginary roots ±iγ0(r1+r2β) at γ=γ0, and (15) has at least a root with positive real part for 0<γ<γ0.(iii)
In the regionD3, for the positive equilibrium E2*, (15) has at least a root with positive real part for any γ>0; for other two positive equilibria E1* and E3*, all roots of (15) have negative real parts for γ>γ0, a pair of purely imaginary roots ±iγ0(r1+r2β) at γ=γ0, and (15) has at least a root with positive real part for 0<γ<γ0.(iv)
In the regionD4, for the positive equilibrium E2*, (15) has at least a root with positive real part for any γ>0; for the positive equilibrium E1*, all roots of (15) have negative real parts for γ>γ0, a pair of purely imaginary roots ±iγ0(r1+r2β) at γ=γ0, and (15) has at least a root with positive real part for 0<γ<γ0; for the positive equilibrium E3*, all roots of (15) have negative real parts for any γ>0.Takingγ as a parameter and letting λ is a root of (15) such that Re(λ(γ0))=0 and Im(λ(γ0))=γ0(r1+r2β). A straight calculation yields the following transversality condition:
(21)dRe(λ)dγ|γ=γ0=Re{-λ+r1+r2β2λ+r1+γ}|γ=γ0=Re{-12-r1+r2β2λ}|γ=γ0=-12<0.From Lemma2 and the transversality condition (21), the following theorem on the stability and bifurcation follows immediately.Theorem 3.
Assuming thatT1,  T2 are defined by (9) and γ0=x*y1′(x*)/((x*)2+α).(i)
In the regionD1, the unique equilibrium E* is asymptotically stable for any γ>0.(ii)
In the regionD2, the unique equilibrium E* is asymptotically stable for γ>γ0 and unstable for 0<γ<γ0, and system (7) undergoes Hopf bifurcation at γ=γ0.(iii)
In the regionD3, the positive equilibrium E2* is unstable for any γ>0; other two positive equilibria E1* and E3* are asymptotically stable for γ>γ0 and unstable for 0<γ<γ0, the and system (7) undergoes Hopf bifurcation near neighborhood of each of these two equilibria at γ=γ0.(iv)
In the regionD4, the positive equilibrium E2* is unstable for any γ>0; the positive equilibrium E1* is asymptotically stable for γ>γ0 and unstable for 0<γ<γ0, and system (7) undergoes Hopf bifurcation near neighborhood of E1* at γ=γ0; the positive equilibrium E3* is asymptotically stable for any γ>0.


### 2.2.   Direction and Stability of Hopf Bifurcation Associated withγ0

It follows from Theorem3 that when (α,β)∈Di,  i=2,3,4, Hopf bifurcation occurs at γ=γ0. When Hopf bifurcation occurs, the corresponding characteristic equation has a pair of purely imaginary roots at γ=γ0 and thus the linearization is not sufficient to determine the stability of the corresponding equilibrium, say, E*. In this subsection, by calculating the focal value we determine the stability of the equilibrium E* at γ=γ0, which is directly associated with the direction and stability of Hopf bifurcation associated with γ0.Assume thatλ1,2=η(γ)±iϖ(γ) are a pair of complex conjugates of (15) such that η(γ0)=0,  ϖ(γ0)=γ0(r1+r2β). Then the eigenvector corresponding to λ1,2 is
(22)(1-r1+η(γ)±iϖ(γ)r2).
Under the following coordinate transformation:
(23)(z1z2)=(10MN)(uv),
where M=-(r1+η(γ))/r2, N=ϖ(γ)/r2, system (7) becomes
(24)(u˙(t)v˙(t))=J(γ)(u(t)v(t))+(F(1)(u,v,γ)F(2)(u,v,γ)),
where
(25)J(γ)=(η(γ)-ϖ(γ)ϖ(γ)η(γ))
with
(26)F(1)(u,v,γ)=(12f20(1)+Mf11(1))u2+Nf11(1)uv+(16f30(1)+12Mf21(1))u3+12Nf21(1)u2v+O(|u,v|4),F(2)(u,v,γ)=-MNF(1)(u,v,γ)+1NF^(2)(u,v,γ),
where
(27)F^(2)(u,v,γ)=(12f20(2)+Mf11(2)+12M2f02(2))u2+(Nf11(2)+MNf02(2))uv+12N2f02(2)v2+(16f30(2)+12Mf21(2)+12M2f12(2))u3+(12Nf21(2)+MNf12(2))u2v+12N2f12(2)uv2+O(|u,v|4).
Find fij(1) and fij(2) at (31). By the formula for the third focal value of a multiple focus in [31], we have
(28)a3(γ0)=116(Fuuu(1)+Fuvv(1)+Fuuv(2)+Fvvv(2))|(0,0,γ0)+116ϖ(γ0)(Fuv(1)(Fuu(1)+Fvv(1))-Fuv(2)(Fuu(2)+Fvv(2))-Fuu(1)Fuu(2)+Fvv(1)Fvv(2))|(0,0,γ0),
where
(29)Fuuu(1)=f30(1)+3Mf21(1),Fuvv(1)=0,Fuuv(2)=-Mf21(1)+f21(2)+2Mf12(2),Fvvv(2)=0,Fuv(1)=Nf11(1),Fuu(1)=f20(1)+2Mf11(1),Fvv(1)=0,Fuv(2)=-Mf11(1)+f11(2)+Mf02(2),Fvv(2)=N2f02(2),Fuu(2)=-MN(12f20(1)+Mf11(1))+1N(12f20(2)+Mf11(2)+12M2f02(2)).It is well known from the result in [31] that if a3(γ0)<0, then the equilibrium E* of system (7) is a stable multiple focus of multiplicity 1 and the corresponding Hopf bifurcation is supercritical, and if a3(γ0)>0, then the equilibrium E* of system (7) is an unstable multiple focus of multiplicity 1 and the corresponding Hopf bifurcation is subcritical.


## 2.1. Existence of Equilibria and Their Stability

From the biological viewpoint, we are only interested in the dynamics of system (7) in the closed first quadrant R+2 in the (x,y) plane. Letting E*(x*,y*) denote the equilibrium of system (7), then (x*,y*) are the positive real roots of the following equations:
(8)1-x=yx2+α,y=βx.We first have the following results on the existence of the positive equilibrium of system (7).Lemma 1.
(i) Ifα≥1/27, then system (7) has a unique positive equilibrium.
(ii) If0<α<1/27, then system (7) has three positive equilibria for β1(α)<β<β2(α), two positive equilibria for β=β1(α) or β=β2(α), and a unique positive equilibrium for β<β1(α) or β>β2(α), where βj(α)=-3Tj2+2Tj-α,  j=1,2, with
(9)T1=1+cos(arccos(-1+54α)/3)6-3sin(arccos(-1+54α)/3)6,T2=1+cos(arccos(-1+54α)/3)6+3sin(arccos(-1+54α)/3)6.Proof.
(i) From (8), we have (1-x)(x2+α)=βx. Let y1(x)=(1-x)(x2+α), and y2(x)=βx. Notice that y1(0)=α>0 and y1(x) is decreasing as the x-value increases for α≥1/3 since y1′(x)=-3x2+2x-α≤0 for α≥1/3. So, when α≥1/3, the two curves y1(x) and y2(x) have only one interaction point for x>0 and β>0.
When0<α<1/3, the curve y1(x) has one minimum at x=m* and one maximum at x=m* and decreases in the interval x∈(0,m*)∪(m*,1) and increases in the interval (m*,m*), where
(10)m*=1-1-3α3,m*=1+1-3α3.
The tangent point of two curves y1(x) and y2(x) is the key point for determining how many intersection points of these two curves. The tangent point of two curves y1(x) and y2(x) can be determined by the roots of the following equations
(11)(1-x)(α+x2)=βx,-3x2+2x-α=β,
which yields
(12)2x3-x2+α=0.
It is easy to verify that the function p(x)=2x3-x2+α has a minimum at x=1/3 for x>0. Noticing that p(1/3)=0 is equivalent to α=1/27, we can deduce that (12) has no positive real roots for 1/27<α<1/3, one positive real root for α=1/27, and two different positive real roots for 0<α<1/27. Thus, when 1/27<α<1/3, the straight line y2(x) is not tangent to the curve y1(x) for any x>0 and β>0. When α=1/27, the straight line y2(x) is tangent to the curve y1(x) only at β=8/27 and x=1/3, as shown in Figure 1(a). Consequently, when 1/27≤α<1/3, the two curves y1(x) and y2(x) have only one interaction point for any x>0 and β>0. This completes the proof of conclusion (i).
(ii) When0<α<1/27, it follows from Shengjin’s theorem [30] that the two positive real roots of (12) are T1 and T2 defined by (9). This, together with (11), implies that two curves y1(x) and y2(x) are tangent at points βj=-3Tj2+2Tj-α,  j=1,2, as shown in Figure 1(b) for α=1/40. This completes the proof of conclusion (ii).Diagram of positive equilibrium. (a) Whenα≥1/27, there exists only one positive equilibrium for any β>0. Here, α=1/27, and β=8/27. (b) When 0<α<1/27, there are three cases: (i) only one positive equilibrium for β<β1 or β>β2, (ii) two positive equilibria for β=β1 or β=β2, and (iii) three positive equilibria for β1<β<β2. Here, α=1/40,  β1=0.2599, and β2=0.2782.
(a)
(b)According to Lemma1, the distribution of the positive equilibria of system (7) in the α-β plane can be illustrated in Figure 2.Bifurcation diagram of positive equilibrium in theα-β plane. System (7) has three positive equilibria in the inner region enclosed by the curves β1 and β2, a unique positive equilibrium outside this region and on the point P(1/27,8/27) and two positive equilibria on the two curves on the boundary curves β1 and β2 except the point P(1/27,8/27).
(a)
(b)The Jacobian matrix of system (7) at the positive equilibrium point E* is
(13)J=(-r1-r2βγ-γ),
where
(14)r1=-x*y1′(x*)(x*)2+α,r2=x*(x*)2+α>0.
Then the characteristic equation at E* is
(15)λ2+(r1+γ)λ+γ(r1+r2β)=0.Lettingβ1* and β2* denote the slopes of the straight line y2(x) crossing the maximum and minimum points of y1(x), respectively, we have
(16)β1*=(1-m*)(α+(m*)2)m*,β2*=(1-m*)(α+m*2)m*,
where m* and m* are defined by (10). When α<1/3, we have β1*<β2*, as shown in Figure 2. The curves β1,  β2,  β1*,  β2* divide the α-β plane into four regions Di,  i=1,2,3,4, as shown in Figure 2. From Lemma 1, we know that system (7) has only one positive equilibrium in the regions D1 and D2, two positive equilibria in the curves β1(α) and β2(α) with α≤1/27, three positive equilibria in the regions D3 and D4. Notice that
(17)r1=-x*y1′(x*)(x*)2+α,r1+r2β=x*(x*)2+α(β-y1′(x*)),
and then by a direct but tedious analysis, we have the following results on the signs of r1 and r1+r2β, which determine the stability of the positive equilibrium.In the regionsD1 and D2, we have
(18)r1{>0,for(α,β)∈D1,=0,for(α,β)lyingonthecurvesβ1*(α)orβ2*(α),<0,for(α,β)∈D2,r1+r2β>0.In the regionsD3 and D4, denote the three positive equilibria by Ei*(xi*,yi*),  i=1,2,3 with x1*<x2*<x3*. Notice that x1*<T1,  T1<x2*<T2,  x3*>T2. Then we have
(19)r1{<0,for(α,β)∈D3,allthreeequilibriaEi*,i=1,2,3,<0,for(α,β)∈D4,theequilibriaE1*andE2*,>0,for(α,β)∈D4,theequilibriaE3*,(20)r1+r2β{<0,forE2*,>0,forE1*andE3*.
In addition, notice that if r1≥0, then r1+γ>0 for any γ>0. However, when r1<0, we can choose γ as a parameter such that r1+γ≥0 for γ≥-r1 and r1+γ<0 for 0<γ<-r1. Thus, by (18), (19), and (20), we have the following results on the distribution of roots of (15).Lemma 2.
Assume thatT1,  T2 are defined by (9) and γ0=x*y1′(x*)/((x*)2+α).(i)
In the regionD1, for the unique equilibrium E*, all roots of (15) have negative real parts for any γ>0.(ii)
In the regionD2, for the unique equilibrium E*, all roots of (15) have negative real parts for γ>γ0, a pair of purely imaginary roots ±iγ0(r1+r2β) at γ=γ0, and (15) has at least a root with positive real part for 0<γ<γ0.(iii)
In the regionD3, for the positive equilibrium E2*, (15) has at least a root with positive real part for any γ>0; for other two positive equilibria E1* and E3*, all roots of (15) have negative real parts for γ>γ0, a pair of purely imaginary roots ±iγ0(r1+r2β) at γ=γ0, and (15) has at least a root with positive real part for 0<γ<γ0.(iv)
In the regionD4, for the positive equilibrium E2*, (15) has at least a root with positive real part for any γ>0; for the positive equilibrium E1*, all roots of (15) have negative real parts for γ>γ0, a pair of purely imaginary roots ±iγ0(r1+r2β) at γ=γ0, and (15) has at least a root with positive real part for 0<γ<γ0; for the positive equilibrium E3*, all roots of (15) have negative real parts for any γ>0.Takingγ as a parameter and letting λ is a root of (15) such that Re(λ(γ0))=0 and Im(λ(γ0))=γ0(r1+r2β). A straight calculation yields the following transversality condition:
(21)dRe(λ)dγ|γ=γ0=Re{-λ+r1+r2β2λ+r1+γ}|γ=γ0=Re{-12-r1+r2β2λ}|γ=γ0=-12<0.From Lemma2 and the transversality condition (21), the following theorem on the stability and bifurcation follows immediately.Theorem 3.
Assuming thatT1,  T2 are defined by (9) and γ0=x*y1′(x*)/((x*)2+α).(i)
In the regionD1, the unique equilibrium E* is asymptotically stable for any γ>0.(ii)
In the regionD2, the unique equilibrium E* is asymptotically stable for γ>γ0 and unstable for 0<γ<γ0, and system (7) undergoes Hopf bifurcation at γ=γ0.(iii)
In the regionD3, the positive equilibrium E2* is unstable for any γ>0; other two positive equilibria E1* and E3* are asymptotically stable for γ>γ0 and unstable for 0<γ<γ0, the and system (7) undergoes Hopf bifurcation near neighborhood of each of these two equilibria at γ=γ0.(iv)
In the regionD4, the positive equilibrium E2* is unstable for any γ>0; the positive equilibrium E1* is asymptotically stable for γ>γ0 and unstable for 0<γ<γ0, and system (7) undergoes Hopf bifurcation near neighborhood of E1* at γ=γ0; the positive equilibrium E3* is asymptotically stable for any γ>0.


## 2.2.   Direction and Stability of Hopf Bifurcation Associated withγ0

It follows from Theorem3 that when (α,β)∈Di,  i=2,3,4, Hopf bifurcation occurs at γ=γ0. When Hopf bifurcation occurs, the corresponding characteristic equation has a pair of purely imaginary roots at γ=γ0 and thus the linearization is not sufficient to determine the stability of the corresponding equilibrium, say, E*. In this subsection, by calculating the focal value we determine the stability of the equilibrium E* at γ=γ0, which is directly associated with the direction and stability of Hopf bifurcation associated with γ0.Assume thatλ1,2=η(γ)±iϖ(γ) are a pair of complex conjugates of (15) such that η(γ0)=0,  ϖ(γ0)=γ0(r1+r2β). Then the eigenvector corresponding to λ1,2 is
(22)(1-r1+η(γ)±iϖ(γ)r2).
Under the following coordinate transformation:
(23)(z1z2)=(10MN)(uv),
where M=-(r1+η(γ))/r2, N=ϖ(γ)/r2, system (7) becomes
(24)(u˙(t)v˙(t))=J(γ)(u(t)v(t))+(F(1)(u,v,γ)F(2)(u,v,γ)),
where
(25)J(γ)=(η(γ)-ϖ(γ)ϖ(γ)η(γ))
with
(26)F(1)(u,v,γ)=(12f20(1)+Mf11(1))u2+Nf11(1)uv+(16f30(1)+12Mf21(1))u3+12Nf21(1)u2v+O(|u,v|4),F(2)(u,v,γ)=-MNF(1)(u,v,γ)+1NF^(2)(u,v,γ),
where
(27)F^(2)(u,v,γ)=(12f20(2)+Mf11(2)+12M2f02(2))u2+(Nf11(2)+MNf02(2))uv+12N2f02(2)v2+(16f30(2)+12Mf21(2)+12M2f12(2))u3+(12Nf21(2)+MNf12(2))u2v+12N2f12(2)uv2+O(|u,v|4).
Find fij(1) and fij(2) at (31). By the formula for the third focal value of a multiple focus in [31], we have
(28)a3(γ0)=116(Fuuu(1)+Fuvv(1)+Fuuv(2)+Fvvv(2))|(0,0,γ0)+116ϖ(γ0)(Fuv(1)(Fuu(1)+Fvv(1))-Fuv(2)(Fuu(2)+Fvv(2))-Fuu(1)Fuu(2)+Fvv(1)Fvv(2))|(0,0,γ0),
where
(29)Fuuu(1)=f30(1)+3Mf21(1),Fuvv(1)=0,Fuuv(2)=-Mf21(1)+f21(2)+2Mf12(2),Fvvv(2)=0,Fuv(1)=Nf11(1),Fuu(1)=f20(1)+2Mf11(1),Fvv(1)=0,Fuv(2)=-Mf11(1)+f11(2)+Mf02(2),Fvv(2)=N2f02(2),Fuu(2)=-MN(12f20(1)+Mf11(1))+1N(12f20(2)+Mf11(2)+12M2f02(2)).It is well known from the result in [31] that if a3(γ0)<0, then the equilibrium E* of system (7) is a stable multiple focus of multiplicity 1 and the corresponding Hopf bifurcation is supercritical, and if a3(γ0)>0, then the equilibrium E* of system (7) is an unstable multiple focus of multiplicity 1 and the corresponding Hopf bifurcation is subcritical.


## 3. Stability and Hopf Bifurcation Analysis for the System with Delay

In this section, we focus on the investigation of the local stability and Hopf bifurcation criteria of the positive equilibriumE*(x*,y*) for the system (6). Letting z1(t)=x(t)-x* and z2(t)=y(t)-y*, we can rewrite the system (6) by Taylor series expression about (x*,y*) as the following system:
(30)z1′(t)=-r1z1(t)-r2z2(t-τ)+∑i+j≥2‍1i!j!fij(1)z1i(t)z2j(t-τ),z2′(t)=βγz1(t)-γz2(t)+∑i+j≥2‍1i!j!fij(2)z1i(t)z2j(t),
where
(31)r1=-x*y1′(x*)(x*)2+α,r2=x*(x*)2+α>0,fij(1)=∂i+jf(1)∂xi∂yj|(x*,y*),fij(2)=∂i+jf(2)∂xi∂yj|(x*,y*),f(1)=x(1-x)-xyx2+α,f(2)=γ(1-yβx)y.


### 3.1. Linear Stability and Delay-Induced Hopf Bifurcation

The characteristic equation atE* takes the form
(32)λ2+(r1+γ)λ+γr1+γr2βe-λτ=0.
When the time delay τ=0, the equilibrium E(x*,y*) is asymptotically stable if
(33)r1+γ>0,r1+r2β>0.
Now for τ>0, assuming that λ=iω(ω>0) is a root of (32), then we have
(34)-ω2+(r1+γ)iω+γr1+γr2βe-iωτ=0.
Separating the real and imaginary parts, we obtain
(35)ω2-γr1=γr2βcos(ωτ),(r1+γ)ω=γr2βsin(ωτ),
which leads to the following fourth degree polynomial equation:
(36)ω4+(r12+γ2)ω2+γ2r12-γ2r22β2=0.
It is easy to see that if
(37)r12-r22β2<0,
then (36) has only one positive root
(38)ω0=-(r12+γ2)+(r12-γ2)2+4γ2r22β22.
Substituting the value of ω0 in (35) and solve for τ, we get
(39)τj=1ω0arccos(ω02-γr1γr2β)+2jπω0,j=0,1,…
then when τ=τj, the characteristic equation (32) has a pair of purely imaginary roots ±iω0.Letλ(τ)=ζ(τ)+iω(τ) be a root of (32) such that the following two conditions hold:
(40)ζ(τj)=0,ω(τj)=ω0.
Substituting λ(τ)=ζ(τ)+iω(τ) into (32) and taking the derivative of resulting expression with respect to τ, we get
(41)dλdτ=λγr2βe-λτ2λ+r1+γ-τγr2βe-λτ
which, together with (35) and (38), leads to
(42)sign{Re[dλdτ]-1}λ=iω0=sign{Re[(2λ+r1+γ)eλτλγr2β]}λ=iω0=sign{2ω02+r12+γ2γ2r22β2}>0.Therefore, the transversality condition holds, and hence Hopf-bifurcation occurs atτ=τj. Because the quartic equation (36) has only one positive root, system (6) has no switching stability. Delay-induced small-amplitude period solution bifurcates from interior equilibrium point when τ passes through its critical magnitude τ=τ0. Here, τ0 is the smallest positive value of τj given in (39).To study the effect of the delay on the stability of the interior equilibrium in regionsDi,  i=1,2,3,4, we just need decide the sign of r1-r2β by (18), (20), and (37). Notice that r1-r2β=-(x*/((x*)2+α))(β+y1′(x*)). So we first determine at which point the slope of curve y1 is negative of the slope of line y2; that is, we need to find the solution of the following equations:
(43)(1-x)(α+x2)=βx,-3x2+2x-α+β=0,
which yields
(44)4x3-3x2+2αx-α=0.
Let q(x)=4x3-3x2+2αx-α. When α≥3/8, q′(x)=12(x-1/4)2-3/4+2α≥0, which implies that the function q(x) is increasing for any x>0. Notice that q(0)=-α<0 and q(1)=1+α>0, and we can deduce that (44) has only one positive real root for α≥3/8. When 0<α<3/8, q′(x)=12(x-(3+9-24α)/12)(x-(3-9-24α)/12) which implies that the function q(x) has a maximum value at x1=(3-9-24α)/12 and a minimal value at x2=(3+9-24α)/12. Furthermore, we can deduce that q(x1)=-1/8+9-24α/24-α/2-α9-24α/9<0 and q(x2)=-1/8-9-24α/24-α/2+α9-24α/9<0. So (44) has also only one positive real root for 0<α<3/8. By direct computation, the positive real root of (44) is
(45)T3=14+112(27+108α+1296α3-27α2+81α)1/3-2α-3/4(27+108α+1296α3-27α2+81α)1/3
for any α>0. Then β3=3T32-2T3+α divides the regions D1 and D4 into two parts, respectively; see Figure 3. When β>β3,  β+y1′(x*)>0, and when β<β3, β+y1′(x*)<0. Therefore,
(46)r1-r2β{<0,β>β3,=0,β=β3,>0,β<β3.Bifurcation diagram of positive equilibrium in theα-β plane. β=β3 divides D1 and D4 into D11,  D12 and D41,  D42, respectively.
(a)
(b)The above results are summarized in the following theorems.Theorem 4.
(i) In the regionD11, for any γ>0, the unique equilibrium E* is asymptotically stable for τ∈[0,τ0) and unstable for τ>τ0, and system (6) undergoes Hopf bifurcation near E* for τ=τj, j=0,1,….
(ii) In the regionD12, for any γ>0, the unique equilibrium E* is asymptotically stable for any τ≥0.Theorem 5.
In the regionD2,(i)
whenγ>γ0, the unique equilibrium E* is asymptotically stable for τ∈[0,τ0) and unstable for τ>τ0, and system (6) undergoes Hopf bifurcation near E* for τ=τj, j=0,1,…;(ii)
when0<γ<γ0, the unique equilibrium E* is unstable for any τ≥0, and system (6) undergoes Hopf bifurcation near E* for τ=τj, j=0,1,….Theorem 6.
In the regionD3,(i)
for anyγ>0, the positive equilibrium E2* is unstable for any τ≥0;(ii)
for other two positive equilibriaE1* and E3*, (a) when γ>γ0, they are asymptotically stable for τ∈[0,τ0) and unstable for τ>τ0, and system (6) undergoes Hopf bifurcation near each of these two equilibria for τ=τj, j=0,1,…; (b) when 0<γ<γ0, they are unstable for any τ≥0, and system (6) undergoes Hopf bifurcation near each of these two equilibria for τ=τj, j=0,1,….Theorem 7.
(i) In the regionD41, for any γ>0, the positive equilibrium E2* is unstable for any τ≥0; for any γ>0, the positive equilibrium E3* is asymptotically stable for τ∈[0,τ0) and unstable for τ>τ0, and system (6) undergoes Hopf bifurcation near E3* for τ=τj, j=0,1,…; for the positive equilibrium E1*, (a) when γ>γ0, E1* is asymptotically stable for τ∈[0,τ0) and unstable for τ>τ0, and system (6) undergoes Hopf bifurcation near E1* for τ=τj, j=0,1,…; (b) when 0<γ<γ0, E1* is unstable for any τ≥0, and system (6) undergoes Hopf bifurcation near E1* for τ=τj, j=0,1,….
(ii) In the regionD42, for any γ>0, the positive equilibrium E2* is unstable for any τ≥0, and system (6) undergoes Hopf bifurcation near E2* for τ=τj, j=0,1,…; the positive equilibrium E1* is asymptotically stable for γ>γ0 and unstable for 0<γ<γ0 for any τ≥0; the positive equilibrium E3* is asymptotically stable for any γ>0 and for any τ≥0.


### 3.2. Global Stability of the Positive Equilibrium

In this section, we give a result on the global stability of the positive equilibrium of system (6).Theorem 8.
Ifβ<α, then the unique positive equilibrium E*(x*,y*) of system (6) is global stability. That is, for any initial values x(t)=φ1(t)≥0,  y(t)=φ2(t)≥0,  t∈[-τ,0], with x(0),y(0)>0, the corresponding solution (x(t),y(t)) of system (6) satisfies
(47)limt→∞x(t)=x*,limt→∞y(t)=y*.Proof.
By the first equation of system (6), we have x˙(t)≤x(t)(1-x(t)), which, together with the comparison principle, implies that limsupt→+∞x(t)≤1. Thus, for any sufficiently small positive number ɛ0, there exists t1>0 such that for any t>t1,
(48)x(t)≤1+ɛ0≜c¯1.
This, together with the second equation of system (6) yields that for t>t1,
(49)y˙(t)≤γ(1-y(t)β(1+ɛ0))y(t).
Again by the comparison principle, we have
(50)liminft→+∞y(t)≤β(1+ɛ0),
which implies that for ɛ0>0 mentioned above, there exists t2>t1>0 such that for any t>t2,
(51)y(t)≤β(1+ɛ0)+ɛ0≜c¯2.
It follows from (51) and the second equation of system (6) that for t>t2+τ,
(52)x˙(t)≥x(t)(1-x(t))-β(1+ɛ0)+ɛ0αx(t),
which leads to
(53)liminft→+∞x(t)≥α-(β(1+ɛ0)+ɛ0)α.
So, there exists t3>t2>0 such that for any t>t3,
(54)x(t)≥α-(β(1+ϵ0)+ɛ0)α-ɛ0≜c_1.
Using this lower bound of x, we obtain
(55)y˙(t)≥γ(1-αβ(α(1-ɛ0)-(β(1+ɛ0)+ɛ0))y(t))y(t),
which, together with the comparison principle, yields
(56)liminft→+∞y(t)≥β(α(1-ɛ0)-(β(1+ϵ0)+ɛ0))α.
Therefore, there exists t4>t3>0 such that for any t>t4,
(57)y(t)≥β(α(1-ɛ0)-(β(1+ϵ0)+ɛ0))α-ɛ0≜c_2.
Clearly, c¯1,  c¯2>0. It follows from (54) and (57) that when α>β, we can choose
(58)0<ɛ0<β(α-β)α+β(α+β+1)
such that c_1,  c_2>0. By (48)–(57), it is easy to verify that
(59)0≥c¯1(1-c¯1)-c¯1c_2c¯12+α,0≥γ(1-c¯2βc¯1)c¯2,0≤c_1(1-c_1)-c_1c¯2c_12+α,0≤γ(1-c_2βc_1)c_2,
which means that (c¯1,c¯2) and (c_1,c_2) are coupled upper and lower solutions of system (6) (see Definition 2.2 in [32]).
Define two sequences of constant vectors(c¯1(m),c¯2(m)) and (c_1(m),c_2(m)),  m=1,2,…, as follows:
(60)c¯1(m)=c¯1(m-1)+1Lc¯1(m-1)(1-c¯1(m-1)-c_2(m-1)(c¯1(m-1))2+α),c¯2(m)=c¯2(m-1)+γLc¯2(m-1)(1-c¯2(m-1)βc¯1(m-1)),c_1(m)=c_1(m-1)+1Lc_1(m-1)(1-c_1(m-1)-c¯2(m-1)(c_1(m-1))2+α),c_2(m)=c_2(m-1)+γLc_2(m-1)(1-c_2(m-1)βc_1(m-1)),
where L is the Lipschitz constant for the vector field of system (6) when (c_1,c_2)≤(x,y)≤(c¯1,c¯2).
Notice that the above two sequences of constant vectors are recursive sequences starting with(c¯1(0),c¯2(0))=(c¯1,c¯2) and (c_1(0),c_2(0))=(c_1,c_2). Then, by the principle of induction, the following monotone property of these sequences holds (see Lemma 2.1 in [33] for more details):
(61)(0,0)<(c_1,c_2)≤(c_1(m),c_2(m))≤(c_1(m+1),c_2(m+1))≤⋯≤(c¯1(m+1),c¯2(m+1))≤(c¯1(m),c¯2(m))≤⋯≤(c¯1,c¯2).
This monotone property implies that there exist positive constants c¯i*>0 and c_i*>0 such that
(62)c_i≤c_i(m)≤c_i(m+1)≤c_i*≤c¯i*≤c¯i(m+1)≤c¯i(m)≤c¯i,(63)limm→∞c¯i(m)=c¯i*,limm→∞c_i(m)=c_i*,i=1,2.
By (60) and (63), we have
(64)c¯1*(1-c¯1*-c_2*(c¯1*)2+α)=0,γc¯2*(1-c¯2*βc¯1*)=0,c_1*(1-c_1*-c¯2*(c_1*)2+α)=0,γc_2*(1-c_2*βc_1*)=0.
From (64), it is easy to verify that
(65)(1-c¯1*)((c¯1*)2+α)=βc_1*,(1-c_1*)((c_1*)2+α)=βc¯1*.
Subtraction of these two equations gives
(66)(c¯1*)2-(c_1*)2+(c_1*)3-(c¯1*)3+α(c_1*-c¯1*)=β(c_1*-c¯1*).
Then we can prove c¯1*=c_1* by contradiction. In fact, if c¯1*≠c_1*, then (66) becomes
(67)-(c¯1*+c_1*)+(c_1*)2+(c¯1*)2+c¯1*c_1*+α=β.
Equation (65) can be written as follows
(68)(1-c¯1*)((c¯1*)2+α)+(1-c_1*)β=β,(1-c¯1*)β+(1-c_1*)((c_1*)2+α)=β.
Solving (68) gives
(69)1-c¯1*=β((c_1*)2+α-β)((c¯1*)2+α)((c_1*)2+α)-β2,1-c_1*=β((c¯1*)2+α-β)((c¯1*)2+α)((c_1*)2+α)-β2.
Substraction of the two equations of (69) and simplification lead to
(70)((c¯1*)2+α)((c_1*)2+α)-β2=β(c¯1*+c_1*).
Addition of the two equations of (69) and using (70) yield
(71)-(c¯1*+c_1*)+(c_1*)2+(c¯1*)2+c¯1*c_1*=2(α-β).
It follows from (67) and (71) that
(72)3(α-β)=0,
which is a contradiction under the assumption α>β. So, the statement c¯1*=c_1* is proved.
By (64), we also have c¯2*=βc¯1*,  c_2*=βc_1*. Therefore, c¯2*=c_2* since c¯1*=c_1*. Letting c1*=c¯1*=c_1* and c2*=c¯2*=c_2*, then (c1*,c2*) is the positive solution of (8). In addition, notice that when α>β, (8) has the unique positive solution (x*,y*). Thus, (c1*,c2*)=(x*,y*). Then from the result due to Pao (see Theorem 2.2 in [33]), the proof of Theorem is complete.


### 3.3. Direction and Stability of Hopf Bifurcation Associated with the Critical Values of Delay

In Section3.1, we have shown that the system (6) undergoes the Hopf bifurcation at critical values τ=τj. In this subsection, we will derive the direction of the Hopf bifurcation and stability of the bifurcating periodic orbits from the interior equilibrium E* of system (6) at τ=τj. We will employ the algorithm of Faria and Magalhães [34] to compute explicitly the normal forms of system (6) on the center manifold.We rescale the time byt↦t/τ to normalize the delay and then introduce the new parameter μ=τ-τj so that μ=0 is the critical value of Hopf singularity. System (30) can be rewritten as
(73)z˙(t)=L(μ)zt+F(zt,μ),
in the phase space C=C([-1,0],ℝ2), where for φ=(φ1,φ2)T∈C, we have
(74)L(μ)(φ)=(τj+μ)(-r1φ1(0)-r2φ2(-1)γβφ1(0)-γφ2(0)),F(φ,μ)=(τj+μ)(∑i+j≥2‍1i!j!fij(1)φ1i(0)φ2j(-1)∑i+j≥2‍1i!j!fij(2)φ1i(0)φ2j(0)).
By the Riesz representation theorem, there exists a function η:[-1,0]→ℝ2×ℝ2 of bounded variation such that
(75)L(μ)φ=∫-10‍[dημ(θ)]φ(θ).Settingω*=ω0τj, then Λ0={iω*,-iω*} and B=diag{iω*,-iω*}. Let P be the generalized eigenspace associated with the eigenvalues in Λ0 and let P* be its dual space. For ℝ2*, the 2-dimensional space of row vectors, define C*=C([0,1];ℝ2*) and the adjoint bilinear form on C*×C as follows:
(76)〈ψ(s),ϕ(θ)〉=ψ(0)ϕ(0)-∫-10‍∫0θ‍ψ(ξ-θ)dημ(θ)ϕ(ξ)dξ,ψ(s)∈C*,ϕ(θ)∈C.
Denote the dual bases of P and P* by Φ and Ψ, respectively. Consider Φ=(ϕ1,ϕ2) with
(77)ϕ1(θ)=eiω*θv,ϕ2(θ)=e-iω*θv¯,
where v is a nonzero solution in ℂ2 to the following equation:
(78)L0(eiω*·I)v=iω*v,
and choose a basis Ψ for the adjoint space P* such that 〈Ψ,Φ〉=I2, where I2 is a 2×2 identical matrix. So Ψ(s)=col(ψ1(s),ψ2(s)) with
(79)ψ1(s)=e-iω*suT,ψ2(s)=eiω*su¯T,0≤s≤1,
where uT∈ℂ2 is the solution to the equation
(80)uTL0(eiω*·I)=iω*uT.
and satisfies 〈ψ1(s),ϕ1(θ)〉=1. By direct computation, we can choose
(81)v=(v1v2)=(1-r1+iω0r2e-iω*),u=(u1u2)=d(-γ+iω0r2e-iω*1),
where d=-r2e-iω*/(r1+γ+2iω0+r1γτj+(r1+γ)iω*-ω0ω*).Following the normal form theory of functional differential equations due to Faria and Magalhães [34] and using a very similar procedure as in [29] but a tedious calculation, we can obtain the following normal form at the critical values τj:
(82)w˙=Bw+(A1w1μA¯1w2μ)+(A2w12w2A¯2w1w22)+O(|w|μ2+|w4|),
where w=(w1,w2)T∈ℂ2. The coefficients A1 and A2 can be determined as follows:
(83)A1=iω0uTv,A2=i2ω*(a20a11-2|a11|2-13|a02|2)+12(a21+b21),
where
(84)a20=τj(b1u1+b2u2),a11=τj(b3u1+b4u2),a02=τj(b¯1u1+b¯2u2)
with
(85)b1=f20(1)v12+2f11(1)v1v2e-iω*,b2=f20(2)v12+2f11(2)v1v2+f02(2)v22,b3=f20(1)|v1|2+2f11(1)Re{v1v¯2eiω*},b4=f20(2)|v1|2+2f11(2)Re{v1v¯2}+f02(2)|v2|2,a21=τju1(f30(1)v1|v1|2+f21(1)(v12v¯2eiω*+2|v1|2v2e-iω*))+τju2(f30(2)v1|v1|2+f21(2)(v12v¯2+2|v1|2v2)+f12(2)(v22v¯1+2|v2|2v1)),b21=τjuT×(d1h11(1)(0)+d¯1h20(1)(0)+d2h11(2)(-1)+d¯2h20(2)(-1)d3h11(1)(0)+d¯3h20(1)(0)+d4h11(2)(0)+d¯4h20(2)(0))
with
(86)d1=f20(1)v1+f11(1)v2e-iω*,d2=f11(1)v1,d3=f20(2)v1+f11(2)v2,d4=f11(2)v1+f02(2)v2,h20(θ)=-1iω*(a20eiω*θv+13a¯02e-iω*θv¯)+e2iω*θ(b1(2iω0+γ)-b2r2e-2iω*(2iω0+γ)(2iω0+r1)+βr2γe-2iω*βγb1+b2(2iω0+r1)(2iω0+γ)(2iω0+r1)+βr2γe-2iω*),h11(θ)=2iω*(a11eiω*θv-a¯11e-iω*θv¯)+(2γb3-2r2b4γ(r1+r2β)2βγb3+2r1b4γ(r1+r2β)).
The normal form (82) can be written in real coordinates w~=(w~1,w~2)2∈ℝ2 through the change of variables w~1=w1-iw2, w~2=w1+iw2. Transforming to polar coordinates w~1=ρcosξ, w~2=ρsinξ, the normal form becomes
(87)ρ˙=K1μρ+K2ρ3+O(μ2ρ+|(ρ,x)|4),ξ˙=-ω*+O(|(ρ,μ)|)
with K1=ReA1, K2=ReA2. It is well known [35] that the sign of K1K2 determines the direction of the bifurcation (supercritical if K1K2<0, subcritical if K1K2>0), and the sign of K2 determines the stability of the nontrivial periodic orbits (stable if K2<0, unstable if K2>0). Thus, if the coefficients of system (4) are given, then we will analyze the direction of the Hopf bifurcation and stability of the bifurcating periodic orbits at τ=τj.


## 3.1. Linear Stability and Delay-Induced Hopf Bifurcation

The characteristic equation atE* takes the form
(32)λ2+(r1+γ)λ+γr1+γr2βe-λτ=0.
When the time delay τ=0, the equilibrium E(x*,y*) is asymptotically stable if
(33)r1+γ>0,r1+r2β>0.
Now for τ>0, assuming that λ=iω(ω>0) is a root of (32), then we have
(34)-ω2+(r1+γ)iω+γr1+γr2βe-iωτ=0.
Separating the real and imaginary parts, we obtain
(35)ω2-γr1=γr2βcos(ωτ),(r1+γ)ω=γr2βsin(ωτ),
which leads to the following fourth degree polynomial equation:
(36)ω4+(r12+γ2)ω2+γ2r12-γ2r22β2=0.
It is easy to see that if
(37)r12-r22β2<0,
then (36) has only one positive root
(38)ω0=-(r12+γ2)+(r12-γ2)2+4γ2r22β22.
Substituting the value of ω0 in (35) and solve for τ, we get
(39)τj=1ω0arccos(ω02-γr1γr2β)+2jπω0,j=0,1,…
then when τ=τj, the characteristic equation (32) has a pair of purely imaginary roots ±iω0.Letλ(τ)=ζ(τ)+iω(τ) be a root of (32) such that the following two conditions hold:
(40)ζ(τj)=0,ω(τj)=ω0.
Substituting λ(τ)=ζ(τ)+iω(τ) into (32) and taking the derivative of resulting expression with respect to τ, we get
(41)dλdτ=λγr2βe-λτ2λ+r1+γ-τγr2βe-λτ
which, together with (35) and (38), leads to
(42)sign{Re[dλdτ]-1}λ=iω0=sign{Re[(2λ+r1+γ)eλτλγr2β]}λ=iω0=sign{2ω02+r12+γ2γ2r22β2}>0.Therefore, the transversality condition holds, and hence Hopf-bifurcation occurs atτ=τj. Because the quartic equation (36) has only one positive root, system (6) has no switching stability. Delay-induced small-amplitude period solution bifurcates from interior equilibrium point when τ passes through its critical magnitude τ=τ0. Here, τ0 is the smallest positive value of τj given in (39).To study the effect of the delay on the stability of the interior equilibrium in regionsDi,  i=1,2,3,4, we just need decide the sign of r1-r2β by (18), (20), and (37). Notice that r1-r2β=-(x*/((x*)2+α))(β+y1′(x*)). So we first determine at which point the slope of curve y1 is negative of the slope of line y2; that is, we need to find the solution of the following equations:
(43)(1-x)(α+x2)=βx,-3x2+2x-α+β=0,
which yields
(44)4x3-3x2+2αx-α=0.
Let q(x)=4x3-3x2+2αx-α. When α≥3/8, q′(x)=12(x-1/4)2-3/4+2α≥0, which implies that the function q(x) is increasing for any x>0. Notice that q(0)=-α<0 and q(1)=1+α>0, and we can deduce that (44) has only one positive real root for α≥3/8. When 0<α<3/8, q′(x)=12(x-(3+9-24α)/12)(x-(3-9-24α)/12) which implies that the function q(x) has a maximum value at x1=(3-9-24α)/12 and a minimal value at x2=(3+9-24α)/12. Furthermore, we can deduce that q(x1)=-1/8+9-24α/24-α/2-α9-24α/9<0 and q(x2)=-1/8-9-24α/24-α/2+α9-24α/9<0. So (44) has also only one positive real root for 0<α<3/8. By direct computation, the positive real root of (44) is
(45)T3=14+112(27+108α+1296α3-27α2+81α)1/3-2α-3/4(27+108α+1296α3-27α2+81α)1/3
for any α>0. Then β3=3T32-2T3+α divides the regions D1 and D4 into two parts, respectively; see Figure 3. When β>β3,  β+y1′(x*)>0, and when β<β3, β+y1′(x*)<0. Therefore,
(46)r1-r2β{<0,β>β3,=0,β=β3,>0,β<β3.Bifurcation diagram of positive equilibrium in theα-β plane. β=β3 divides D1 and D4 into D11,  D12 and D41,  D42, respectively.
(a)
(b)The above results are summarized in the following theorems.Theorem 4.
(i) In the regionD11, for any γ>0, the unique equilibrium E* is asymptotically stable for τ∈[0,τ0) and unstable for τ>τ0, and system (6) undergoes Hopf bifurcation near E* for τ=τj, j=0,1,….
(ii) In the regionD12, for any γ>0, the unique equilibrium E* is asymptotically stable for any τ≥0.Theorem 5.
In the regionD2,(i)
whenγ>γ0, the unique equilibrium E* is asymptotically stable for τ∈[0,τ0) and unstable for τ>τ0, and system (6) undergoes Hopf bifurcation near E* for τ=τj, j=0,1,…;(ii)
when0<γ<γ0, the unique equilibrium E* is unstable for any τ≥0, and system (6) undergoes Hopf bifurcation near E* for τ=τj, j=0,1,….Theorem 6.
In the regionD3,(i)
for anyγ>0, the positive equilibrium E2* is unstable for any τ≥0;(ii)
for other two positive equilibriaE1* and E3*, (a) when γ>γ0, they are asymptotically stable for τ∈[0,τ0) and unstable for τ>τ0, and system (6) undergoes Hopf bifurcation near each of these two equilibria for τ=τj, j=0,1,…; (b) when 0<γ<γ0, they are unstable for any τ≥0, and system (6) undergoes Hopf bifurcation near each of these two equilibria for τ=τj, j=0,1,….Theorem 7.
(i) In the regionD41, for any γ>0, the positive equilibrium E2* is unstable for any τ≥0; for any γ>0, the positive equilibrium E3* is asymptotically stable for τ∈[0,τ0) and unstable for τ>τ0, and system (6) undergoes Hopf bifurcation near E3* for τ=τj, j=0,1,…; for the positive equilibrium E1*, (a) when γ>γ0, E1* is asymptotically stable for τ∈[0,τ0) and unstable for τ>τ0, and system (6) undergoes Hopf bifurcation near E1* for τ=τj, j=0,1,…; (b) when 0<γ<γ0, E1* is unstable for any τ≥0, and system (6) undergoes Hopf bifurcation near E1* for τ=τj, j=0,1,….
(ii) In the regionD42, for any γ>0, the positive equilibrium E2* is unstable for any τ≥0, and system (6) undergoes Hopf bifurcation near E2* for τ=τj, j=0,1,…; the positive equilibrium E1* is asymptotically stable for γ>γ0 and unstable for 0<γ<γ0 for any τ≥0; the positive equilibrium E3* is asymptotically stable for any γ>0 and for any τ≥0.


## 3.2. Global Stability of the Positive Equilibrium

In this section, we give a result on the global stability of the positive equilibrium of system (6).Theorem 8.
Ifβ<α, then the unique positive equilibrium E*(x*,y*) of system (6) is global stability. That is, for any initial values x(t)=φ1(t)≥0,  y(t)=φ2(t)≥0,  t∈[-τ,0], with x(0),y(0)>0, the corresponding solution (x(t),y(t)) of system (6) satisfies
(47)limt→∞x(t)=x*,limt→∞y(t)=y*.Proof.
By the first equation of system (6), we have x˙(t)≤x(t)(1-x(t)), which, together with the comparison principle, implies that limsupt→+∞x(t)≤1. Thus, for any sufficiently small positive number ɛ0, there exists t1>0 such that for any t>t1,
(48)x(t)≤1+ɛ0≜c¯1.
This, together with the second equation of system (6) yields that for t>t1,
(49)y˙(t)≤γ(1-y(t)β(1+ɛ0))y(t).
Again by the comparison principle, we have
(50)liminft→+∞y(t)≤β(1+ɛ0),
which implies that for ɛ0>0 mentioned above, there exists t2>t1>0 such that for any t>t2,
(51)y(t)≤β(1+ɛ0)+ɛ0≜c¯2.
It follows from (51) and the second equation of system (6) that for t>t2+τ,
(52)x˙(t)≥x(t)(1-x(t))-β(1+ɛ0)+ɛ0αx(t),
which leads to
(53)liminft→+∞x(t)≥α-(β(1+ɛ0)+ɛ0)α.
So, there exists t3>t2>0 such that for any t>t3,
(54)x(t)≥α-(β(1+ϵ0)+ɛ0)α-ɛ0≜c_1.
Using this lower bound of x, we obtain
(55)y˙(t)≥γ(1-αβ(α(1-ɛ0)-(β(1+ɛ0)+ɛ0))y(t))y(t),
which, together with the comparison principle, yields
(56)liminft→+∞y(t)≥β(α(1-ɛ0)-(β(1+ϵ0)+ɛ0))α.
Therefore, there exists t4>t3>0 such that for any t>t4,
(57)y(t)≥β(α(1-ɛ0)-(β(1+ϵ0)+ɛ0))α-ɛ0≜c_2.
Clearly, c¯1,  c¯2>0. It follows from (54) and (57) that when α>β, we can choose
(58)0<ɛ0<β(α-β)α+β(α+β+1)
such that c_1,  c_2>0. By (48)–(57), it is easy to verify that
(59)0≥c¯1(1-c¯1)-c¯1c_2c¯12+α,0≥γ(1-c¯2βc¯1)c¯2,0≤c_1(1-c_1)-c_1c¯2c_12+α,0≤γ(1-c_2βc_1)c_2,
which means that (c¯1,c¯2) and (c_1,c_2) are coupled upper and lower solutions of system (6) (see Definition 2.2 in [32]).
Define two sequences of constant vectors(c¯1(m),c¯2(m)) and (c_1(m),c_2(m)),  m=1,2,…, as follows:
(60)c¯1(m)=c¯1(m-1)+1Lc¯1(m-1)(1-c¯1(m-1)-c_2(m-1)(c¯1(m-1))2+α),c¯2(m)=c¯2(m-1)+γLc¯2(m-1)(1-c¯2(m-1)βc¯1(m-1)),c_1(m)=c_1(m-1)+1Lc_1(m-1)(1-c_1(m-1)-c¯2(m-1)(c_1(m-1))2+α),c_2(m)=c_2(m-1)+γLc_2(m-1)(1-c_2(m-1)βc_1(m-1)),
where L is the Lipschitz constant for the vector field of system (6) when (c_1,c_2)≤(x,y)≤(c¯1,c¯2).
Notice that the above two sequences of constant vectors are recursive sequences starting with(c¯1(0),c¯2(0))=(c¯1,c¯2) and (c_1(0),c_2(0))=(c_1,c_2). Then, by the principle of induction, the following monotone property of these sequences holds (see Lemma 2.1 in [33] for more details):
(61)(0,0)<(c_1,c_2)≤(c_1(m),c_2(m))≤(c_1(m+1),c_2(m+1))≤⋯≤(c¯1(m+1),c¯2(m+1))≤(c¯1(m),c¯2(m))≤⋯≤(c¯1,c¯2).
This monotone property implies that there exist positive constants c¯i*>0 and c_i*>0 such that
(62)c_i≤c_i(m)≤c_i(m+1)≤c_i*≤c¯i*≤c¯i(m+1)≤c¯i(m)≤c¯i,(63)limm→∞c¯i(m)=c¯i*,limm→∞c_i(m)=c_i*,i=1,2.
By (60) and (63), we have
(64)c¯1*(1-c¯1*-c_2*(c¯1*)2+α)=0,γc¯2*(1-c¯2*βc¯1*)=0,c_1*(1-c_1*-c¯2*(c_1*)2+α)=0,γc_2*(1-c_2*βc_1*)=0.
From (64), it is easy to verify that
(65)(1-c¯1*)((c¯1*)2+α)=βc_1*,(1-c_1*)((c_1*)2+α)=βc¯1*.
Subtraction of these two equations gives
(66)(c¯1*)2-(c_1*)2+(c_1*)3-(c¯1*)3+α(c_1*-c¯1*)=β(c_1*-c¯1*).
Then we can prove c¯1*=c_1* by contradiction. In fact, if c¯1*≠c_1*, then (66) becomes
(67)-(c¯1*+c_1*)+(c_1*)2+(c¯1*)2+c¯1*c_1*+α=β.
Equation (65) can be written as follows
(68)(1-c¯1*)((c¯1*)2+α)+(1-c_1*)β=β,(1-c¯1*)β+(1-c_1*)((c_1*)2+α)=β.
Solving (68) gives
(69)1-c¯1*=β((c_1*)2+α-β)((c¯1*)2+α)((c_1*)2+α)-β2,1-c_1*=β((c¯1*)2+α-β)((c¯1*)2+α)((c_1*)2+α)-β2.
Substraction of the two equations of (69) and simplification lead to
(70)((c¯1*)2+α)((c_1*)2+α)-β2=β(c¯1*+c_1*).
Addition of the two equations of (69) and using (70) yield
(71)-(c¯1*+c_1*)+(c_1*)2+(c¯1*)2+c¯1*c_1*=2(α-β).
It follows from (67) and (71) that
(72)3(α-β)=0,
which is a contradiction under the assumption α>β. So, the statement c¯1*=c_1* is proved.
By (64), we also have c¯2*=βc¯1*,  c_2*=βc_1*. Therefore, c¯2*=c_2* since c¯1*=c_1*. Letting c1*=c¯1*=c_1* and c2*=c¯2*=c_2*, then (c1*,c2*) is the positive solution of (8). In addition, notice that when α>β, (8) has the unique positive solution (x*,y*). Thus, (c1*,c2*)=(x*,y*). Then from the result due to Pao (see Theorem 2.2 in [33]), the proof of Theorem is complete.


## 3.3. Direction and Stability of Hopf Bifurcation Associated with the Critical Values of Delay

In Section3.1, we have shown that the system (6) undergoes the Hopf bifurcation at critical values τ=τj. In this subsection, we will derive the direction of the Hopf bifurcation and stability of the bifurcating periodic orbits from the interior equilibrium E* of system (6) at τ=τj. We will employ the algorithm of Faria and Magalhães [34] to compute explicitly the normal forms of system (6) on the center manifold.We rescale the time byt↦t/τ to normalize the delay and then introduce the new parameter μ=τ-τj so that μ=0 is the critical value of Hopf singularity. System (30) can be rewritten as
(73)z˙(t)=L(μ)zt+F(zt,μ),
in the phase space C=C([-1,0],ℝ2), where for φ=(φ1,φ2)T∈C, we have
(74)L(μ)(φ)=(τj+μ)(-r1φ1(0)-r2φ2(-1)γβφ1(0)-γφ2(0)),F(φ,μ)=(τj+μ)(∑i+j≥2‍1i!j!fij(1)φ1i(0)φ2j(-1)∑i+j≥2‍1i!j!fij(2)φ1i(0)φ2j(0)).
By the Riesz representation theorem, there exists a function η:[-1,0]→ℝ2×ℝ2 of bounded variation such that
(75)L(μ)φ=∫-10‍[dημ(θ)]φ(θ).Settingω*=ω0τj, then Λ0={iω*,-iω*} and B=diag{iω*,-iω*}. Let P be the generalized eigenspace associated with the eigenvalues in Λ0 and let P* be its dual space. For ℝ2*, the 2-dimensional space of row vectors, define C*=C([0,1];ℝ2*) and the adjoint bilinear form on C*×C as follows:
(76)〈ψ(s),ϕ(θ)〉=ψ(0)ϕ(0)-∫-10‍∫0θ‍ψ(ξ-θ)dημ(θ)ϕ(ξ)dξ,ψ(s)∈C*,ϕ(θ)∈C.
Denote the dual bases of P and P* by Φ and Ψ, respectively. Consider Φ=(ϕ1,ϕ2) with
(77)ϕ1(θ)=eiω*θv,ϕ2(θ)=e-iω*θv¯,
where v is a nonzero solution in ℂ2 to the following equation:
(78)L0(eiω*·I)v=iω*v,
and choose a basis Ψ for the adjoint space P* such that 〈Ψ,Φ〉=I2, where I2 is a 2×2 identical matrix. So Ψ(s)=col(ψ1(s),ψ2(s)) with
(79)ψ1(s)=e-iω*suT,ψ2(s)=eiω*su¯T,0≤s≤1,
where uT∈ℂ2 is the solution to the equation
(80)uTL0(eiω*·I)=iω*uT.
and satisfies 〈ψ1(s),ϕ1(θ)〉=1. By direct computation, we can choose
(81)v=(v1v2)=(1-r1+iω0r2e-iω*),u=(u1u2)=d(-γ+iω0r2e-iω*1),
where d=-r2e-iω*/(r1+γ+2iω0+r1γτj+(r1+γ)iω*-ω0ω*).Following the normal form theory of functional differential equations due to Faria and Magalhães [34] and using a very similar procedure as in [29] but a tedious calculation, we can obtain the following normal form at the critical values τj:
(82)w˙=Bw+(A1w1μA¯1w2μ)+(A2w12w2A¯2w1w22)+O(|w|μ2+|w4|),
where w=(w1,w2)T∈ℂ2. The coefficients A1 and A2 can be determined as follows:
(83)A1=iω0uTv,A2=i2ω*(a20a11-2|a11|2-13|a02|2)+12(a21+b21),
where
(84)a20=τj(b1u1+b2u2),a11=τj(b3u1+b4u2),a02=τj(b¯1u1+b¯2u2)
with
(85)b1=f20(1)v12+2f11(1)v1v2e-iω*,b2=f20(2)v12+2f11(2)v1v2+f02(2)v22,b3=f20(1)|v1|2+2f11(1)Re{v1v¯2eiω*},b4=f20(2)|v1|2+2f11(2)Re{v1v¯2}+f02(2)|v2|2,a21=τju1(f30(1)v1|v1|2+f21(1)(v12v¯2eiω*+2|v1|2v2e-iω*))+τju2(f30(2)v1|v1|2+f21(2)(v12v¯2+2|v1|2v2)+f12(2)(v22v¯1+2|v2|2v1)),b21=τjuT×(d1h11(1)(0)+d¯1h20(1)(0)+d2h11(2)(-1)+d¯2h20(2)(-1)d3h11(1)(0)+d¯3h20(1)(0)+d4h11(2)(0)+d¯4h20(2)(0))
with
(86)d1=f20(1)v1+f11(1)v2e-iω*,d2=f11(1)v1,d3=f20(2)v1+f11(2)v2,d4=f11(2)v1+f02(2)v2,h20(θ)=-1iω*(a20eiω*θv+13a¯02e-iω*θv¯)+e2iω*θ(b1(2iω0+γ)-b2r2e-2iω*(2iω0+γ)(2iω0+r1)+βr2γe-2iω*βγb1+b2(2iω0+r1)(2iω0+γ)(2iω0+r1)+βr2γe-2iω*),h11(θ)=2iω*(a11eiω*θv-a¯11e-iω*θv¯)+(2γb3-2r2b4γ(r1+r2β)2βγb3+2r1b4γ(r1+r2β)).
The normal form (82) can be written in real coordinates w~=(w~1,w~2)2∈ℝ2 through the change of variables w~1=w1-iw2, w~2=w1+iw2. Transforming to polar coordinates w~1=ρcosξ, w~2=ρsinξ, the normal form becomes
(87)ρ˙=K1μρ+K2ρ3+O(μ2ρ+|(ρ,x)|4),ξ˙=-ω*+O(|(ρ,μ)|)
with K1=ReA1, K2=ReA2. It is well known [35] that the sign of K1K2 determines the direction of the bifurcation (supercritical if K1K2<0, subcritical if K1K2>0), and the sign of K2 determines the stability of the nontrivial periodic orbits (stable if K2<0, unstable if K2>0). Thus, if the coefficients of system (4) are given, then we will analyze the direction of the Hopf bifurcation and stability of the bifurcating periodic orbits at τ=τj.


## 4. Numerical Simulations

In this section, we present some numerical simulations to verify and extend our theoretical analysis proved in Section3 by using the Simulink for Matlab. To investigate the effect of delay, we are interested in the cases when (α,β)∈D11, (α,β)∈D2, and (α,β)∈D3. We would also like to mention that for the case of delay, the initial values should be given by x(t)=φ1(t),  y(t)=φ2(t) for t∈[-τ,0]. However, in the Simulink for Matlab, we only need to specify the initial values as x(0) and y(0). The initial values for the case of delay can be selected automatically by the default function of Simulink.(i) Takingα=0.2 and β=2 such that (α,β)∈D11, system (6) has a unique positive equilibrium E*(0.0946,0.1892). Without delay, the positive equilibrium E* is asymptotically stable for any γ>0. Let γ=0.2 and then by (39), we have the critical value τ0≐1.2508 of delay and it follows from Theorem 4 that E* is asymptotically stable for 0≤τ<τ0≐1.2508 and unstable for τ>τ0≐1.2508, as shown in Figures 4(a) and 4(b) for τ=0.5. Furthermore, by the procedure of Section 3.3, we get K1≐0.0882>0, and K2≐5.2428>0, which implies that the Hopf bifurcation associated with this critical value τ0≐1.2508 is subcritical. So, there exists a sufficiently small positive real number ɛ0 such that for each τ∈(τ0-ɛ0,τ0), system (6) has an unstable periodic solution near the positive equilibrium E*(0.0946,0.1892). Taking τ=1.15<τ0, Figures 4(c) and 4(d) show the existence of the subcritical Hopf bifurcation. In this case, a small-amplitude unstable periodic solution, coming from the Hopf bifurcation, and a large-amplitude stable periodic solution coexist. Taking τ=1.3>τ0, the positive equilibrium E* becomes unstable and numerical simulation also shows the existence of large-amplitude stable periodic solution (see Figures 4(e) and 4(f)).The trajectory graphs and phase portrait of system (6) for (α,β)∈D11 with α=0.2,  β=2, and γ=0.2. (a)-(b) When 0≤τ<τ0, the equilibrium E* is stable. (c)-(d) Simulations of subcritical Hopf bifurcation. Starting from the initial value (x(0),y(0))=(0.0615,0.2897), the trajectory slowly gets away from a small-amplitude unstable periodic solution and finally converges to a large-amplitude stable periodic solution. Starting from the initial value (x(0),y(0))=(0.0615,0.2895), the trajectory converges to E*. (e)-(f) When τ>τ0, the equilibrium E* is unstable with a large-amplitude stable periodic solution.
(a)
(b)
(c)
(d)
(e)
(f)(ii) Takingα=0.2 and β=1.2 such that (α,β)∈D2, system (6) has a unique positive equilibrium E*(0.1578,0.1894). Without delay, the parameter γ affects the stability of the positive equilibrium E*. By Lemma 2, we get γ0≐0.0287. The positive equilibrium E* is stable for γ>γ0, and unstable for γ<γ0 and system (7) undergoes Hopf bifurcation near the positive equilibrium E* at γ=γ0. By (28), we have a3(γ0)=0.3822>0, which means that the Hopf bifurcation occurring at γ0 is subcritical. So, there exists a sufficiently small positive real number ɛ1 such that for any γ∈(γ0,γ0+ɛ1), system (7) has an unstable periodic orbit near the positive equilibrium E*. Taking γ=0.02<γ0, Figures 5(a) and 5(b) illustrate the instability of the positive equilibrium E* and the existence of a large-amplitude stable periodic orbit. Taking γ=0.04 larger than and closer to γ0, Figures 5(c) and 5(d) illustrate the stability of the positive equilibrium E* and the existence of subcritical Hopf bifurcation. In this case, a small-amplitude unstable periodic solution, coming from the Hopf bifurcation, and a large-amplitude stable periodic solution coexist. Taking γ=0.06 larger than and away from γ0, the positive equilibrium E* is stable and there are no Hopf bifurcating periodic orbits, as shown in Figures 5(e) and 5(f).The trajectory graphs and phase portrait of system (7) without delay for (α,β)∈D2 with α=0.2 and β=1.2. (a)-(b) γ=0.02<γ0. The positive equilibrium E* is unstable with a large-amplitude stable periodic orbit. The initial value is (0.1578,0.18935). (c)-(d) γ=0.04>γ0. The positive equilibrium E* is stable, and a small-amplitude unstable periodic solution and a large-amplitude stable periodic solution coexist. The initial values are (0.1578,0.21532) and (0.1578,0.21533). (e)-(f) γ=0.06>γ0. The positive equilibrium E* is stable and there are no Hopf bifurcating periodic orbits. The initial value is (0.5,0.3).
(a)
(b)
(c)
(d)
(e)
(f)To investigate the effect of delay, we fixα,  β, and γ. For (α,β)=(0.2,1.2)∈D2 and γ=0.04∈(γ0,γ0+ɛ1), the positive equilibrium E*(0.1578,0.1894) is stable and there exists an unstable periodic orbit near the positive equilibrium E* (see Figures 5(c) and 5(d)). The stability of the positive equilibrium E* changes with the increase of delay. By (39), the critical value of delay is τ0≐0.3346. By the procedure of Section 3.3, we have K1=0.0056, and K2=1.1923, which implies that the Hopf bifurcation occurring at τ0≐0.3346 is subcritical. So, the positive equilibrium E* is stable with the coexistence of small-amplitude stable and large-amplitude stable periodic orbits for 0≤τ<τ0≐0.3346 and unstable for τ>τ0≐0.3346. Figures 6(a) and 6(b) and Figures 6(c) and 6(d) are numerical simulations for τ=0.2∈[0,τ0) and τ=0.4>τ0, respectively. For (α,β)=(0.2,1.2)∈D2 and γ=0.06>γ0+ɛ1, the positive equilibrium E*(0.1578,0.1894) is stable and there are no periodic orbits near the positive equilibrium E* (see Figures 5(e) and 5(f)). In this case, the critical value of delay is τ0≐0.6207. The positive equilibrium E* is stable for 0≤τ<τ0≐0.6207 and unstable for τ>τ0≐0.6207. Figures 7(a) and 7(b), Figures 7(c) and 7(d), and Figures 7(e) and 7(f) are numerical simulations for τ=0.04<τ0≐0.6207, τ=0.05<τ0≐0.6207, and τ=0.65>τ0≐0.6207, respectively. It follows from the procedure of Section 3.3 that K1=0.0154 and K2=2.1776, which implies that the Hopf bifurcation occurring at τ0≐0.6207 is also subcritical. Thus, there exists a positive real number ɛ2 such that for each τ∈(τ0-ɛ2,τ0), there exists an unstable periodic orbit, as shown in Figures 7(c) and 7(d). For τ>τ0≐0.6207, the the positive equilibrium E* becomes unstable and there exists a large-amplitude stable periodic orbit, as shown in Figures 7(e) and 7(f) for τ=0.65>τ0.The trajectory graphs and phase portrait of system (6) for (α,β)∈D2 with α=0.2,  β=1.2, and γ=0.04. (a)-(b) For τ=0.2∈[0,τ0), the positive equilibrium E* is stable with small-amplitude unstable and large-amplitude stable periodic orbits. The initial values are (0.1578,0.21532) and (0.1578,0.21533). (c)-(d) For τ=0.4>τ0, the positive equilibrium E* is unstable with a large-amplitude stable periodic orbit. The initial value is (0.1578,0.18935).
(a)
(b)
(c)
(d)The trajectory graphs and phase portrait of system (6) for (α,β)∈D2 with α=0.2,  β=1.2, and γ=0.06. (a)-(b) For τ=0.04<τ0, the positive equilibrium E* is stable. The initial values are (0.5,0.3). (c)-(d) For τ=0.05<τ0, the positive equilibrium E* is stable with an unstable periodic orbit. The initial values are (0.4,0.258) and (0.5,0.3). (e)-(f) For τ=0.65>τ0, the positive equilibrium E* is unstable with an stable periodic orbit. The initial values are (0.1578,0.18935) and (0.5,0.3).
(a)
(b)
(c)
(d)
(e)
(f)(iii) Takingα=0.01 and β=0.24 such that (α,β)∈D3, system (6) has three positive equilibria E1*(0.0492,0.0118),  E2*(0.3244,0.0779), and E3*(0.6264,0.1503). The positive equilibrium E2* is always unstable for any γ>0 and τ≥0. But the parameters γ and τ affect the stability of E1* and E3*. Without delay, the critical values of γ for E1* and E3* are γ01=0.3215 and γ03=0.1024. The positive equilibrium E1* (resp., E3*) is asymptotically stable for γ>γ01 (resp., γ>γ03) and unstable for γ<γ01 (resp., γ<γ03). Taking γ=0.1<min{γ01,γ03}, all of these three equilibria are unstable. There is a large-amplitude stable periodic orbit enclosing these three equilibria, as shown in Figures 8(a) and (b). For γ03<γ<γ01, the positive equilibria E1* and E2* are still unstable, but the positive equilibrium E3* becomes stable, as shown in Figures 8(c) and 8(d) with γ=0.3.The trajectory graphs and phase portrait of system (7) without delay for (α,β)∈D3 with α=0.01 and β=0.24. (a)-(b) γ=0.1<min{γ01,γ03}. The initial value is (0.626,0.15). All of the three equilibria are unstable and there exists a stable periodic orbit. (c)-(d) γ03<γ=0.3<γ01. The initial value is (0.049,0.0116). The positive equilibria E1* and E2* are unstable and the positive equilibrium E3* is stable.
(a)
(b)
(c)
(d)Whenγ>max{γ01,γ03}, the positive equilibria E1* and E3* both become asymptotically stable. Taking γ=0.5, Figures 9(a) and 9(b) illustrate this theoretical result. The critical values of the delay for these two positive equilibria E1* and E3* are τ01≐0.3782 and τ03≐2.3129, respectively. Thus the positive equilibrium E1* is asymptotically stable for τ<τ01 and unstable for τ>τ01, and the positive equilibrium E3* is asymptotically stable for τ<τ03 and unstable for τ>τ03. According to the procedure of Section 3.3, we obtain that K1≐0.0851>0, and K2≐49.9838>0 for E1* and τ01, and K1≐0.1197>0, and K2≐13.9616>0 for E3* and τ03. So, the delay-induced Hopf bifurcations for these two positive equilibria are both subcritical. Taking τ=0.37<min{τ01,τ03}, the positive equilibrium E1* and E3* are both asymptotically stable. Near the positive equilibria E1*, the subcritical Hopf bifurcation occurs. Figures 9(c) and 9(d) illustrate the existence of this subcritical Hopf bifurcation and show that there is an orbit connecting the unstable periodic orbit to the stable equilibrium E3*. When τ>max{τ01,τ03}, the positive equilibria E1* and E3* both becomes unstable. There is a large-amplitude stable periodic orbit enclosing the three equilibria, as shown in Figures 9(e) and 9(f) for τ=2.5.The trajectory graphs and phase portrait of system (6) for (α,β)∈D3 with α=0.01,  β=0.24, and γ=0.5. (a)-(b) When τ=0, the positive equilibria E1* and E3* are both asymptotically stable. The initial values are (0.2,0.05) and (0.6,0.05). (c)-(d) When τ=0.37<min{τ01,τ03}, the positive equilibria E1* and E3* are both asymptotically stable and there exists a subcritical Hopf bifurcation near E1*. Here, the initial values is (0.042,0.013). (e)-(f) When τ=2.5>max{τ01,τ03}. All of the three equilibria become unstable and there exists a large-amplitude stable periodic orbit enclosing the three equilibria. Here, the initial value is (0.6,0.15).
(a)
(b)
(c)
(d)
(e)
(f)


## 5. Conclusions

In this paper, we have considered a predator-prey system with the hunting delay. By analysing the concerned characteristic equations, local asymptotic stability of the various positive equilibria of the system with or without delay is discussed. The dynamics of the system can be described in the plane of the parametersα and β. For the system without delay, the bistability (existence of two stable equilibria) occurs in some parameter regions. We have shown that the ratio (γ=s/r) of the intrinsic growth rates of the predator and prey is the key parameter which has an important influence on the dynamics of the system. In the region where the system has only one positive equilibrium and the ratio γ can affect the stability of this positive, the ratio γ can induce large-amplitude periodic oscillations of both species and there are two limit cycles (one is stable and the other is unstable) near the critical value of the ratio via subcritical Hopf bifurcation. In the region where the system has three positive equilibria, the dynamics of the system can evolve from a heteroclinic orbit connecting one positive equilibrium to a large-amplitude limit cycle enclosing three equilibria (see Figures 8(a) and 8(b)), a heteroclinic orbit connecting two equilibria (see Figures 8(c) and 8(d)), two stable equilibria with the increase of the ratio γ (see Figures 9(a) and 9(b)).When the hunting delay is introduced by the system, its effect on the dynamics of all positive equilibria is investigated. We first show that whenα<β (or equivalent to kmn<ar), the unique positive equilibrium is global stability for any nonnegative hunting delay. Then the absolute investigate the influence of the delay on the local stability of the positive equilibria. For the case when the system has the unique positive equilibrium, depending on the choice of the parameters α and β, the delay can induce the emergence of two limit cycles with one stable subcritical Hopf bifraction and another unstable subcritical Hopf bifurcation (see Figures 4(c) and 4(d)) and induce the emerge of semistable limit cycle which is inward unstable and outward stable (see Figures 7(c) and 7(d)). For the case when the system has the three positive equilibria, delay can induce a heteroclinic orbit connecting a limit cycle to a stable equilibrium (see Figures 9(c) and 9(d)). The analytical expressions that determine the properties of bifurcating periodic solutions are also obtained by using the normal form theory and the centre manifold theorem. It is also worthy to mention that for the case without delay, the size of limit cycle increases with the decreasing magnitude of the ratio of the intrinsic growth rates of the predator and prey, but for the case with delay, the size of limit cycle increases with the increasing magnitude of the hunting delay.



---

*Source: 152459-2013-05-08.xml*