Bifurcation Analysis for a Predator-Prey Model with Time Delay and Delay-Dependent Parameters

Abstract

A class of stage-structured predator-prey model with time delay and delay-dependent parameters is considered. Its linear stability is investigated and Hopf bifurcation is demonstrated. Using normal form theory and center manifold theory, some explicit formulae for determining the stability and the direction of the Hopf bifurcation periodic solutions bifurcating from Hopf bifurcations are obtained. Finally, numerical simulations are performed to verify the analytical results.

1. Introduction

Over the past decade, a great many predator-prey models have been developed to describe the interaction between predator and prey. Their dynamical phenomena have been extensively studied because of the wide application in the field of biomathematics. In particular, the appearance of a cycle bifurcating from the equilibrium of an ordinary or a delayed predator-prey model with a single parameter, which is known as a Hopf bifurcation, has attracted much attention due to its theoretical and practical significance [1–5]. But most of the research literature on these models are only connected with parameters which are independent of time delay; thus, the corresponding characteristic equations are easy to deal with. While in most applications of delay differential equations in population dynamics, the need of incorporation of a time delay is often the result of existence of some stage structure [6–8]. Indeed, every population goes through some distinct life stages [9, 10]. Since the through-stage survival rate is often a function of time delay, it is easy to conceive that these models will inevitably involve some delay-dependent parameters. Thus, the corresponding characteristic equations dependent on the delay 𝜏 become more complicated. In view of the fact that it is often difficult to analytically study models with delay-dependent parameters even if only a single discrete delay is present, we resort to the help of computer programs.

In 2008, Wang et al. [11] introduced and investigated the following predator-prey interaction model with time delay and delay-dependent parameters: 𝑥̇(𝑡)=𝑥(𝑡)[𝑎−𝑏𝑥(𝑡)]−𝑐𝑥²(𝑡)𝑦(𝑡)1+𝜛𝑥²(𝑡),𝑦̇(𝑡)=𝑐𝑥²(𝑡−𝜏)𝑦(𝑡−𝜏)1+𝜛𝑥²()𝑒𝑡−𝜏^−𝑑𝜏−𝑑𝑦(𝑡), (1.1) where 𝑥(𝑡) and 𝑦(𝑡) stand for prey and predator density at time 𝑡, respectively. 𝑎,𝑏,𝑐,𝜛 are real positive parameters and the time delay 𝜏 is a positive constant. Wang et al. [11] obtained the conditions that guarantee the system asymptotically stable and permanent. For more knowledge about the model, one can see [11].

It is well known that time delays which occur in the interaction between predator-prey will affect the stability of a model by creating instability, oscillation, and chaos phenomena. Based on the discussion above, the main purpose of this paper is to investigate the stability and the properties of Hopf bifurcation of the model (1.1) which involves some delay-dependent parameters. Recently, there are few papers on the topic that involves some delay-dependent parameters, for example, Liu and Zhang [12] investigated the stability and Hopf bifurcation of the following SIS model with nonlinear birth rate: 𝐼̇(𝑡)=𝛽(𝑁(𝑡)−𝐼(𝑡))𝐼(𝑡)𝑁(𝑡)−𝐼(𝑡),𝑁̇(𝑡)=()𝑑+𝜀+𝛾𝑃𝑁𝑡−𝜏1+𝑞𝑁³()𝑒𝑡−𝜏^−𝑑1𝜏−𝑑𝑁(𝑡)−𝜀𝐼(𝑡). (1.2) Jiang and Wei [13] studied the stability and Hopf bifurcation of the following SIR model: 𝑆̇(𝑡)=𝜇−𝜇𝑆(𝑡)−𝜙𝐼(𝑡)𝑆(𝑡)(𝑡)()𝑒1+𝐼+𝛾𝐼𝑡−𝜏^−𝜇𝜏,𝐼̇(𝑡)=𝜙𝐼(𝑡)𝑆(𝑡)(𝑡)−𝐼(𝑡).1+𝐼𝜇+𝛾 (1.3) It worth pointing out that Liu and Zhang [12] investigated the Hopf bifurcation of system (1.2) by choosing 𝑝 (not delay 𝜏) as the bifurcation parameters and Jiang and Wei [13] studied the Hopf bifurcation of system (1.3) by choosing 𝜙 (not delay 𝜏) as the bifurcation parameters. In this paper, we will investigate the Hopf bifurcation by regarding the delay 𝜏 as the bifurcation parameter which is different from the papers [12, 13]. To the best of our knowledge, it is the first time to deal with the stability and Hopf bifurcation of system (1.1).

This paper is organized as follows. In Section 2, the stability of the equilibrium and the existence of Hopf bifurcation at the equilibrium are studied. In Section 3, the direction of Hopf bifurcation and the stability and periodic of bifurcating periodic solutions on the center manifold are determined. In Section 4, numerical simulations are carried out to illustrate the validity of the main results. Some main conclusions are drawn in Section 5.

2. Stability of the Equilibrium and Local Hopf Bifurcations

Throughout this paper, we assume that the following condition(H1)𝑐𝑒^−𝑑𝜏>𝑑𝜛,𝑎²(𝑐𝑒^−𝑑𝜏−𝑑𝜛)>𝑏²𝑑 holds.

The hypothesis (H1) implies that system (1.1) has a unique positive equilibrium 𝐸_∗(𝑥^∗,𝑦^∗), where 𝑥^∗=𝑑𝑐𝑒^−𝑑𝜏−𝑑𝜛,𝑦^∗=𝑎−𝑏𝑥^∗1+𝜛𝑥^∗2𝑐𝑥^∗. (2.1) The linearized system of (1.1) around 𝐸_∗(𝑥^∗,𝑦^∗) takes the form 𝑑𝑥(𝑡)𝑑𝑡=𝑎₁𝑥(𝑡)−𝑏₁𝑦(𝑡),𝑑𝑦(𝑡)(𝑡)𝑑𝑡=−𝑑𝑦+𝑐₁𝑥(𝑡−𝜏)+𝑑₁𝑦(𝑡−𝜏), (2.2) where 𝑎₁=𝑎−2𝑏𝑥^∗−2𝑐𝑥^∗𝑦^∗1+𝜛𝑥^∗2+2𝑥^∗3𝑦^∗𝜛(1+𝜛𝑥^∗)²,𝑏₁=−𝑐𝑥^∗21+𝜛𝑥^∗2,𝑐₁=𝑒^−𝑑𝜏2𝑐𝑥^∗𝑦^∗1+𝜛𝑥^∗2−2𝑥^∗3𝑦^∗𝜛(1+𝜛𝑥^∗)²,𝑑₁=𝑒^−𝑑𝜏𝑐𝑥^∗21+𝜛𝑥^∗2. (2.3) The associated characteristic equation of (2.2) is 𝑃(𝜆,𝜏)+𝑄(𝜆,𝜏)𝑒^−𝜆𝜏=0, (2.4) where 𝑃(𝜆,𝜏)=𝜆²+𝑑−𝑎₁𝜆−𝑎₁𝑄()𝑑𝑑,𝜆,𝜏=−₁𝜆−𝑎₁𝑑₁−𝑏₁𝑐₁. (2.5) When 𝜏=0, then (2.4) becomes 𝜆²+𝑑−𝑎₁−𝑑⁰₁𝜆+𝑏₁𝑐⁰₁−𝑎₁𝑑−𝑎₁𝑑⁰₁=0, (2.6) where 𝑐⁰₁=2𝑐𝑥^∗𝑦^∗1+𝜛𝑥^∗2−2𝑥^∗3𝑦^∗𝜛(1+𝜛𝑥^∗)²,𝑑⁰₁=𝑐𝑥^∗1+𝜛𝑥^∗2. (2.7) It is easy to obtain the following result.

Lemma 2.1. If the condition
   (H2)𝑑−𝑎₁−𝑑⁰₁>0,b₁𝑐⁰₁−𝑎₁𝑑−𝑎₁𝑑⁰₁>0,
holds, then the positive equilibrium 𝐸_∗(𝑥^∗,𝑦^∗) of system (1.1) is asymptotically stable.

In the following, one investigates the existence of purely imaginary roots 𝜆=𝑖𝜔(𝜔>0) of (2.4). Equation (2.4) takes the form of a second-degree exponential polynomial in 𝜆, which some of the coefficients of P and Q depend on 𝜏. Beretta and Kuang [14] established a geometrical criterion which gives the existence of purely imaginary roots of a characteristic equation with delay-dependent coefficients. In order to apply the criterion due to Beretta and Kuang [14], one needs to verify the following properties for all 𝜏∈[0,𝜏_max), where 𝜏_max is the maximum value which 𝐸_∗(𝑥^∗,𝑦^∗) exists.(a)𝑃(0,𝜏)+𝑄(0,𝜏)≠0; (b)𝑃(𝑖𝜔,𝜏)+𝑄(𝑖𝜔,𝜏)≠0; (c)limsup{|𝑄(𝜆,𝜏)/𝑃(𝜆,𝜏)|∶|𝜆|→∞,Re𝜆≥0}<1; (d)𝐹(𝜔,𝜏)=|𝑃(𝑖𝜔,𝜏)|²−|𝑄(𝑖𝜔,𝜏)|² has a finite number of zeros;(e)Each positive root 𝜔(𝜏) of 𝐹(𝜔,𝜏)=0 is continuous and differentiable in 𝜏 whenever it exists. Here, 𝑃(𝜆,𝜏) and 𝑄(𝜆,𝜏) are defined as in (2.5), respectively.

Let 𝜏∈[0,𝜏_max). It is easy to see that 𝑃(0,𝜏)+𝑄(0,𝜏)=−𝑎₁𝑑+𝑎₁𝑑₁+𝑏₁𝑐₁≠0, (2.8) which implies that (a) is satisfied, and (b) 𝑃(𝑖𝜔,𝜏)+𝑄(𝑖𝜔,𝜏)=−𝜔²+𝑖𝜔𝑑−𝑎₁−𝑎₁𝑑−𝑖𝜔𝑑₁+𝑎₁𝑑₁+𝑏₁𝑐₁=−𝜔²−𝑎₁𝑑+𝑎₁𝑑₁+𝑏₁𝑐₁+𝑖𝜔𝑑−𝑎₁−𝑑₁≠0. (2.9) From (2.4), one has lim_{|𝜆|→+∞}𝑄(𝜆,𝜏)𝑃(𝜆,𝜏)=lim_{|𝜆|→+∞}−𝑑₁𝜆+𝑎₁𝑑₁−𝑏₁𝑐₁𝜆²+𝑑−𝑎₁𝜆−𝑎₁𝑑=0. (2.10) Therefore, (c) follows. Let 𝐹 be defined as in (d). From |𝑃(𝑖𝜔,𝜏)|²=𝜔²+𝑎₁𝑑²+𝑑−𝑎₁²𝜔²=𝜔⁴+𝑑²+𝑎²₁𝜔²+𝑎²₁𝑑²,|𝑄(𝑖𝜔,𝜏)|²=𝑑²₁𝜔²+𝑎₁𝑑₁+𝑏₁𝑐₁², (2.11) one obtain 𝐹(𝜔,𝜏)=𝜔⁴+𝑑²+𝑎²₁−𝑑²₁𝜔²+𝑎₁𝑑²−𝑎₁𝑑₁+𝑏₁𝑐₁². (2.12) Obviously, property (d) is satisfied, and by implicit function theorem, (e) is also satisfied.

Now let 𝜆=𝑖𝜔(𝜔>0) be a root of (2.4). Substituting it into (2.4) and separating the real and imaginary parts yields 𝑎₁𝑑₁+𝑏₁𝑐₁cos𝜔𝜏−𝑑₁𝜔sin𝜔𝜏=𝜔²+𝑎₁𝑑𝑑,₁𝜔cos𝜔𝜏+𝑎₁𝑑₁+𝑏₁𝑐₁sin𝜔𝜏=𝑑−𝑎₁𝜔. (2.13) From (2.13), it follows that 𝜔sin𝜔𝜏=−²+𝑎₁𝑑𝑑₁𝜔−𝑑−𝑎₁𝜔𝑎₁𝑑₁+𝑏₁𝑐₁𝑑²₁𝜔²+𝑎₁𝑑₁+𝑏₁𝑐₁²,𝜔cos𝜔𝜏=²+𝑎₁𝑑𝑎₁𝑑₁+𝑏₁𝑐₁+𝑑−𝑎₁𝜔𝑑₁𝜔𝑑²₁𝜔²+𝑎₁𝑑₁+𝑏₁𝑐₁². (2.14) By the definitions of 𝑃 and 𝑄 as in (2.5), respectively, and applying the property (a), then (2.14) can be written as 𝑃()sin𝜔𝜏=Im𝑖𝜔,𝜏𝑄(),𝑃()𝑖𝜔,𝜏cos𝜔𝜏=−Re𝑖𝜔,𝜏𝑄(),𝑖𝜔,𝜏 (2.15) which yields |𝑃(𝑖𝜔,𝜏)|²=|𝑄(𝑖𝜔,𝜏)|². Assume that 𝐼∈𝑅₊₀ is the set where 𝜔(𝜏) is a positive root of 𝐹(𝜔,𝜏)=|𝑃(𝑖𝜔,𝜏)|²−|𝑄(𝑖𝜔,𝜏)|², (2.16) and for 𝜏∈𝐼,𝜔(𝜏) is not definite. Then for all 𝜏 in 𝐼,𝜔(𝜏) satisfied 𝐹(𝜔,𝜏)=0. The polynomial function 𝐹 can be written as 𝐹(𝜔,𝜏)=ℎ𝜔²,𝜏, (2.17) where ℎ is a second degree polynomial, defined by ℎ(𝑧,𝜏)=𝑧²+𝑑²+𝑎²₁−𝑑²₁𝑧+𝑎₁𝑑²−𝑎₁𝑑₁+𝑏₁𝑐₁². (2.18) It is easy to see that ℎ(𝑧,𝜏)=𝑧²+𝑑²+𝑎²₁−𝑑²₁𝑧+𝑎₁𝑑²−𝑎₁𝑑₁−𝑏₁𝑐₁²=0 (2.19) has only one positive real root if the following condition (H3) holds:(H3)𝑎₁𝑑²<(𝑎₁𝑑₁+𝑏₁𝑐₁)².

One denotes this positive real root by 𝑧₊. Hence, (2.17) has only one positive real root 𝜔=√𝑧₊. Since the critical value of 𝜏 and 𝜔(𝜏) are impossible to solve explicitly, so one will use the procedure described in Beretta and Kuang [14]. According to this procedure, one defines 𝜃(𝜏)∈[0,2𝜋) such that sin𝜃(𝜏) and cos𝜃(𝜏) are given by the righthand sides of (2.14), respectively, with 𝜃(𝜏) given by (2.19). This define 𝜃(𝜏) in a form suitable for numerical evaluation using standard software. And the relation between the argument 𝜃 and 𝜔𝜏 in (2.18) for 𝜏>0 must be 𝜔𝜏=𝜃+2𝑛𝜋,  𝑛=1,2,….

Hence, one can define the maps: 𝜏_𝑛∶𝐼→𝑅₊₀ given by 𝜏_𝑛(𝜏)𝜃(𝜏)∶=+2𝑛𝜋𝜔(𝜏),𝜏_𝑛>0,𝑛=0,1,2,…, (2.20) where a positive root 𝜔(𝜏) of 𝐹(𝜔,𝜏)=0 exists in I. Let us introduce the functions 𝑆_𝑛(𝜏)∶𝐼→𝑅,𝑆_𝑛(𝜏)𝜃(𝜏)=𝜏−+2𝑛𝜋𝜔(𝜏),𝑛=0,1,2,…, (2.21) which are continuous and differentiable in 𝜏. Thus, one gives the following theorem which is due to Beretta and Kuang [14].

Theorem 2.2. Assume that 𝜔(𝜏) is a positive root of (2.4) defined for 𝜏∈𝐼,𝐼⊆𝑅₊₀, and at some 𝜏₀∈𝐼,𝑆_𝑛(𝜏₀)=0 for some 𝑛∈𝑁₀. Then, a pair of simple conjugate pure imaginary roots 𝜆=±𝑖𝜔 exists at 𝜏=𝜏₀ which crosses the imaginary axis from left to right if 𝛿(𝜏₀)>0 and crosses the imaginary axis from right to left if 𝛿(𝜏₀)<0, where 𝛿(𝜏₀)=sign[𝐹^′_𝜔(𝜔𝜏₀,𝜏₀)]sign[(𝑑𝑆_𝑛(𝜏))/𝑑𝜏|_𝜏=𝜏0].

Applying Lemma 2.1 and the Hopf bifurcation theorem for functional differential equation [5], we can conclude the existence of a Hopf bifurcation as stated in the following theorem.

Theorem 2.3. For system (1.1), if (H1)–(H3) hold, then there exists s 𝜏₀∈𝐼 such that the positive equilibrium 𝐸_∗(𝑥^∗,𝑦^∗) is asymptotically stable for 0≤𝜏<𝜏₀ and becomes unstable for 𝜏 staying in some right neighborhood of 𝜏₀, with a Hopf bifurcation occurring when 𝜏=𝜏₀.

3. Direction and Stability of the Hopf Bifurcation

In the previous section, we obtained some conditions which guarantee that the stage-structured predator-prey model with time delay undergoes the Hopf bifurcation at some values of 𝜏=𝜏₀. In this section, we will derive the explicit formulae determining the direction, stability, and period of these periodic solutions bifurcating from the positive equilibrium 𝐸_∗(𝑥^∗,𝑦^∗) at these critical value of 𝜏, by using techniques from normal form and center manifold theory [15]. Throughout this section, we always assume that system (1.1) undergoes Hopf bifurcation at the positive equilibrium 𝐸_∗(𝑥^∗,𝑦^∗) for 𝜏=𝜏₀, and then ±𝑖𝜔₀ is corresponding purely imaginary roots of the characteristic equation at the equilibrium 𝐸_∗(𝑥^∗,𝑦^∗).

For convenience, let 𝜏=𝜏₀+𝜇,𝜇∈𝑅. Then 𝜇=0 is the Hopf bifurcation value of (1.1). Thus, one will study Hopf bifurcation of small amplitude periodic solutions of (1.1) from the positive equilibrium point 𝐸_∗(𝑥^∗,𝑦^∗) for 𝜇 close to 0.

Let 𝑢₁(𝑡)=𝑥(𝑡)−𝑥^∗, 𝑢₂(𝑡)=𝑦(𝑡)−𝑦^∗, 𝑥_𝑖(𝑡)=𝑢_𝑖(𝜏𝑡), (𝑖=1,2),𝜏=𝜏₀+𝜇, then system (1.1) can be transformed into an functional differential equation (FDE) in (𝐶=𝐶(−1,0],𝑅²) as𝑑𝑢𝑑𝑡=𝐿_𝜇𝑢_𝑡+𝑓𝜇,𝑢_𝑡, (3.1) where 𝑢(𝑡)=(𝑥₁(𝑡),𝑥₂(𝑡))^𝑇∈𝑅² and 𝐿_𝜇∶𝐶→𝑅,𝑓∶𝑅×𝐶→𝑅 are given, respectively, by 𝐿_𝜇𝜙=𝜏₀+𝜇𝐵𝜙(0)+𝜏₀+𝜇𝐺𝜙(−1), (3.2) where 𝑎𝐵=¹−𝑏₁0−𝑑𝑐,𝐶=00₁𝑑₁,𝑎₁=𝑎−2𝑏𝑥^∗−2𝑐𝑥^∗𝑦^∗1+𝜛𝑥^∗2+2𝑥^∗3𝑦^∗𝜛(1+𝜛𝑥^∗)²,𝑏₁=−𝑐𝑥^∗1+𝜛𝑥^∗2,𝑐₁=𝑒^−𝑑𝜏2𝑐𝑥^∗𝑦^∗1+𝜛𝑥^∗2−2𝑥^∗3𝑦^∗𝜛(1+𝜛𝑥^∗)²,𝑑₁=𝑒^−𝑑𝜏𝑐𝑥^∗1+𝜛𝑥^∗2𝑓=𝑓𝜇,𝜙¹𝑓𝜇,𝜙₂,𝜇,𝜙 (3.3) where 𝑓₁𝜇,𝜙=𝜏₀+𝜇𝑚₁𝜙²₁(0)+𝑚₂𝜙₁(0)𝜙₂(0)+𝑚₃𝜙³₁(0)+𝑚₄𝜙²₁(0)𝜙₂(0),𝑓+h.o.t.₂𝜇,𝜙=𝜏₀+𝜇𝑛₁𝜙²₁(−1)+𝑛₂𝜙₁(−1)𝜙₂(−1)+𝑛₃𝜙³₁(−1)+𝑚₄𝜙²₁(−1)𝜙₂(−1)+h.o.t., (3.4) where 𝑚₁=4𝑥^∗2𝑦^∗𝜛(1+𝜛𝑥^∗)²−𝑐𝑦^∗1+𝜛𝑥^∗2−𝑐𝑥^∗2𝑦^∗4𝜛²𝑥^∗−𝜛𝑥^∗2−𝜛(1+𝜛𝑥^∗)⁴,𝑚₂=2𝑐𝜛𝑥^∗31+𝜛𝑥^∗22−2𝑐𝑥^∗1+𝜛𝑥^∗2,𝑚₃=2𝑐𝜛𝑥^∗𝑦^∗1+𝜛𝑥^∗22−𝑐𝑥^∗2𝑦^∗𝑔6−2𝑐𝑥^∗𝑦^∗4𝜛²𝑥^∗−𝜛²𝑥^∗2−𝜛1+𝜛𝑥^∗24,𝑚₄=4𝑐𝜛𝑥^∗21+𝜛𝑥^∗22−𝑐1+𝜛𝑥^∗−𝑐𝑥^∗24𝜛²𝑥^∗−𝜛²𝑥^∗2−𝜛1+𝜛𝑥^∗24,𝑛₁=𝑒^−𝑑𝜏0𝑐𝑦^∗1+𝜛𝑥^∗2+𝑐𝑥^∗2𝑦^∗4𝜛²𝑥^∗−𝜛𝑥^∗2−𝜛(1+𝜛𝑥^∗)⁴−4𝑥^∗2𝑦^∗𝜛(1+𝜛𝑥^∗)²,𝑛₂=𝑒^−𝑑𝜏02𝑐𝑥^∗1+𝜛𝑥^∗2−2𝑐𝜛𝑥^∗31+𝜛𝑥^∗22,𝑛₃=𝑒^−𝑑𝜏0𝑐𝑥^∗2𝑦^∗𝑔₁6+2𝑐𝑥^∗𝑦^∗4𝜛²𝑥^∗−𝜛²𝑥^∗2−𝜛1+𝜛𝑥^∗24−2𝑐𝜛𝑥^∗𝑦^∗1+𝜛𝑥^∗22,𝑛₄=𝑒^−𝑑𝜏0𝑐1+𝜛𝑥^∗+𝑐𝑥^∗24𝜛²𝑥^∗−𝜛²𝑥^∗2−𝜛1+𝜛𝑥^∗24−4𝑐𝜛𝑥^∗21+𝜛𝑥^∗22, (3.5) where 𝑔=18𝜛²𝑥^∗1+𝜛²𝑥^∗23𝜛²𝑥^∗4−4𝜛²𝑥^∗3+2𝜛𝑥^∗2−4𝜛𝑥^∗+11+𝜛𝑥^∗28−8𝜛²𝑥^∗+4𝑥^∗3−3𝜛𝑥^∗2−4𝜛1+𝜛𝑥^∗24. (3.6)

Clearly, 𝐿_𝜇 is a linear continuous operator from 𝐶 to 𝑅². By the Riesz representation theorem, there exists a matrix function with bounded variation components 𝜂(𝜃,𝜇),𝜃∈[−1,0] such that 𝐿_𝜇𝜙=⁰₋₁𝑑𝜂𝜃,𝜇𝜙(𝜃),for𝜙∈𝐶. (3.7) In fact, we can choose 𝜂𝜃,𝜇=𝜏₀𝑎+𝜇¹−𝑏₁𝛿(𝜃)−𝜏0−𝑑⁰𝑐+𝜇00₁𝑑₁𝛿(𝜃+1), (3.8) where 𝛿 is the Dirac delta function.

For 𝜙∈𝐶([−1,0],𝑅²), define 𝐴𝜇(𝜃)𝜙=𝑑𝜙∫𝑑𝜃,−1≤𝜃<0,₀₋₁𝜙(𝑠)𝑅𝜇𝑓𝑑𝜂𝑠,𝜇,𝜃=0,𝜙=0,−1≤𝜃<0,𝜇,𝜙,𝜃=0. (3.9) Then (1.1) is equivalent to the abstract differential equation 𝑢̇_𝑡=𝐴𝜇𝑢_𝑡+𝑅𝜇𝑢_𝑡, (3.10) where 𝑢=(𝑥₁,𝑥₂)^𝑇, 𝑢_𝑡(𝜃)=𝑢(𝑡+𝜃), 𝜃∈[−1,0].

For 𝜓∈𝐶([0,1],(𝑅²)^∗), define 𝐴^∗𝜓(𝑠)=−𝑑𝜓(𝑠)(],∫𝑑𝑠,𝑠∈0,1₀₋₁𝑑𝜂^𝑇(𝑡,0)𝜓(−𝑡),𝑠=0. (3.11)

For 𝜙∈𝐶([−1,0],𝑅²) and 𝜓∈𝐶([0,1],(𝑅²)^∗), define the bilinear form 𝜓,𝜙=𝜓(0)𝜙(0)−⁰₋₁^𝜃_𝜉=0𝜓(𝜉−𝜃)𝑑𝜂(𝜃)𝜙(𝜉)𝑑𝜉, (3.12) where 𝜂(𝜃)=𝜂(𝜃,0). We have the following result on the relation between the operators 𝐴=𝐴(0) and 𝐴^∗.

Lemma 3.1. 𝐴=𝐴(0) and 𝐴^∗ are adjoint operators.

The proof is easy from (3.12), so we omit it.

By the discussions in Section 2, we know that ±𝑖𝜔₀𝜏₀ are eigenvalues of 𝐴(0), and they are also eigenvalues of 𝐴^∗ corresponding to 𝑖𝜔₀𝜏₀ and −𝑖𝜔₀𝜏₀, respectively. We have the following result.

Lemma 3.2. The vector 𝑞(𝜃)=1,𝛾^𝑇𝑒^{𝑖𝜔0𝜏0𝜃},𝜃∈[−1,0], (3.13) where 𝛾=𝑖𝜔₀−𝑎₁𝑏₁ (3.14) is the eigenvector of 𝐴(0) corresponding to the eigenvalue 𝑖𝜔₀𝜏₀, and 𝑞^∗(𝑠)=𝐷1,𝛾^∗𝑒^{𝑖𝜔0𝜏0𝑠},𝑠∈[0,1], (3.15) where 𝛾^∗=−𝑖𝜔₀𝑎₁𝑐₁𝑒^{−𝑖𝜔𝜏0} (3.16) is the eigenvector of 𝐴^∗ corresponding to the eigenvalue −𝑖𝜔₀𝜏₀; moreover, ⟨𝑞^∗(𝑠),𝑞(𝜃)⟩=1, where 𝐷=1+𝛾𝛾^∗+𝑐₁𝛾^∗𝑒^{𝑖𝜔0𝜏0}+𝑑₁𝛾𝛾^∗𝑒^{𝑖𝜔0𝜏0}. (3.17)

Proof. Let 𝑞(𝜃) be the eigenvector of 𝐴(0) corresponding to the eigenvalue 𝑖𝜔₀𝜏₀ and 𝑞^∗(𝑠) be the eigenvector of 𝐴^∗ corresponding to the eigenvalue −𝑖𝜔₀𝜏₀, namely, 𝐴(0)𝑞(𝜃)=𝑖𝜔₀𝜏₀𝑞(𝜃) and 𝐴^∗𝑞^∗𝑇(𝑠)=−𝑖𝜔₀𝜏₀𝑞^∗𝑇(𝑠). From the definitions of 𝐴(0) and 𝐴^∗, we have 𝐴(0)𝑞(𝜃)=𝑑𝑞(𝜃)/𝑑𝜃 and A^∗𝑞^∗𝑇(𝑠)=−𝑑𝑞^∗𝑇(𝑠)/𝑑𝑠. Thus, 𝑞(𝜃)=𝑞(0)𝑒^{𝑖𝜔0𝜏0𝜃} and 𝑞^∗(𝑠)=𝑞^∗(0)𝑒^{𝑖𝜔0𝜏0𝑠}. In addition, 𝑑𝜂(⁰₋₁𝜃)𝑞(𝜃)=𝜏₀𝐵𝑞(0)+𝜏₀𝐺𝑞(−1)=𝜏₀𝐴(0)𝑞(0)=𝑖𝜔₀𝜏₀𝑞(0). (3.18) That is, 𝑖𝜔₀−𝑎₁−𝑏₁−𝑐₁𝑒^{−𝑖𝜔0𝜏0}𝑖𝜔₀+𝑑−𝑑₁𝑒^{−𝑖𝜔0𝜏0}𝑞(0)=00. (3.19) Therefore, we can easily obtain 𝛾=𝑖𝜔₀−𝑎₁𝑏₁. (3.20) And so 𝑞(0)=1,𝑖𝜔₀−𝑎₁𝑏₁^𝑇, (3.21) and hence 𝑞(𝜃)=1,𝑖𝜔₀−𝑎₁𝑏₁^𝑇𝑒^{𝑖𝜔0𝜏0𝜃}. (3.22) On the other hand, ⁰₋₁𝑞^∗(−𝑡)𝑑𝜂(𝑡)=𝜏₀𝐵^𝑇𝑞^∗𝑇(0)+𝜏₀𝐺^𝑇𝑞^∗𝑇(−1)=𝜏₀𝐴^∗𝑞^∗𝑇(0)=−𝑖𝜔₀𝜏₀𝑞^∗𝑇(0). (3.23) Namely, −𝑖𝜔₀−𝑎₁−𝑐₁𝑒^{−𝑖𝜔0𝜏0}𝑏₁−𝑖𝜔₀+𝑑−𝑑₁𝑒^{−𝑖𝜔0𝜏0}𝑞^∗𝑇(0)=00. (3.24) Therefore, we can easily obtain 𝛾^∗=−𝑖𝜔₀𝑎₁𝑐₁𝑒^{−𝑖𝜔0𝜏0}, (3.25) and so 𝑞^∗(0)=1,−𝑖𝜔₀𝑎₁𝑐₁𝑒^{−𝑖𝜔𝜏0}, (3.26) and hence 𝑞^∗(𝑠)=1,−𝑖𝜔₀𝑎₁𝑐₁𝑒^{−𝑖𝜔𝜏0}𝑒^{𝑖𝜔0𝜏0𝑠}. (3.27) In the sequel, one will verify that ⟨𝑞^∗(𝑠),𝑞(𝜃)⟩=1. In fact, from (3.12), we have 𝑞^∗(𝑠),𝑞(𝜃)=𝐷1,𝛾^∗1,𝛾^𝑇−⁰₋₁⁰_𝜉=0𝐷1,𝛾^∗𝑒^{−𝑖𝜔0(𝜉−𝜃)}𝑑𝜂(𝜃)1,𝛾^𝑇𝑒^{𝑖𝜔0𝜉}=𝑑𝜉𝐷1+𝛾𝛾^∗−⁰₋₁1,𝛾^∗𝜃𝑒^{𝑖𝜔0𝜃}𝑑𝜂(𝜃)1,𝛾^𝑇=𝐷1+𝛾𝛾^∗−1,𝛾^∗−𝜏₀𝐶𝑒^{−𝑖𝜔0𝜏0}1,𝛾^𝑇=𝐷1+𝛾𝛾^∗+𝑐₁𝛾^∗𝑒^{−𝑖𝜔0𝜏0}+𝑑₁𝛾𝛾^∗𝑒^{−𝑖𝜔0𝜏0}=1. (3.28)
Next, we use the same notations as those in Hassard et al. [15], and we first compute the coordinates to describe the center manifold 𝐶₀ at 𝜇=0. Let 𝑢_𝑡 be the solution of (1.1) when 𝜇=0.
Define 𝑧(𝑡)=𝑞^∗,𝑢_𝑡,𝑊(𝑡,𝜃)=𝑢_𝑡(𝜃)−2Re𝑧(𝑡)𝑞(𝜃), (3.29) on the center manifold 𝐶₀, and we have 𝑊(𝑡,𝜃)=𝑊(𝑧(𝑡),𝑧(𝑡),𝜃), (3.30) where 𝑊(𝑧(𝑡),𝑧(𝑡),𝜃)=𝑊(𝑧,𝑧)=𝑊₂₀𝑧²2+𝑊₁₁𝑧𝑧+𝑊₀₂𝑧²2+⋯, (3.31) and 𝑧 and 𝑧 are local coordinates for center manifold 𝐶₀ in the direction of 𝑞^∗ and 𝑞^∗. Noting that 𝑊 is also real if 𝑢_𝑡 is real, we consider only real solutions. For solutions 𝑢_𝑡∈𝐶₀ of (1.1),𝑧̇(𝑡)=𝑞^∗(𝑠),𝑢̇_𝑡=𝑞^∗(𝑠),𝐴(0)𝑢_𝑡+𝑅(0)𝑢_𝑡=𝑞^∗(𝑠),𝐴(0)𝑢_𝑡+𝑞^∗(𝑠),𝑅(0)𝑢_𝑡=𝐴^∗𝑞^∗(𝑠),𝑢_𝑡+𝑞^∗(0)𝑅(0)𝑢_𝑡−⁰₋₁^𝜃_𝜉=0𝑞^∗(𝜉−𝜃)𝑑𝜂(𝜃)𝐴(0)𝑅(0)𝑢_𝑡(𝜉)=𝑑𝜉𝑖𝜔₀𝑞^∗(𝑠),𝑢_𝑡+𝑞^∗(0)𝑓0,𝑢_𝑡(𝜃)_def=𝑖𝜔₀𝑧(𝑡)+𝑞^∗(0)𝑓₀(𝑧(𝑡),𝑧(𝑡)). (3.32) That is, 𝑧̇(𝑡)=𝑖𝜔₀𝑧+𝑔(𝑧,𝑧), (3.33) where 𝑔(𝑧,𝑧)=𝑔₂₀𝑧²2+𝑔₁₁𝑧𝑧+𝑔₀₂𝑧²2+𝑔₂₁𝑧²𝑧2+⋯. (3.34) Hence, we have 𝑔(𝑧,𝑧)=𝑞^∗(0)𝑓₀(𝑧,𝑧)=𝑓0,𝑢_𝑡=𝐷1,𝛾^∗𝑓₁0,𝑢_𝑡,𝑓₂0,𝑢_𝑡^𝑇, (3.35) where 𝑓₁0,𝑢_𝑡=𝜏₀𝑚₁𝑥²_𝑡(0)+𝑚₂𝑥_𝑡(0)𝑦_𝑡(0)+𝑚₃𝑥³_𝑡(0)+𝑚₄𝑥²_𝑡(0)𝑦_𝑡(0),𝑓+h.o.t.₂0,𝑢_𝑡=𝜏₀𝑛₁𝑥²_𝑡(−1)+𝑛₂𝑥_𝑡(−1)𝑦_𝑡(−1)+𝑛₃𝑥³_𝑡(−1)+𝑚₄𝑥²_𝑡(−1)𝑦_𝑡(−1)+h.o.t.. (3.36) Noticing 𝑢_𝑡(𝜃)=(𝑥_𝑡(𝜃),𝑦_𝑡(𝜃))^𝑇=𝑊(𝑡,𝜃)+𝑧𝑞(𝜃)+𝑧𝑞(𝜃) and 𝑞(𝜃)=(1,𝛾)^𝑇𝑒^{𝑖𝜔0𝜏0𝜃}, we have𝑥_𝑡(0)=𝑧+𝑧+𝑊⁽¹⁾₂₀(0)𝑧²2+𝑊⁽¹⁾₁₁(0)𝑧𝑧+𝑊⁽¹⁾₀₂(0)𝑧²2𝑦+⋯,_𝑡(0)=𝛾𝑧+𝛾𝑧+𝑊⁽²⁾₂₀(0)𝑧²2+𝑊⁽²⁾₁₁(0)𝑧𝑧+𝑊⁽²⁾₀₂(0)𝑧²2𝑥+⋯,_𝑡(−1)=𝑒^{−𝑖𝜔0𝜏0}𝑧+𝑒^{𝑖𝜔0𝜏0}𝑧+𝑊⁽¹⁾₂₀()𝑧−1²2+𝑊⁽¹⁾₁₁(−1)𝑧𝑧+𝑊⁽¹⁾₀₂(−1)𝑧²2𝑦+⋯,_𝑡(−1)=𝛾𝑒^{−𝑖𝜔0𝜏0}𝛽𝑧+𝛾𝑒^{𝑖𝜔0𝜏0}𝑧+𝑊⁽²⁾₂₀()𝑧−1²2+𝑊⁽²⁾₁₁(−1)𝑧𝑧+𝑊⁽²⁾₀₂(−1)𝑧²2+⋯. (3.37) From (3.34) and (3.35), we have𝑔(𝑧,𝑧)=𝑞^∗(0)𝑓₀(𝑧,𝑧)=𝐷𝑓₁0,𝑥_𝑡+𝛾^∗𝑓₂0,𝑥_𝑡=𝐷𝜏₀𝑚₁+𝑚₂𝛾+𝛾^∗𝑛₁+𝑛₂𝛾𝑧²+𝐷𝜏₀2𝑚₁+𝑚₂𝛾+𝛾^∗2𝑛₁+𝑛₂𝛾+𝛾^∗𝑧𝑧+𝐷𝜏₀𝑚₁+𝑚₂𝛾+𝛾^∗𝑛₁+𝑛₂𝛾𝑧²+𝐷𝜏₀𝑚₁𝑊⁽¹⁾₂₀(0)+2𝑊⁽¹⁾₁₁(0)+𝑚₂𝛾𝑊⁽¹⁾₂₀(0)2+𝑊⁽²⁾₂₀(0)+𝛾𝑊⁽¹⁾₁₁(0)+𝑊⁽²⁾₁₁(0)+3𝑚₃+𝑚₄𝛾+2𝛾+𝑛₁𝑊⁽¹⁾₂₀(−1)+2𝑊⁽¹⁾₁₁(−1)+𝑛₂𝛾𝑊⁽¹⁾₂₀(−1)2+𝑊⁽²⁾₂₀(−1)+𝛾𝑊⁽¹⁾₁₁(−1)+𝑊⁽²⁾₁₁(−1)+3𝑛₃+𝑛₄𝛾+2𝛾𝑧²𝑧+h.o.t, (3.38) and we obtain 𝑔₂₀=2𝐷𝜏₀𝑚₁+𝑚₂𝛾+𝛾^∗𝑛₁+𝑛₂𝛾,𝑔₁₁=𝐷𝜏₀2𝑚₁+𝑚₂𝛾+𝛾^∗+𝛾^∗2𝑛₁+𝑛₂𝛾+𝛾^∗,𝑔₀₂=2𝐷𝜏₀𝑚₁+𝑚₂𝛾+𝛾^∗𝑛₁+𝑛₂𝛾,𝑔₂₁=2𝐷𝜏₀𝑚₁𝑊⁽¹⁾₂₀(0)+2𝑊⁽¹⁾₁₁(0)+𝑚₂𝛾𝑊⁽¹⁾₂₀(0)2+𝑊⁽²⁾₂₀(0)+𝛾𝑊⁽¹⁾₁₁(0)+𝑊⁽²⁾₁₁(0)+3𝑚₃+𝑚₄𝛾+2𝛾+𝑛₁𝑊⁽¹⁾₂₀(−1)+2𝑊⁽¹⁾₁₁(−1)+𝑛₂𝛾𝑊⁽¹⁾₂₀(−1)2+𝑊⁽²⁾₂₀(−1)+𝛾𝑊⁽¹⁾₁₁(−1)+𝑊⁽²⁾₁₁(−1)+3𝑛₃+𝑛₄𝛾+2𝛾. (3.39) For unknown 𝑊⁽¹⁾₂₀(0),𝑊⁽¹⁾₂₀(−1),𝑊⁽¹⁾₁₁(0),𝑊⁽²⁾₁₁(−1),𝑊⁽²⁾₁₁(0),𝑊⁽²⁾₁₁(−1), (3.40) in 𝑔₂₁, we still need to compute them.
From (3.10) and (3.29), we have 𝑊^′=𝑊−2Re𝑞^∗(0)𝑓𝑞(𝜃),−1≤𝜃<0,𝑊−2Re𝑞^∗(0)𝑓𝑞(𝜃)+𝑓,𝜃=0._def=𝐴𝑊+𝐻(𝑧,𝑧,𝜃), (3.41) where 𝐻(𝑧,𝑧,𝜃)=𝐻₂₀(𝜃)𝑧²2+𝐻₁₁(𝜃)𝑧𝑧+𝐻₀₂(𝜃)𝑧²2+⋯. (3.42) Comparing the coefficients, we obtain 𝐴−2𝑖𝜔₀𝜏₀𝑊₂₀=−𝐻₂₀(𝜃),𝐴𝑊₁₁(𝜃)=−𝐻₁₁(𝜃), (3.43)(3.44) and we know that for 𝜃∈[−1,0), 𝐻(𝑧,𝑧,𝜃)=−𝑞^∗(0)𝑓₀𝑞(𝜃)−𝑞^∗(0)𝑓₀𝑞(𝜃)=−𝑔(𝑧,𝑧)𝑞(𝜃)−𝑔(𝑧,𝑧)𝑞(𝜃). (3.45) Comparing the coefficients of (3.42) with (3.45) gives that 𝐻₂₀(𝜃)=−𝑔₂₀𝑞(𝜃)−𝑔₀₂𝑞(𝜃),𝐻₁₁(𝜃)=−𝑔₁₁𝑞(𝜃)−𝑔₁₁𝑞(𝜃). (3.46)(3.47) From (3.43) and (3.46) and the definition of 𝐴, we get 𝑊̇₂₀(𝜃)=2𝑖𝜔₀𝜏₀𝑊₂₀(𝜃)+𝑔₂₀𝑞(𝜃)+𝑔₀₂𝑞(𝜃). (3.48) Noting that 𝑞(𝜃)=𝑞(0)𝑒^{𝑖𝜔0𝜏0𝜃}, we have 𝑊₂₀(𝜃)=𝑖𝑔₂₀𝜔₀𝜏₀𝑞(0)𝑒^{𝑖𝜔0𝜏0𝜃}+𝑖𝑔₀₂3𝜔₀𝜏₀𝑞(0)𝑒^{−𝑖𝜔0𝜏0𝜃}+𝐸₁𝑒^{2𝑖𝜔0𝜏0𝜃}, (3.49) where 𝐸₁=(𝐸₁⁽¹⁾,𝐸₁⁽²⁾)^𝑇 is a constant vector.
Similarly, from (3.44) and (3.47) and the definition of 𝐴, we have 𝑊̇₁₁(𝜃)=𝑔₁₁𝑞(𝜃)+𝑔₁₁𝑞(𝜃),𝑊₁₁(𝜃)=−𝑖𝑔₁₁𝜔₀𝜏₀𝑞(0)𝑒^{𝑖𝜔0𝜏0𝜃}+𝑖𝑔₁₁𝜔₀𝜏₀𝑞(0)𝑒^{−𝑖𝜔0𝜃}+𝐸₂, (3.50)(3.51) where 𝐸₂=(𝐸₂⁽¹⁾,𝐸₂⁽²⁾)^𝑇 is a constant vector.
In what follows, one will seek appropriate 𝐸₁, 𝐸₂ in (3.49) and (3.51), respectively. It follows from the definition of 𝐴 and (3.46) and (3.47) that 𝑑𝜂(⁰₋₁𝜃)𝑊₂₀(𝜃)=2𝑖𝜔₀𝜏₀𝑊₂₀(0)−𝐻₂₀(0),₀₋₁𝑑𝜂(𝜃)𝑊₁₁(𝜃)=−𝐻₁₁(0), (3.52)(3.53) where 𝜂(𝜃)=𝜂(0,𝜃).
From (3.43), we have 𝐻₂₀(0)=−𝑔₂₀𝑞(0)−𝑔₀₂𝑞(0)+𝜏₀𝑀₁,𝑀₂^𝑇, (3.54) where 𝑀₁=𝑚₁+𝑚₂𝑀𝛾,₂=𝑛₁+𝑛₂𝛾. (3.55) From (3.44), we have 𝐻₁₁(0)=−𝑔₁₁𝑞(0)−𝑔₁₁(0)𝑞(0)+𝜏₀𝑁₁,𝑁₂^𝑇, (3.56) where 𝑁₁=2𝑚₁+𝑚₂𝛾+𝛾,𝑁₂=2𝑛₁+𝑛₂𝛾+𝛾. (3.57) Noting that 𝑖𝜔₀𝜏₀𝐼−⁰₋₁𝑒^{𝑖𝜔0𝜏0𝜃}𝑑𝜂(𝜃)𝑞(0)=0,−𝑖𝜔₀𝜏₀𝐼−⁰₋₁𝑒^{−𝑖𝜔0𝜏0𝜃}𝑑𝜂(𝜃)𝑞(0)=0, (3.58) and substituting (3.49) and (3.54) into (3.52), we have 2𝑖𝜔₀𝜏₀𝐼−⁰₋₁𝑒^{2𝑖𝜔0𝜏0𝜃}𝑑𝜂(𝜃)𝐸₁=𝜏₀𝑀₁,𝑀₂^𝑇. (3.59) That is, 2𝑖𝜔₀𝜏₀𝐼−𝜏₀𝐵−𝜏₀𝐺𝑒^{−2𝑖𝜔0𝜏0}𝐸₁=𝜏₀𝑀₁,𝑀₂^𝑇, (3.60) then 2𝑖𝜔₀−𝑎₁𝑏₁𝑐₁𝑒^{−2𝑖𝜔0𝜏0}2𝑖𝜔₀+𝑑−𝑑₁𝑒^{−2𝑖𝜔0𝜏0}𝐸₁⁽¹⁾𝐸₁⁽²⁾=𝑚₁+𝑚₂𝛾𝑛₁+𝑛₂𝛾. (3.61) Hence, 𝐸₁⁽¹⁾=Δ₁₁Δ₁,𝐸₁⁽²⁾=Δ₁₂Δ₁, (3.62) where Δ₁=det2𝑖𝜔₀−𝑎₁𝑏₁𝑐₁𝑒^{−2𝑖𝜔0𝜏0}2𝑖𝜔₀+𝑑−𝑑₁𝑒^{−2𝑖𝜔0𝜏0},Δ₁₁𝑚=det¹+𝑚₂𝛾𝑏₁𝑛₁+𝑛₂𝛾2𝑖𝜔₀+𝑑−𝑑₁𝑒^{−2𝑖𝜔0𝜏0},Δ₁₂=det2𝑖𝜔₀−𝑎₁𝑚₁+𝑚₂𝛾𝑐₁𝑒^{−2𝑖𝜔0𝜏0}𝑛₁+𝑛₂𝛾. (3.63) Similarly, substituting (3.51) and (3.56) into (3.53), we have 𝑑𝜂(⁰₋₁𝜃)𝐸₂=𝜏₀𝑁₁,𝑁₂^𝑇. (3.64) Then, (𝐵+𝐺)𝐸₂=−𝑁₁,−𝑁₂^𝑇. (3.65) That is, 𝑎₁−𝑏₁𝑐₁𝐷₁𝐸−𝑑𝐸²⁽¹⁾₂⁽²⁾=−𝑁₁−𝑁₂. (3.66) Hence, 𝐸₂⁽¹⁾=Δ₂₁Δ₂,𝐸₂⁽²⁾=Δ₂₂Δ₂, (3.67) where Δ₂𝑎=det¹−𝑏₁𝑐₁𝑑₁,Δ−𝑑₂₁=det−2𝑚₁−𝑚₂𝛾+𝛾−𝑏₁−2𝑛₁−𝑛₂𝛾+𝛾𝑑₁,Δ−𝑑₂₂𝑎=det¹−2𝑚₁−𝑚₂𝛾+𝛾𝑐₁−2𝑛₁−𝑛₂𝛾+𝛾. (3.68) From (3.49) and (3.51), we can calculate 𝑔₂₁ and derive the following values: 𝑐₁(0)=𝑖2𝜔₀𝜏₀𝑔₂₀𝑔₁₁−2𝑔₁₁²−𝑔₀₂²3+𝑔₂₁2,𝜇₂𝑐=−Re₁(0)Re𝜆^′𝜏₀,𝛽₂=2Re𝑐₁(0),𝑇₂𝑐=−Im₁(0)+𝜇₂Im𝜆^′𝜏₀𝜔₀𝜏₀. (3.69) These formulae give a description of the Hopf bifurcation periodic solutions of (1.1) at 𝜏=𝜏₀ on the center manifold. From the discussion above, we have the following result.

Theorem 3.3. The periodic solution is supercritical (subcritical) if 𝜇₂>0 (𝜇₂<0). The bifurcating periodic solutions are orbitally asymptotically stable with asymptotical phase (unstable) if 𝛽₂<0   (𝛽₂>0). The periods of the bifurcating periodic solutions increase (decrease) if 𝑇₂>0(𝑇₂<0).

4. Numerical Examples

In this section, we present some numerical results to verify the analytical predictions obtained in the previous section. As an example, we consider the following special case of system (1.1) with the parameters 𝑎=0.5,𝑏=0.5,𝑐=2,𝜛=2,and𝑑=0.2. Then system (1.1) becomes 𝑥̇(𝑡)=𝑥(𝑡)[1−0.5𝑥(𝑡)]−2𝑥²(𝑡)𝑦(𝑡)1+3𝑥²(𝑡),𝑦̇(𝑡)=2𝑥²(𝑡−𝜏)𝑦(𝑡−𝜏)1+3𝑥²()e𝑡−𝜏^−0.2𝜏−0.2𝑦(𝑡), (4.1) which has a positive equilibrium 𝐸_∗(𝑥^∗,𝑦^∗)=(1.4928,1.3217). By some complicated computation by means of Matlab 7.0, we get only one critical values of the delay 𝜏₀≈1.7355,𝜆^′(𝜏₀)≈0.2035−0.5423𝑖. Thus, we derive 𝑐₁(0)≈−1.3122−5.0131𝑖,𝜇₂≈0.6177,𝛽₂≈−3.3326,𝑇₂≈9.3042. We obtain that the conditions indicated in Theorem 2.3 are satisfied. Furthermore, it follows that 𝜇₂>0 and 𝛽₂<0. Thus, the positive equilibrium 𝐸_∗(𝑥^∗,𝑦^∗) is stable when 𝜏<𝜏₀ which is illustrated by the computer simulations (see Figures 1(a)–1(d)). When 𝜏 passes through the critical value 𝜏₀, the positive equilibrium 𝐸_∗(𝑥^∗,𝑦^∗) loses its stability and a Hopf bifurcation occurs, that is, a family of periodic solutions bifurcations from the positive equilibrium 𝐸_∗(𝑥^∗,𝑦^∗). Since 𝜇₂>0 and 𝛽₂<0, the direction of the Hopf bifurcation is 𝜏>𝜏₀, and these bifurcating periodic solutions from 𝐸_∗(𝑥^∗,𝑦^∗) at 𝜏₀ are stable, which are depicted in Figures 2(a)–2(d).

(a)

(b)

(c)

(d)

Figure 1 

The time histories and phase portrait of system (4.1) with the following parameters: 𝑎=0.5,𝑏=0.5,𝑐=2,𝜛=2,𝑑=0.2, and 𝜏=1.7<𝜏0≈1.7355. The positive equilibrium 𝐸∗(𝑥∗,𝑦∗)=(0.4928,1.3217) is asymptotically stable. The initial value is (0.2, 2).

(a)

(b)

(c)

(d)

Figure 2 

The time histories and phase portrait of system (4.1) with the following parameters: 𝑎=0.5,𝑏=0.5,𝑐=2,𝜛=2,𝑑=0.2, and 𝜏=1.8>𝜏0≈1.7355. Hopf bifurcation occurs from the positive equilibrium 𝐸∗(𝑥∗,𝑦∗)=(0.5014,1.3108). The initial value is (0.2, 2).

5. Conclusions

In this paper, the main object is to investigate the local stability and Hopf bifurcation and also to study the stability of bifurcating periodic solutions and some formulae for determining the direction of Hopf bifurcation for a stage-structured predator-prey model with time delay and delay. By choosing the delay as a bifurcation parameter, It is shown that under certain condition, the positive equilibrium 𝐸_∗(𝑥^∗,𝑦^∗) of system (1.1) is asymptotically stable for all 𝜏∈[0,𝜏₀) and unstable for 𝜏>𝜏₀ and under another condition; when the delay 𝜏 increases, the equilibrium loses its stability and a sequence of Hopf bifurcations occur at the positive equilibrium 𝐸_∗(𝑥^∗,𝑦^∗), that is, a family of periodic orbits bifurcate from the positive equilibrium 𝐸_∗(𝑥^∗,𝑦^∗). At the same time, using the normal form theory and the center manifold theorem, the direction of Hopf bifurcation and the stability of the bifurcating periodic orbits are discussed. Finally, numerical simulations are carried out to validate the theorems obtained.

Acknowledgments

This work is supported by The National Natural Science Foundation of China (no. 11261010), Soft Science and Technology Program of Guizhou Province (no. 2011LKC2030), Natural Science and Technology Foundation of Guizhou Province (J (2012) 2100), Governor Foundation of Guizhou Province ((2012)53), and Doctoral Foundation of Guizhou University of Finance and Economics (2010).