Bogdanov-Takens and Triple Zero Bifurcations of a Delayed Modified Leslie-Gower Predator Prey System

Abstract

A delayed modified Leslie-Gower predator prey system with nonlinear harvesting is considered. The existence conditions that an equilibrium is Bogdanov-Takens (BT) or triple zero singularity of the system are given. By using the center manifold reduction, the normal form theory, and the formulae developed by Xu and Huang, 2008 and Qiao et al., 2010, the normal forms and the versal unfoldings for this singularity are presented. The Hopf bifurcation of the system at another interior equilibrium is analyzed by taking delay (small or large) as bifurcation parameter.

1. Introduction

For a more detailed study on the properties of the predator prey systems, the multiple bifurcations for some systems (ODE) with more interior equilibria are investigated by many authors, see [1–4] for example. To deal with this type of systems, the more difficult problem is how to obtain the normal form of the system at its degenerate equilibrium, that is, BT bifurcation.

When introducing the time delay into this type of systems, using the methods developed by [5], the authors in [6, 7] have researched the BT bifurcation of some predator prey systems (DDE). Their results show that time delay may have an effect or not on the BT bifurcation.

Recently, papers [8–10] have considered the triple zero bifurcation of some delay differential equations, depending on the parameters in the original system; some interesting bifurcation results are obtained. But we find that there are few results about the triple zero bifurcation for predator prey systems.

Summarizing the above references, we will consider the following predator prey system with Michaelis-Menten type (nonlinear) prey harvesting: ̇𝑥𝑥=𝑟𝑥1−𝐾−𝑎₁𝑥𝑦−𝑛+𝑥𝑞𝐸𝑥𝑘₁𝐸+𝑘₂𝑥,𝑦̇=𝑠₁𝑦𝑎1−²𝑦(𝑡−𝜏)(),𝑛+𝑥𝑡−𝜏 (1) where 𝑥 and 𝑦 denote the prey and predator populations, respectively. 𝜏 represents the negative feedback of the predator’s density. For complete reason, we give the biological meaning of the parameters, one can be seeing them in [11]. 𝑟 and 𝐾 are intrinsic growth rate and environmental carrying capacity for the prey, respectively. 𝑎₁ is the maximum value of the per capita reduction rate of prey, 𝑛 measures the extent to which the environment provides protection to prey and predator respectively, 𝑠₁ measures the growth rate of the predator species, and 𝑠₁𝑎₂ is the maximum value of the per capita reduction rate of predator. 𝑝(𝑥)=𝑞𝐸𝑥/(𝑘₁𝐸+𝑘₂𝑥) represents Michaelis-Menten type harvesting, 𝑞 is the catchability coefficient, 𝐸 is the effort applied to harvest the prey species, and 𝑘₁ and 𝑘₂ are suitable constants.

The authors in [12, 13] have studied system (1) without prey harvesting 𝑝(𝑥), respectively, in [12], the global stability and persistence of the system are investigated. In [13], by using the Hopf bifurcation theorem and taking the delay as a parameter of bifurcation for small and large cases, the existence of the bifurcated limit cycle around a boundary equilibrium or an interior equilibrium is mainly considered.

For system (1) with 𝜏=0, the authors in [11] have given detailed analysis about the existence of the multiple bifurcations (including BT bifurcation) depending on the parameters of the system.

For computation simplicity, we first rescale system (1).

Let 𝑡=𝑟𝑡, 𝑥(𝑡)=𝑥(𝑡)/𝐾, 𝑦(𝑡)=𝑎₁𝑦(𝑡)/𝐾, and  𝜏=𝑟𝜏; then dropping the bars we obtain 𝑥̇=𝑥1−𝑥−𝛼𝑦−ℎ𝑚+𝑥,̇()𝑐+𝑥𝑦=𝜌𝑦1−𝛽𝑦𝑡−𝜏(),𝑚+𝑥𝑡−𝜏 (2) where 𝛼=1/𝑟, 𝛽=𝑎₂/𝑎₁, 𝑚=𝑛/𝐾, 𝜌=𝑠₁/𝑟, ℎ=𝑞𝐸/𝑟𝑘₂𝐾, and 𝑐=𝑘₁𝐸/𝑘₂𝐾.

System (2) with initial conditions is 𝜙(𝜗)𝜓(𝜗)[),𝜙(0)𝜓(0)⩾0,⩾0,𝜗∈−𝜏,0>0,>0, (3) where 𝜙(𝜗) and 𝜓(𝜗) are all continuous bounded functions in the interval [−𝜏,0).

From [11], we know that (2) has interior equilibrium 𝐸=(𝑥,𝑦) if ℎ=ℎ^∗=14𝛼1−𝛽+𝑐²=(𝑐+𝑥)²,𝛼𝛽+𝑐<1, (4) where 𝑥=(1/2)(1−𝑐−𝛼/𝛽), 𝑦=(𝑚+𝑥)/𝛽 and interior equilibrium 𝐸^∗=(𝑥^∗,𝑦^∗) if ℎ<ℎ^∗, where 𝑥^∗=(1/2)(1−𝑐−𝛼/𝛽+√Δ), 𝑦^∗=(𝑚+𝑥^∗)/𝛽, and Δ=(𝛼/𝛽+𝑐−1)²−4𝑐(𝛼/𝛽+ℎ/𝑐−1).

In this paper, for system (2), we will mainly consider the BT and triple zero bifurcations at 𝐸 and the Hopf bifurcation at 𝐸^∗. It is easy to see that this system is with six parameters which will let our work become more challenging. When dealing with the BT and triple bifurcations of the delay systems, the core problem is to change delay systems as ordinary differential systems (ODEs).

The concrete organization of the paper is as follows: in Section 2, we will give the conditions under which the equilibrium 𝐸 is a BT singularity, and a universal unfolding will be exhibited; in Section 3, when 𝐸 is a triple zero, the universal unfolding will be presented, and in Section 4, some Hopf bifurcation results at 𝐸^∗ will be obtained.

2. Bogdanov-Takens Bifurcation

System (2) also can be written as 𝑥̇=𝜏𝑥1−𝑥−𝛼𝑦−ℎ𝑚+𝑥,̇()𝑐+𝑥𝑦=𝜏𝜌𝑦1−𝛽𝑦𝑡−1().𝑚+𝑥𝑡−1 (5) Linearizing system (5) at 𝐸 yields the following linear system 𝑥̇₁=𝜏𝛼𝑥𝛽(𝑚+𝑥)𝑥₁(𝑡)−𝜏𝛼𝑥𝑚+𝑥𝑥₂(𝑡),𝑥̇₂=𝜏𝜌𝛽𝑥₁(𝑡−1)−𝜏𝜌𝑥₂(𝑡−1), (6) the corresponding characteristic equation is 𝐹(𝜆)=𝜆²+𝜏𝜌𝑒^−𝜆−𝜏𝛼𝑥𝛽(𝑚+𝑥)𝜆. (7) Evidently, 𝜆=0 is a double zero eigenvalue if 1𝜏≠𝜌𝛼,𝜌=𝑥𝛽(𝑚+𝑥); (8)𝜆=0 is a triple zero eigenvalue if 1𝜏=𝜌𝛼,𝜌=𝑥𝛽(𝑚+𝑥). (9) It is easy to prove that in the above two cases the rest eigenvalues all have negative real parts.

Under the conditions (4) and (8), let 𝑥₁=𝑥−𝑥, and let 𝑥₂=𝑦−𝑦; then the Taylor expansion system (5) at 𝐸 is 𝑥̇₁=𝜏𝜌𝑥₁(𝑡)−𝜏𝛽𝜌𝑥₂(𝑡)+_{𝑖+𝑗≥2}1𝑓𝑖!𝑗!^(1)𝑖𝑗𝑥^𝑖₁(𝑡)𝑥^𝑗₂(𝑡),𝑥̇₂=𝜏𝜌𝛽𝑥₁(𝑡−1)−𝜏𝜌𝑥₂(𝑡−1)+_{𝑖+𝑗+𝑘≥2}1𝑓𝑖!𝑗!𝑘!^{(2)𝑖𝑗𝑘}𝑥^𝑖₁(𝑡−1)𝑥^𝑗₂(𝑡−1)𝑥^𝑘₂(𝑡), (10) where 𝑖,𝑗,𝑘≥0,  𝑓⁽¹⁾_𝑖𝑗=(𝜕^𝑖+𝑗𝑓⁽¹⁾/𝜕^𝑖𝑥𝜕^𝑗𝑦)|_(𝑥,𝑦), 𝑓⁽²⁾_𝑖𝑗𝑘=(𝜕^{𝑖+𝑗+𝑘}𝑓⁽²⁾/𝜕^𝑖𝑥𝜕^𝑗𝑦𝜕^𝑘𝑦₀)|_(𝑥,𝑦), 𝑓⁽¹⁾=𝜏𝑥(1−𝑥−𝛼𝑦/(𝑚+𝑥)−ℎ/(𝑐+𝑥)), 𝑓⁽²⁾=𝜏𝑦₀(1−𝛽𝑦/(𝑚+𝑥)).

In the following, we first give the normal form of the system (10) at the singularity (0,0). Reference [6], we first rewrite system (10) as 𝑋̇(𝑡)=𝐿(𝑋_𝑡), here 𝑋(𝑡)=(𝑥₁(𝑡),𝑥₂(𝑡)),  𝐿(𝜙)=𝐿^𝜙1_𝜙⁽⁻¹⁾₂₍₀₎, and 𝜙=(𝜙₁,𝜙₂). By the normal form theory developed by Faria and Magalhaes [5], one can obtain the center manifold of this system at the origin which is two-dimensional and system can be reduced to an ODE in the plane.

Define 𝐴₀ to be the infinitesimal generator of system. Consider 𝜆={0} and let 𝑃 denote the invariant space of 𝐴₀ associated with the eigenvalue 𝜆=0, using the formal adjoint theory in [5], the phase space 𝐶₁ can be decomposed by 𝜆 as 𝐶₁=𝑃⊕𝑄. Let Φ and Ψ be the bases for 𝑃 and 𝑃^∗, respectively, and be let them normalized such that ⟨Ψ,Φ⟩=𝐼, Φ̇=Φ𝐽,and Ψ̇=−𝐽Ψ, where Φ and Ψ are 2×2 matrices, where 𝐽=⁰¹₀₀.

Next, we will find the Φ(𝜃) and Ψ(𝑠) based on the techniques developed by [14].

Lemma 1 (see Xu and Huang [14]). The bases of 𝑃 and their dual space 𝑃^∗ have the following representations: 𝑃=spanΦ,Φ(𝜃)=𝜑₁(𝜃),𝜑₂(𝜃)𝑃,−1≤𝜃≤0,_∗=spanΨ,Ψ(𝑠)=col𝜓₁(𝑠),𝜓₂(𝑠),0≤𝑠≤1, (11) where 𝜑₁(𝜃)=𝜑⁰₁∈𝑅^𝑛⧵{0}, 𝜑₂(𝜃)=𝜑⁰₂+𝜑⁰₁𝜃,𝜑⁰₂∈𝑅^𝑛, and 𝜓₂(𝑠)=𝜓⁰₂∈𝑅^𝑛∗⧵{0}, 𝜓₁(𝑠)=𝜓⁰₁−𝑠𝜓⁰₂, 𝜓⁰₁∈𝑅^𝑛∗, which satisfy(1)(𝐴+𝐵)𝜑⁰₁=0, (2)(𝐴+𝐵)𝜑⁰₂=(𝐵+𝐼)𝜑⁰₁, (3)𝜓⁰₂(𝐴+𝐵)=0, (4)𝜓⁰₁(𝐴+𝐵)=𝜓⁰₂(𝐵+𝐼), (5)𝜓⁰₂𝜑⁰₂−(1/2)𝜓⁰₂𝐵𝜑⁰₁+𝜓⁰₂𝐵𝜑⁰₂=1, (6)𝜓⁰₁𝜑⁰₂−(1/2)𝜓⁰₁𝐵𝜑⁰₁+𝜓⁰₁𝐵𝜑⁰₂+(1/6)𝜓⁰₂𝐵𝜑⁰₁−(1/2)𝜓⁰₂𝐵𝜑⁰₂=0.

By system (6), we know that 𝐴=𝜏𝛼𝑥𝛽(𝑚+𝑥)−𝜏𝛼𝑥𝑚+𝑥=,00𝜏𝜌−𝜏𝛽𝜌00𝐵=00𝜏𝜌𝛽.−𝜏𝜌 (12)

By Lemma 1, we have Φ(𝜃)=𝛽𝛽2−𝜏𝜌−2²12𝜏𝜌𝜏𝜌−1+𝛽𝜃2−𝜏𝜌2≜𝜏𝜌−1+𝜃𝛽𝑚₁+𝛽𝜃1𝑚₂Ψ(𝑠)=−+𝜃,−1≤𝜃≤0,𝜏𝜌𝛽1−𝜏𝜌𝑠1−𝜏𝜌𝑠𝜏𝜌−1𝜏𝜌𝛽1−𝜏𝜌𝜏𝜌𝜏𝜌−1,0≤𝑠≤1, (13) where 𝑚₁=𝛽[2−(𝜏𝜌−2)²]/2𝜏𝜌(𝜏𝜌−1), 𝑚₂=(2−𝜏𝜌)/2(𝜏𝜌−1), Φ(𝜃)𝑧=𝛽𝑚₁+𝛽𝜃1𝑚₂𝑧+𝜃𝑧¹₂=𝛽𝑧₁+𝑚₁+𝛽𝜃𝑧₂𝑧₁+𝑚₂+𝜃𝑧₂, (14) therefore, 𝜙₁(0)=𝛽𝑧₁+𝑚₁𝑧₂,𝜙₂(0)=𝑧₁+𝑚₂𝑧₂,𝜙₁(−1)=𝛽𝑧₁+𝑚₁−𝛽𝑧₂,𝜙₂(−1)=𝑧₁+𝑚₂−1𝑧₂. (15)

Taking a similar theory in [6], system (6) in the center manifold is 𝑧̇=𝐽𝑧+Ψ(0)𝐹(Φ𝑧), (16) where 𝐹𝜙=𝜏_{𝑖+𝑗≥2}1𝑓𝑖!𝑗!^(1)𝑖𝑗𝜙₁(0)^𝑖𝜙₂(0)^𝑗𝜏𝜌_{𝑖+𝑗+𝑘≥2}1𝑓𝑖!𝑗!𝑘!^{(2)𝑖𝑗𝑘}𝜙^𝑖₁(−1)𝜙^𝑗₂(−1)𝜙^𝑘₂(0). (17)

Using (10) and (15), system (16) also can be written as 𝑧̇₁=𝑧₂+𝛽𝜏𝜌−1𝜏𝜌(𝑚+𝑥)𝑧²₂𝑧̇+h.o.t,₂=𝑔₁𝑧²₁+𝑔₂𝑧₁𝑧₂+𝑔₃𝑧²₂+h.o.t, (18) where 𝑔₁=−𝜏²𝛽𝜌𝑥/(1−𝜏𝜌)(𝑐+𝑥), 𝑔₂=𝜏𝛼𝑚/(1−𝜏𝜌)(𝑚+𝑥)²−2𝑚₁𝜌𝜏²𝑥/(1−𝜏𝜌)(𝑐+𝑥), 𝑔₃=(𝜏𝜌/𝛽(1−𝜏𝜌))[𝛼𝑚₁𝑚/𝜌(𝑚+𝑥)²−𝜏𝑥𝑚²₁/(𝑐+𝑥)]+𝛽/(𝑚+𝑥). Through a series of transformation, system (18) becomes 𝑧̇₁=𝑧₂𝑧̇+h.o.t,₂=𝑑₁𝑧²₁+𝑑₂𝑧²₂+h.o.t, (19) where 𝑑₁=𝑔₁ and 𝑑₂=𝑔₂. If 𝑑₁𝑑₂≠0, we have the following theorem.

Theorem 2. Let (4), (8), and 𝜏≠𝛽𝑚(𝑐+𝑥)/2𝑚₁𝑥²(𝑚+𝑥) hold. Then, the equilibrium 𝐸 of system (5) is a BT singularity.

Next, we are interested in giving a versal unfolding for system (5) at BT singularity. Choosing ℎ and 𝜌 as bifurcation parameters and incorporating ℎ+𝜆₁ and 𝜌+𝜆₂ to system (5), where 𝜆₁ and 𝜆₂ vary in a small neighborhood of (0,0), we obtain 𝑥̇=𝜏𝑥1−𝑥−𝛼𝑦−𝑚+𝑥ℎ+𝜆₁,̇𝑐+𝑥𝑦=𝜏𝜌+𝜆₂𝑦()1−𝛽𝑦𝑡−1().𝑚+𝑥𝑡−1 (20)

Let 𝑦₁=𝑥−𝑥 and 𝑦₂=𝑦−𝑦. Then, system (20) becomes 𝑦̇₁𝜏=−𝑥𝜆₁𝑐+𝑥𝛼+𝜏𝑥𝛽(𝑚+𝑥)+𝑥𝜆₁ℎ−𝜆₁𝑐+𝑥𝑦₁(𝑡)−𝜏𝛼𝑥𝑚+𝑥𝑦₂(𝑡)+𝜏𝛼𝑚𝛽(𝑚+𝑥)²+𝜏𝑐𝜆₁−ℎ𝑥ℎ(𝑐+𝑥)𝑦²₁(𝑡)−𝜏𝛼𝑚(𝑚+𝑥)²𝑦₁(𝑡)𝑦₂(𝑡)𝑦̇+h.o.t,₂=𝜏𝜌^∗+𝜆₂×1𝛽𝑦₁(𝑡−1)−𝑦₂()+1𝑡−1𝑚+𝑥𝑦₁(𝑡−1)𝑦₂(𝑡)+1𝑚+𝑥𝑦₁(𝑡−1)𝑦₂()−1𝑡−1𝛽(𝑚+𝑥)𝑦²₁()−𝛽𝑡−1𝑚+𝑥𝑦₂(𝑡)𝑦₂(𝑡−1)+h.o.t. (21) Then, system (21) can be decomposed as 𝑧̇=𝐵₁+𝐵₂𝑧+Ψ₀𝐺(),̇Φ𝑧+𝑢𝑢=𝐴_𝑄′𝑢+(𝐼−𝜋)𝑋₀×𝐵₀+𝐵₂Φ₀𝑧+𝑢(0)+𝐺(Φ𝑧+𝑢), (22) where 𝑋₀(𝜃)=𝐵𝐼,𝜃=00,−1≤𝜃<0,₀=−𝜏𝑥𝜆₁𝑐+𝑥0,𝐵₁=Ψ(0)𝐵₀=01𝜏𝜌𝛽1−𝜏𝜌𝜏𝜌−𝜏𝜏𝜌−1𝑥𝜆₁𝑐+𝑥0=0−𝜏²𝜌𝑥𝜆₁𝛽(𝑐+𝑥),𝐺𝜙=1−𝜏𝜌𝜏_{𝑖+𝑗≥2}1𝑔𝑖!𝑗!^(1)𝑖𝑗𝜙₁(0)^𝑖𝜙₂(0)^𝑗𝜏𝜌+𝜆_{2𝑖+𝑗+𝑘≥2}1𝑔𝑖!𝑗!𝑘!^{(2)𝑖𝑗𝑘}𝜙^𝑖₁(−1)𝜙^𝑗₂(−1)𝜙^𝑘₂(0),𝑔⁽¹⁾_𝑖𝑗=𝜕^𝑖+𝑗𝑔⁽¹⁾𝜕^𝑖𝑥𝜕^𝑗𝑦_{(𝑥,𝑦,𝜆1)},𝑔⁽²⁾_𝑖𝑗𝑘=𝜕^{𝑖+𝑗+𝑘}𝑓⁽²⁾𝜕^𝑖𝑥𝜕^𝑗𝑦𝜕^𝑘𝑦_0(𝑥,𝑦),𝑔⁽¹⁾=𝑥1−𝑥−𝛼𝑦−ℎ𝑚+𝑥^∗+𝜆₁,𝑔𝑐+𝑥₍₂₎=𝑦₀1−𝛽𝑦,𝐵𝑚+𝑥₂=𝜆01+²𝜌𝜏²𝜌𝑐𝜆₁ℎ𝜏𝜏𝜌−1²𝜌𝜏𝜌−1𝑐𝑚₁𝜆₁+𝜆ℎ𝛽².𝜏𝜌 (23)

Then, the normal form of system (21) at 𝐸 is 𝑧̇=𝐵₁+𝐵₂𝑧+Ψ(0)𝐺(Φ𝑧); that is, 𝑧̇₁=𝑧₂+𝜆₂𝜌𝑧₂+𝛽𝜌+𝜆₂𝜏𝜌−1𝜏𝜌²(𝑚+𝑥)𝑧²₂𝑧̇+h.o.t,₂𝜏=−²𝑥𝜌𝜆₁𝛽(𝑐+𝑥)+𝜏1−𝜏𝜌²𝑐𝜌𝜆₁ℎ𝑧𝜏𝜌−1¹+𝜏²𝜌𝜏𝜌−1𝑐𝑚₁𝜆₁+𝜆ℎ𝛽²𝑧𝜏𝜌²+𝑃₁𝑧²₁+𝑃₂𝑧₁𝑧₂+𝑃₃𝑧²₂+h.o.t, (24) where 𝑃₁=𝜏²𝛽𝜌𝑐𝜆₁−ℎ𝑥ℎ1−𝜏𝜌(𝑐+𝑥),𝑃₂=𝜏𝛼𝑚1−𝜏𝜌(𝑚+𝑥)²+2𝑚₁𝜌𝜏²𝑐𝜆₁−ℎ𝑥ℎ1−𝜏𝜌(𝑐+𝑥),𝑃₃=𝜏𝜌𝛽1−𝜏𝜌𝛼𝑚₁𝑚𝜌(𝑚+𝑥)²+𝜏𝑐𝜆₁−ℎ𝑥𝑚²₁ℎ(𝑐+𝑥)+𝛽𝜌+𝜆₂𝜌(𝑚+𝑥). (25)

Following the normal form formula in Kuznetsov [15], system (24) can be reduced to 𝑧̇₁=𝑧₂𝑧̇+h.o.t,₂=𝛾₁+𝛾₂𝑧₁+𝑑₁𝑧²₁+𝑑₂𝑧₁𝑧₂+h.o.t, (26) where 𝛾₁𝜏=−²𝑥𝜌𝜆₁𝛽(𝑐+𝑥,𝛾)1−𝜏𝜌₂=2𝜏²𝑥𝑐+𝑥2𝛽²𝑥²−𝜏𝑐𝑚₁𝛼𝑥²(𝑐+𝑥)𝜌𝑚𝛽²(𝑐+𝑥)−2𝑚₁𝜏𝛼𝑥³+𝜌³𝜏³𝑚₁𝑚𝛽²−𝜏𝑚²₁𝛼𝑥³(𝑐+𝑥)𝛼𝑥𝛽²1−𝜏𝜌²(𝑐+𝑥)𝜆₁−2𝛼𝜏²𝑥³𝜌1−𝜏𝜌𝑚𝛽(𝑐+𝑥)−2𝑚₁𝜏𝑥²(𝑚+𝑥𝜆)². (27)

Then, system (5) exists in the following bifurcation curves in a small neighborhood of the origin in the (𝜆₁,𝜆₂) plane.

Theorem 3. Let (4), (8), and 𝜏≠𝛽𝑚(𝑐+𝑥)/2𝑚₁𝑥²(𝑚+𝑥) hold. System (5) admits the following bifurcations:(i)a saddle-node bifurcation curve SN={(𝜆₁,𝜆₂);𝛾₁=(1/4𝑑₁)𝛾²₂};(ii)a Hopf bifurcation curve 𝐻={(𝜆₁,𝜆₂);𝜆₂=0,𝛾₂<0};(iii)a homoclinic bifurcation curve 𝐻𝐿={(𝜆₁,𝜆₂);𝛾₁=−(6/25𝑑₁)𝛾²₂,𝛾₂<0}.

3. Triple-Zero Bifurcation

From Section 2, we know that under the conditions (4) and (9) the equilibrium 𝐸 of system (5) is a triple zero singularity. In the following reference, from the work of [8, 9] we will give the triple zero bifurcation at 𝐸.

Let 𝑥₁=𝑥−𝑥 and 𝑥₂=𝑦−𝑦; then, system (5) becomes 𝑥̇₁=𝜏1−2𝑥−𝛼𝑚𝛽(𝑚+𝑥)−ℎ𝑐(𝑐+𝑥)²𝑥₁(𝑡)−𝜌𝛽𝑥₂(𝑡)+_{𝑖+𝑗≥2}1𝑓𝑖!𝑗!^(1)𝑖𝑗𝑥^𝑖₁(𝑡)𝑥^𝑗₂(𝑡),𝑥̇₂1=𝜏𝜌𝛽𝑥₁(𝑡−1)−𝑥₂(𝑡−1)+_{𝑖+𝑗+𝑘≥2}1𝑓𝑖!𝑗!𝑘!^{(2)𝑖𝑗𝑘}𝑥^𝑖₁(𝑡−1)𝑥^𝑗₂(𝑡−1)𝑥^𝑘₂(𝑡). (28)

To determine a versal unfolding for the original system (28) at 𝐸, we choose 𝜏, 𝜌, and ℎ as bifurcation parameters, and let them become 𝜏+𝜇₁,𝜌+𝜇₂, and ℎ+𝜇₃, respectively, where 𝜇₁,𝜇₂, and 𝜇₃ vary in a small neighborhood of (0,0,0); then (28) can be written as 𝑥̇₁=𝜏+𝜇₁×𝜌−𝑐𝜇₃ℎ𝑥₁(𝑡)−𝜌𝛽𝑥₂(𝑡)+𝛼𝑚𝛽(𝑚+𝑥)²+𝑐𝜇₃−ℎ𝑥ℎ(𝑐+𝑥)𝑥²₁(𝑡)−𝛼𝑚(𝑚+𝑥)²𝑥₁(𝑡)𝑥₂(𝑡𝑥̇)+h.o.t,₂=𝜏+𝜇₁𝜌+𝜇₂×1𝛽𝑥₁(𝑡−1)−𝑥₂()+1𝑡−1𝑚+𝑥𝑥₁(𝑡−1)𝑥₂(𝑡)+1𝑚+𝑥𝑥₁(𝑡−1)𝑥₂()−1𝑡−1𝛽(𝑚+𝑥)𝑥²₁()−𝛽𝑡−1𝑚+𝑥𝑥₂(𝑡)𝑥₂(𝑡−1)+h.o.t. (29)

Next, we need to find the expressions of Φ(𝜃) and Ψ(𝑠) based on the techniques developed by [9].

Lemma 4 (see Qiao et al. [9]). The bases of 𝑃 and their dual space 𝑃^∗ have the following representations: 𝑃=spanΦ,Φ(𝜃)=𝜑₁(𝜃),𝜑₂(𝜃),𝜑₃(𝜃),𝑃−1≤𝜃≤0,_∗=spanΨ,Ψ(𝑠)=col𝜓₁(𝑠),𝜓₂(𝑠),𝜓₃(𝑠),0≤𝑠≤1, (30) where 𝜑₁(𝜃)=𝜑⁰₁∈𝑅^𝑛⧵{0}, 𝜑₂(𝜃)=𝜑⁰₂+𝜑⁰₁𝜃, 𝜑₃(𝜃)=𝜑⁰₃𝜑⁰₂𝜃+(1/2)𝜑⁰₁𝜃², 𝜑⁰₂,𝜑⁰₃∈𝑅^𝑛, and 𝜓₃(𝑠)=𝜓⁰₃∈𝑅^𝑛∗⧵{0}, 𝜓₂(𝑠)=𝜓⁰₂−𝑠𝜓⁰₃, 𝜓₁(𝑠)=𝜓⁰₁−𝑠𝜓⁰₂+(1/2)𝑠²𝜓⁰₃, 𝜓⁰₁,𝜓⁰₂∈𝑅^𝑛∗, which satisfy (1)(𝐴+𝐵)𝜑⁰₁=0, (2)(𝐴+𝐵)𝜑⁰₂=(𝐵+𝐼)𝜑⁰₁, (3)(𝐴+𝐵)𝜑⁰₃=(𝐵+𝐼)𝜑⁰₂−(1/2)𝐵𝜑⁰₁, (4)𝜓⁰₃(𝐴+𝐵)=0, (5)𝜓⁰₂(𝐴+𝐵)=𝜓⁰₃(𝐵+𝐼), (6)𝜓⁰₁(𝐴+𝐵)=𝜓⁰₂(𝐵+𝐼)−(1/2)𝜓⁰₃𝐵, (7)𝜓⁰₃(𝐵+𝐼)𝜑⁰₃−(1/2)𝜓⁰₃𝐵𝜑⁰₂+(1/6)𝜓⁰₃𝐵𝜑⁰₁=1, (8)𝜓⁰₂(𝐵+𝐼)𝜑⁰₃−(1/2)𝜓⁰₂𝐵𝜑⁰₂+(1/6)𝜓⁰₂𝐵𝜑⁰₁−(1/2)𝜓⁰₃𝐵𝜑⁰₃+(1/6)𝜓⁰₃𝐵𝜑⁰₂−(1/24)𝜓⁰₃𝐵𝜑⁰₁=0, (9)  𝜓⁰₁(𝐵+𝐼)𝜑⁰₃−(1/2)𝜓⁰₁𝐵𝜑⁰₂+(1/6)𝜓⁰₁𝐵𝜑⁰₁−(1/2)𝜓⁰₂𝐵𝜑⁰₃+(1/6)𝜓⁰₂𝐵𝜑⁰₂−(1/24)𝜓⁰₂𝐵𝜑⁰₁ + (1/6)𝜓⁰₃𝐵𝜑⁰₃−(1/24)𝜓⁰₃𝐵𝜑⁰₂ + (1/120)𝜓⁰₃𝐵𝜑⁰₁=0.

For system (29), one can see that 𝐴+𝐵=𝜏𝜌11−𝛽𝛽=1−11−𝛽𝛽−1. (31) By Lemma 4, we can obtain Φ(𝜃)=1𝛽𝛽𝜃𝛽1+2𝜃²11−1+𝜃1−𝜃+2𝜃²Ψ(𝑠)=−,−1≤𝜃≤0,17−18𝛽2𝑠+𝑠3𝛽𝛽²35+182𝑠3−𝑠²2−3𝛽2𝑠𝛽−232+2𝑠𝛽,−20≤𝑠≤1, (32) such that ⟨Ψ,Φ⟩=𝐼, Φ̇=Φ𝐽, and Ψ̇=−𝐽Ψ, where 𝐽 is given by 𝐽=⁰¹⁰⁰⁰¹₀₀₀.

Using the similar methods as used in [9], system (29) can be rewritten as 𝑥̇(𝑡)=𝐴𝜇𝑥(𝑡)+𝐵𝜇𝑥(𝑡−1)+𝐹𝑥(𝑡),𝑥(𝑡−1),𝜇, (33) where 𝜇=(𝜇₁,𝜇₂,𝜇₃)∈𝑅³ is a parameter vector, 𝑥∈𝑅². Let 𝐹𝑥(𝑡)()=(𝐴(𝜆)(0𝑥(𝑡)+𝐵𝜇(0)𝑥()𝑥(𝑡)(),,𝑥𝑡−1,𝜇−𝐴))−𝐵𝑡−1+𝐹,𝑥𝑡−1,𝜇 (34) then, we can expand 𝐹(𝑥(𝑡), 𝑥(𝑡−1), 𝜇)=𝐴₁𝜇₁𝑥(𝑡) + 𝐴₂𝜇₂𝑥(𝑡) + 𝐴₃𝜇₃𝑥(𝑡) + 𝐵₁𝜇₁𝑥(𝑡−1) + 𝐵₂𝜇₂𝑥(𝑡−1) + 𝐵₃𝜇₃𝑥(𝑡−1) + 𝐸₁𝑥₁(𝑡)𝑥(𝑡−1) + 𝐸₂𝑥₂(𝑡)𝑥(𝑡−1) + 𝐹₁𝑥₁(𝑡)𝑥(𝑡) + 𝐹₂𝑥₂(𝑡)𝑥(𝑡) + 𝐺₁𝑥₁(𝑡−1)𝑥(𝑡−1) + 𝐺₂𝑥₂(𝑡−1)𝑥(𝑡−1) + h.o.t, where 𝐴₁=𝜌−𝛽𝜌00,𝐴₂=0000,𝐴₃=−𝑐𝜏ℎ0,𝐵00₁=𝜌00𝛽−𝜌,𝐵₂=𝜏00𝛽−𝜏,𝐵₃=,𝐸0000₁=0000,𝐸₂=100𝑚+𝑥−𝛽𝑚+𝑥,𝐹₁=𝜏𝛼𝑚𝛽(𝑚+𝑥)²−𝜏𝑥𝑐+𝑥−𝜏𝛼𝑚2(𝑚+𝑥)²,𝐹00₂=−𝜏𝛼𝑚2(𝑚+𝑥)²0,𝐺00₁=−100𝛽(𝑚+𝑥)12(𝑚+𝑥),𝐺₂=1002(𝑚+𝑥)0,𝜑⁰₁=𝛽1,𝜑⁰₂=0−1,𝜑⁰₃=𝛽1,𝜓⁰₁=−1718𝛽35,𝜓18₀₂=2−23𝛽3,𝜓⁰₃=2𝛽−2. (35)

Following the formula of Theorem 3.1 in [9], the normal form with versal unfolding of system (29) on the center manifold takes the following form: 𝑧̇₁=𝑧₂,𝑧̇₂=𝑧₃,𝑧̇₃=𝑙₁𝑧₁+𝑙₂𝑧₂+𝑙₃𝑧₃+𝜂₁𝑧²₁+𝜂₂𝑧²₂+𝜂₃𝑧₁𝑧₂+𝜂₄𝑧₁𝑧₃+h.o.t, (36) where 𝑙₁=−2𝑐𝜏ℎ𝜇₃,𝑙₂=−4𝜏𝜇₂−2𝑐𝜏𝜇3ℎ³,𝑙₃=2𝜌𝜇₁+4𝜏3𝜇₂−19𝑐𝜏𝜇18ℎ³,𝜂₁=−2𝜏𝑥𝛽𝑐+𝑥,𝜂₂=17𝜏𝑥𝛽9(𝑐+𝑥)+2𝜏𝛼𝑚3(𝑚+𝑥)²,𝜂₃=−4𝜏𝑥𝛽3(𝑐+𝑥)+2𝜏𝛼𝑚(𝑚+𝑥)²,𝜂₄=2𝜏𝛼𝑚3(𝑚+𝑥)²−19𝜏𝑥𝛽9(𝑐+𝑥). (3) Referring [8], we know that if 𝜕𝑙₁,𝑙₂,𝑙₃𝜕𝜇₁,𝜇₂,𝜇_3𝜇=0≠0,𝜂₁𝜂₂𝜂₃𝜂₄≠0, (38) that is, 2𝜌≠𝑥²(𝑚+𝑥)3𝑚(𝑐+𝑥),𝜌≠19𝑥²(𝑚+𝑥)6𝑚(𝑐+𝑥), (39) for the unfolding normal form (36) then there exist the following results.

Theorem 5. Let (4), (9), and (39) hold. For the parameters 𝜇₁,𝜇₂, and 𝜇₃ are sufficiently small,(i)system (36) undergoes a transcritical bifurcation at the origin on the curve 𝑇=𝜇₁,𝜇₂,𝜇₃∶𝜇₃=0, (40)(ii)system (36) undergoes a Hopf bifurcation at the origin on the curve 𝐻₁=𝜇₁,𝜇₂,𝜇₃∶2𝛼𝑥𝛽(𝑚+𝑥)𝜇₁=19𝑐𝜏18(𝑐+𝑥)²𝜇₃−4𝜏3𝜇₂+3𝑐𝜇₃6(𝑐+𝑥)²𝜇₂+𝑐𝜇₃,𝜇₂𝑐>−6(𝑐+𝑥)²𝜇₃, (41)(iii)system (36) undergoes a Hopf-bifurcation at the nontrivial equilibrium point on the curve 𝐻₂=𝜇₁,𝜇₂,𝜇₃∶𝑙₃=𝜂₄𝜂₁−𝜂₁𝜂₃𝑙₁−𝜂₁𝑙₂𝑙₁,𝜂₁𝜂₃𝑙₁−𝜂₁𝑙₂>0, (42)(iv)system (36) undergoes a Bogdanov-Takens bifurcation at the origin on the curve 𝐵=𝜇₁,𝜇₂,𝜇₃∶𝜇₃=0,𝜇₂𝑐=−6(𝑐+𝑥)²𝜇₃, (43)(v)system (36) undergoes a zero-Hopf bifurcation at the origin on the curve 𝐻₃=𝜇₁,𝜇₂,𝜇₃∶𝜇₃=0,𝜇₁=2𝜏𝛽(𝑚+𝑥)3𝛼𝑥𝜇₂,𝜇₂>0. (44)

4. Hopf Bifurcation

The Jacobian matrix of system (2) at 𝐸^∗ takes the following form: 𝐽𝐸^∗=𝑥^∗−1+𝛼𝑦^∗(𝑚+𝑥^∗)²+ℎ(𝑐+𝑥^∗)²−𝛼𝑥^∗𝑚+𝑥^∗𝜌𝛽𝑒^−𝜆𝜏−𝜌𝑒^−𝜆𝜏, (45) then, the characteristic equation is 𝐹(𝜆)=𝜆²+𝜆𝑥^∗𝛼1−𝛽(𝑚+𝑥^∗)−ℎ(𝑐+𝑥^∗)²+𝜌𝜆+𝜌𝑥^∗ℎ1−(𝑐+𝑥^∗)²𝑒^−𝜆𝜏=𝜆²+𝜆𝑥^∗√Δ𝑐+𝑥^∗−𝛼𝛽(𝑚+𝑥^∗)+𝜌𝜆+𝜌𝑥^∗√Δ𝑐+𝑥^∗𝑒^−𝜆𝜏. (46) When 𝜏=0, by Routh-Hurwitz criterion, all roots of (46) have negative real part if 𝜌>𝑥^∗𝛼𝛽(𝑚+𝑥^∗)−√Δ𝑐+𝑥^∗. (47) Using the formula in Theorem 2.4 of [13] or Lemma 2.2 of [16], we have the following theorem.

Theorem 6. Let (47) hold. Then, there exist 𝑤₀=√22𝜌²−𝐴²+𝜌²−𝐴²²+4𝐵²,𝜏_𝑘=1𝑤₀𝑤arccos𝐵−𝐴𝜌²₀𝐵²+𝜌²𝑤²₀+2𝑘𝜋𝑤₀,𝑘=0,1,… (48) such that the positive equilibrium 𝐸^∗ is asymptotically stable for 0⩽𝜏<𝜏₀, moreover, system (2) undergoes Hopf bifurcation at 𝐸^∗ for 𝜏=𝜏_𝑘, where 𝐴=𝜌+𝑥^∗(√Δ/(𝑐+𝑥^∗)−𝛼/𝛽(𝑚+𝑥^∗)) and 𝐵=𝜌𝑥^∗√Δ/(𝑐+𝑥^∗).

In the following, using the Hopf bifurcation theorem for a retarded differential system introduced by [17], the Hopf bifurcations at 𝐸^∗ for small delay and large delay are presented.

Using the same methods as the ones used in [13], for small delay, let 𝑒^−𝜆𝜏≃1−𝜆𝜏, together with (46), the following bifurcation results can be obtained.

Theorem 7. Let 0<𝜏<𝜏=1/𝜌, 𝜌<𝛼𝑥^∗/𝛽(𝑚+𝑥^∗), 𝜏_𝑠=𝐴/𝐵 and 𝑤_𝑠=𝐵/(1−𝜌𝜏_𝑠). Then, there exists 𝜖_𝑠>0 such that for each 0⩽𝜖<𝜖_𝑠, system (2) near 𝐸^∗ has a family or periodic solutions 𝛾_𝑠(𝜖) with period 𝑇_𝑠=𝑇_𝑠(𝜖) for 𝜏=𝜏(𝜖) such that 𝛾_𝑠(0)=𝐸^∗, 𝑇_𝑠(0)=2𝜋/𝑤_𝑠, and 𝜏(0)=𝜏_𝑠.

For large delay, by [13] we have the following results.

Theorem 8. Let (47) hold. Then there exists 𝜖₀>0 such that for each 0⩽𝜖<𝜖₀ system (2) near 𝐸^∗ has a family or periodic solutions 𝛾_𝑙(𝜖) with period 𝑇_𝑙=𝑇_𝑙(𝜖) for 𝜏=𝜏(𝜖) such that 𝛾_𝑙(0)=𝐸^∗, 𝑇_𝑙(0)=2𝜋/𝑤₀ and 𝜏(0)=𝜏₀.

Remark. Because the proofs of Theorems 7 and 8 are the same as the proofs of Theorems 3.4 and 3.5 in [13], we omit them here.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This paper is supported by NSFC (11226142), the Foundation of Henan Educational Committee (2012A110012), the Foundation of Henan Normal University (2011QK04, 2012PL03), and the Scientific Research Foundation for Ph.D. of Henan Normal University (no. 1001).