Stability and Hopf Bifurcation Analysis for a Stage-Structured Predator-Prey Model with Discrete and Distributed Delays

Abstract

We propose a three-dimensional stage-structured predatory-prey model with discrete and distributed delays. By use of a new variable, the original three-dimensional system transforms into an equivalent four-dimensional system. Firstly, we study the existence and local stability of positive equilibrium of the new system. And, by choosing the time delay 𝜏 as a bifurcation parameter, we show that Hopf bifurcation may occur as the time delay 𝜏 passes through some critical values. Secondly, by use of normal form theory and central manifold argument, we establish the direction and stability of Hopf bifurcation. At last, some simple discussion is presented.

1. Introduction

Since the pioneering theoretical works by Lotka [1] and Volterra [2], there were a lot of authors who studied all kinds of predator-prey models modeled by ordinary differential equations (ODEs) [3–5]. To reflect that the dynamical behavior of the models depends on the past history of the system, it is often necessary to incorporate time delays into the models. Therefore, a more realistic predator-prey model should be described by delayed differential equations (DDEs) [6–26]. Some of them investigated discrete delays [6–20]; others were about distributed delays [21–24]; and both discrete and distributed delays were studied in [25]. In general, delay differential equations exhibit more complicated dynamics on stability, periodic structure, bifurcation, and so on [26].

In the natural world, many individuals have a life story that takes them through two stages, immature and mature. The predator only catches the mature prey, as the immature preys are protected by their eggshells or refuge. Some predator-prey models with stage structure were investigated in [27–33]. Motivated by [25, 27, 31] and the references cited therein, in the present paper, we will consider the following stage-structured predator-prey model with discrete and distributed delay: 𝑥^′₁(𝑡)=𝑟𝑥₂(𝑡)−𝑑₁𝑥₁(𝑡)−𝑟𝑒^−𝑑1𝜏𝑥₂(),𝑥𝑡−𝜏_′2(𝑡)=𝑟𝑒^−𝑑1𝜏𝑥₂(𝑡−𝜏)−𝑑₂𝑥₂(𝑡)−𝑏₁𝑥₂(𝑡)𝑦(𝑡),𝑦^′(𝑡)=𝑦(𝑡)−𝑑₃+𝑏₂^𝑡_−∞𝐹(𝑡−𝑠)𝑥₂(𝑠)d𝑠−𝛼𝑦(𝑡), (1) where 𝑥₁(𝑡), 𝑥₂(𝑡), and 𝑦(𝑡) can be interpreted as the population densities of the immature prey, mature prey, and predator at time 𝑡, respectively. 𝑟 denotes the birth rate of the prey population; 𝑑₁, 𝑑₂, and 𝑑₃ denote the death rate of the immature prey, mature prey, and the predator; 𝛼 is the density-depended death rate of the predator; 𝑏₁ denotes the per capita per unit time predation rate of the predator; the term 𝑏₂∫^𝑡_−∞𝐹(𝑡−𝑠)𝑥₂(𝑠)d𝑠 is the conversion rate from prey to predator, and the distributed delay may interpret as digest delay. The function 𝐹(𝑠) is called the delayed kernel that is a nonnegative bounded function defined on [0,∞). Following the ideas of Cushing et al. [34], we define 𝐹(𝑡) as the following weak kernel function: 𝐹(𝑡)=𝛽𝑒^−𝛽𝑡,𝛽>0. (2)

Next, we define a new variable: 𝑢(𝑡)=^𝑡_−∞𝛽𝑒^{−𝛽(𝑡−𝑠)}𝑥₂(𝑠)d𝑠. (3) Then by use of linear chain trick technique, system (1) can be transformed into the following equivalent system: 𝑥^′₁(𝑡)=𝑟𝑥₂(𝑡)−𝑑₁𝑥₁(𝑡)−𝑟𝑒^−𝑑1𝜏𝑥₂(),𝑥𝑡−𝜏_′2(𝑡)=𝑟𝑒^−𝑑1𝜏𝑥₂(𝑡−𝜏)−𝑑₂𝑥₂(𝑡)−𝑏₁𝑥₂(𝑡)𝑦(𝑡),𝑢^′(𝑡)=𝛽𝑥₂(𝑡)(𝑡),𝑦−𝛽𝑢_′(𝑡)=𝑦(𝑡)−𝑑₃+𝑏₂𝑢(𝑡)−𝛼𝑦(𝑡). (4)

The organization of this paper is as follows: In Section 2, we will get the conditions for the existence and stability of positive equilibrium of system (4). The occurring condition for Hopf bifurcation is also obtained. In Section 3, by use of normal form theory and central manifold argument, we illustrate the direction and stability of Hopf bifurcation. In Section 4, we give some brief discussion.

2. Stability of Positive Equilibrium and Existence of Hopf Bifurcation

In this section, we will firstly investigate the existence and stability of positive equilibrium of system (4) then study the effect of time delay on the system (4); that is, we will choose 𝜏 as bifurcating parameter to analyze Hopf bifurcation.

Theorem 1. There exists a unique positive equilibrium 𝐸^∗ for system (4), if assumption(H1)𝑟𝑒^−𝑑2𝜏−𝑑₂>0 holds. And 𝐸^∗=(𝑥^∗₁,𝑥^∗₂,𝑢^∗,𝑦^∗), with 𝑥^∗₁=𝑟1−𝑒^−𝑑1𝜏𝛼𝑟𝑒^−𝑑1𝜏−𝛼𝑑₂+𝑏₁𝑏₃𝑏₁𝑏₂𝑑₁,𝑥^∗₂=𝑢^∗=𝛼𝑟𝑒^−𝑑1𝜏−𝑑₂+𝑏₁𝑏₃𝑏₁𝑏₂,𝑦^∗=𝑟𝑒^−𝑑1𝜏−𝑑₂𝑏₁. (5)

Linearizing system (4) at 𝐸^∗, we get 𝑥^′₁(𝑡)=−𝑑₁𝑥₁(𝑡)+𝑟𝑥₂(𝑡)−𝑟𝑒^−𝑑1𝜏𝑥₂(),𝑥𝑡−𝜏_′2(𝑡)=−𝑑₂−𝑏₁𝑦^∗𝑥₂(𝑡)+𝑟𝑒^−𝑑1𝜏𝑥₂(𝑡−𝜏)−𝑏₁𝑥^∗₂𝑦(𝑡),𝑢^′(𝑡)=𝛽𝑥₂(𝑡)(𝑡),𝑦−𝛽𝑢_′(𝑡)=𝑏₂𝑦^∗𝑢(𝑡)−𝛼𝑦^∗𝑦(𝑡), (6) and the characteristic equation for system (6) takes the form 𝜆⁴+ℎ₁𝜆³+ℎ₂𝜆²+ℎ₃𝜆+ℎ₄+ℎ₅𝜆³+ℎ₆𝜆²+ℎ₇𝜆+ℎ₈𝑒^−𝜆𝜏=0, (7) where ℎ₁=𝑑₁+𝑑₂+𝑏₁𝑦^∗+𝛽+𝛼𝑦^∗,ℎ₂=𝑑₁𝑑₂+𝑏₁𝑦^∗+𝛽+𝛼𝑦^∗𝑑₁+𝑑₂+𝑏₁𝑦^∗+𝛼𝛽𝑦^∗,ℎ₃=𝑑₁𝑑₂+𝑏₁𝑦^∗𝛽+𝛼𝑦^∗+𝛼𝛽𝑦^∗𝑑₁+𝑑₂+𝑏₁𝑦^∗,ℎ₄=𝛼𝛽𝑑₁𝑦^∗𝑑₂+𝑏₁𝑦^∗+𝛽𝑏₁𝑏₂𝑥^∗₂𝑦^∗,ℎ₅=−𝑟𝑒^−𝑑1𝜏,ℎ₆=−𝛽+𝛼𝑦^∗+𝑑₁𝑟𝑒^−𝑑1𝜏,ℎ₇=−𝛼𝛽𝑦^∗+𝑑₁𝛽+𝛼𝑑₁𝑦^∗𝑟𝑒^−𝑑1𝜏,ℎ₈=−𝛼𝛽𝑑₁𝑦^∗𝑟𝑒^−𝑑1𝜏. (8)

Note that when 𝜏=0, (7) becomes 𝜆⁴+ℎ₁+ℎ₅𝜆³+ℎ₂+ℎ₆𝜆²+ℎ₃+ℎ₇𝜆+ℎ₄+ℎ₈=0. (9) It is easy to confirm that ℎ₁+ℎ₅=𝑑₁+𝛽+𝛼𝑦^∗ℎ>0,₂+ℎ₆=𝛽+𝛼𝑦^∗𝑑₁+𝛼𝛽𝑦^∗ℎ>0,₃+ℎ₇=𝛼𝛽𝑦^∗𝑑₁ℎ>0,₄+ℎ₈=𝛼𝛽𝑏₁𝑏₂𝑥^∗₂𝑦^∗ℎ>0,₁+ℎ₅ℎ₂+ℎ₆−ℎ₃+ℎ₇ℎ>0,₁+ℎ₅ℎ₂+ℎ₆ℎ₃+ℎ₇−ℎ₁+ℎ₅ℎ₄+ℎ₈−ℎ₃+ℎ₇²>0. (10) Thus, by the Routh-Hurwitz criterion we know that all the roots of (9) have negative real parts, which means that the positive equilibrium 𝐸^∗ is locally asymptotically stable for 𝜏=0.

Next, we will consider the case for 𝜏>0. Suppose that there is a pure imaginary root 𝜆=𝑖𝜔, 𝜔>0. Then we get 𝜔⁴−ℎ₁𝑖𝜔³−ℎ₂𝜔²+ℎ₃𝑖𝜔+ℎ₄+−ℎ₅𝑖𝜔³−ℎ₆𝜔²+ℎ₇𝑖𝜔+ℎ₈(cos𝜔𝜏−𝑖sin𝜔𝜏)=0. (11)

Separating the real and imaginary parts, we have ℎ₆𝜔²−ℎ₈cos𝜔𝜏+ℎ₅𝜔³−ℎ₇𝜔sin𝜔𝜏=𝜔⁴−ℎ₂𝜔²+ℎ₄,ℎ₆𝜔²−ℎ₈sin𝜔𝜏+−ℎ₅𝜔³+ℎ₇𝜔cos𝜔𝜏=ℎ₁𝜔³−ℎ₃𝜔. (12) Incorporating sin²𝜔𝜏+cos²𝜔𝜏=1, we have 𝜔⁸+𝑓₁𝜔⁶+𝑓₂𝜔⁴+𝑓₃𝜔²+𝑓₄=0, (13) where 𝑓₁=ℎ²₁−2ℎ₂−ℎ²₅,𝑓₂=2ℎ₄+ℎ²₂−2ℎ₁ℎ₃−ℎ²₆+2ℎ₅ℎ₇,𝑓₃=ℎ²₃−2ℎ₂ℎ₄+2ℎ₆ℎ₈−ℎ²₇,𝑓₄=ℎ²₄−ℎ²₈. (14)

Denote 𝑧=𝜔². Then (13) becomes 𝑧⁴+𝑓₁𝑧³+𝑓₂𝑧²+𝑓₃𝑧+𝑓₄=0. (15) Let 𝐺(𝑧)=𝑧⁴+𝑓₁𝑧³+𝑓₂𝑧²+𝑓₃𝑧+𝑓₄. (16) Then the following assumption holds true.(H2) Equation (15) has at least one positive real root.

In fact, if all the parameters of system (4) are given, it is easy to calculate the root of (15). Since lim_𝑧→∞𝐺(𝑧)=+∞, we conclude that if 𝑓₄<0, then (15) has at least one positive real root. Without loss of generality, we assume that (15) has four positive root, defined by 𝑧₁, 𝑧₂, 𝑧₃, 𝑧₄, respectively. Then (13) has four positive roots as 𝜔₁=√𝑧₁,𝜔₂=√𝑧₂,𝜔₃=√𝑧₃,𝜔₄=√𝑧₄. (17) From (12), we obtain sin𝜔𝜏=ℎ₅𝜔⁷+ℎ₁ℎ₆−ℎ₇−ℎ₂ℎ₅𝜔⁵+ℎ₂ℎ₇+ℎ₄ℎ₅−ℎ₁ℎ₈−ℎ₃ℎ₆𝜔³+ℎ₃ℎ₈−ℎ₄ℎ₇𝜔×ℎ²₅𝜔⁶+ℎ²₆−2ℎ₅ℎ₇𝜔⁴+ℎ²₇−2ℎ₆ℎ₈𝜔²+ℎ²₈⁻¹,cos𝜔𝜏=ℎ₆−ℎ₁ℎ₅𝜔⁶+ℎ₁ℎ₇+ℎ₃ℎ₅−ℎ₂ℎ₆−ℎ₈𝜔⁴+ℎ₄ℎ₆+ℎ₂ℎ₈−ℎ₃ℎ₇𝜔²+ℎ₄ℎ₈×ℎ²₅𝜔⁶+ℎ²₆−2ℎ₅ℎ₇𝜔⁴+ℎ²₇−2ℎ₆ℎ₈𝜔²+ℎ²₈⁻¹. (18)(19)

Denote 𝑒₁=ℎ²₅,𝑒₂=ℎ²₆−2ℎ₅ℎ₇,𝑒₃=ℎ²₇−2ℎ₆ℎ₈,𝑒₄=ℎ²₈,𝑒₅=ℎ₅,𝑒₆=ℎ₁ℎ₆−ℎ₇−ℎ₂ℎ₅,𝑒₇=ℎ₂ℎ₇+ℎ₄ℎ₅−ℎ₁ℎ₈−ℎ₃ℎ₆,𝑒₈=ℎ₃ℎ₈−ℎ₄ℎ₇,𝑒₉=ℎ₆−ℎ₁ℎ₅,𝑒₁₀=ℎ₁ℎ₇+ℎ₃ℎ₅−ℎ₂ℎ₆−ℎ₈,𝑒₁₁=ℎ₄ℎ₆+ℎ₂ℎ₈−ℎ₃ℎ₇,𝑒₁₂=−ℎ₄ℎ₈. (20) Then cos𝜔𝜏 can be written as 𝑒cos𝜔𝜏=⁹𝜔⁶+𝑒₁₀𝜔⁴+𝑒₁₁𝜔²+𝑒₁₂𝑒₁𝜔⁶+𝑒₂𝜔⁴+𝑒₃𝜔²+𝑒₄, (21) from which we can get 𝜏_𝑘^(𝑗)=1𝜔_𝑘𝑒arccos⁹𝜔⁶_𝑘+𝑒₁₀𝜔⁴_𝑘+𝑒₁₁𝜔²_𝑘+𝑒₁₂𝑒₁𝜔⁶_𝑘+𝑒₂𝜔⁴_𝑘+𝑒₃𝜔²_𝑘+𝑒₄+2𝑗𝜋,𝑘=1,2,3,4;𝑗=0,1,2,…. (22)

Thus, ±𝑖𝜔_𝑘 is a pair of purely imaginary root of (7). Define 𝜏₀=𝜏_𝑘⁽⁰⁾₀=min_{{}𝑘∈1,2,3,4}𝜏_𝑘⁽⁰⁾,𝜔₀=𝜔, (23)

In order to obtain the main result, it is necessary to make the following assumption:(H3)Re(𝑑𝜆/𝑑𝜏)|_𝜏=𝜏0≠0.

Taking the derivative of 𝜆 with respect to 𝜏 in (7), it is easy to obtain 4𝜆³+3ℎ₁𝜆²+2ℎ₂𝜆+ℎ₃𝑑𝜆+𝑑𝜏3ℎ₅𝜆²+2ℎ₆𝜆+ℎ₇𝑒^−𝜆𝜏𝑑𝜆−ℎ𝑑𝜏₅𝜆³+ℎ₆𝜆²+ℎ₇𝜆+ℎ₈𝑒^−𝜆𝜏𝜆+𝜏𝑑𝜆𝑑𝜏=0, (24) and it may be rewritten as 𝑑𝜆=ℎ𝑑𝜏⁵𝜆⁴+ℎ₆𝜆³+ℎ₇𝜆²+ℎ₈𝜆𝑒^−𝜆𝜏×4𝜆³+3ℎ₁𝜆²+2ℎ₂𝜆+ℎ₃+3ℎ₅𝜆²+2ℎ₆𝜆+ℎ₇𝑒^−𝜆𝜏−ℎ₅𝜆³+ℎ₆𝜆²+ℎ₇𝜆+ℎ₈𝜏𝑒^{−𝜆𝜏−1}, (25) or equivalently, we have 𝑑𝜆𝑑𝜏⁻¹=4𝜆³+3ℎ₁𝜆²+2ℎ₂𝜆+ℎ₃+3ℎ₅𝜆²+2ℎ₆𝜆+ℎ₇𝑒^−𝜆𝜏−ℎ₅𝜆³+ℎ₆𝜆²+ℎ₇𝜆+ℎ₈𝜏𝑒^−𝜆𝜏×ℎ₅𝜆⁴+ℎ₆𝜆³+ℎ₇𝜆²+ℎ₈𝜆𝑒^{−𝜆𝜏−1}=4𝜆³+3ℎ₁𝜆²+2ℎ₂𝜆+ℎ₃+3ℎ₅𝜆²+2ℎ₆𝜆+ℎ₇𝑒^−𝜆𝜏×ℎ₅𝜆⁴+ℎ₆𝜆³+ℎ₇𝜆²+ℎ₈𝜆⁻¹−𝜏𝜆. (26)

Taking 𝜆=𝑖𝜔 into the above equation, we get 𝑑𝜆𝑑𝜏⁻¹=−4𝑖𝜔³−3ℎ₁𝜔²+2ℎ₂𝑖𝜔+ℎ₃+−3ℎ₅𝜔²+2ℎ₆𝑖𝜔+ℎ₇(×ℎcos𝜔𝜏−𝑖sin𝜔𝜏)₅𝜔⁴−ℎ₆𝑖𝜔³−ℎ₇𝜔²+ℎ₈𝑖𝜔⁻¹−𝜏=𝑖𝜔−3ℎ₁𝜔²+ℎ₃−3ℎ₅𝜔²−ℎ₇cos𝜔𝜏+2ℎ₆×ℎsin𝜔𝜏₅𝜔⁴−ℎ₇𝜔²+𝑖ℎ₈𝜔−ℎ₆𝜔³⁻¹+𝑖−4𝜔³+2ℎ₂𝜔+2ℎ₆+𝜔cos𝜔𝜏3ℎ₅𝜔²−ℎ₇×ℎsin𝜔𝜏₅𝜔⁴−ℎ₇𝜔²+𝑖ℎ₈𝜔−ℎ₆𝜔³⁻¹−𝜏.𝑖𝜔 (27)

Denote 𝑄=ℎ₅𝜔⁴−ℎ₇𝜔²²+ℎ₈𝜔−ℎ₆𝜔³²>0. (28) Then we have 𝑄Re𝑑𝜆𝑑𝜏⁻¹=−3ℎ₁𝜔²+ℎ₃−3ℎ₅𝜔²−ℎ₇cos𝜔𝜏+2ℎ₆×ℎsin𝜔𝜏₅𝜔⁴−ℎ₇𝜔²+−4𝜔³+2ℎ₂𝜔+2ℎ₆𝜔cos𝜔𝜏+3ℎ₅𝜔²−ℎ₇×ℎsin𝜔𝜏₈𝜔−ℎ₆𝜔³. (29)

Note that SignRe𝑑𝜆∣𝑑𝜏^𝜏=𝜏0=SignRe𝑑𝜆𝑑𝜏⁻¹∣_𝜏=𝜏0. (30) Now, we can employ a result from [35] to analyze (7).

Lemma 2 (see [35]). Consider the exponential polynomial 𝑃𝜆,𝑒^−𝜆𝜏1,𝑒^−𝜆𝜏2,…,𝑒^{−𝜆𝜏𝑚}=𝜆^𝑛+𝑝₁⁽⁰⁾𝜆^𝑛−1+⋯+𝑝⁽⁰⁾_𝑛−1𝜆+𝑝_𝑛⁽⁰⁾+𝑝₁⁽¹⁾𝜆^𝑛−1+⋯+𝑝⁽¹⁾_𝑛−1𝜆+𝑝_𝑛⁽¹⁾𝑒^−𝜆𝜏1+⋯+𝑝₁^(𝑚)𝜆^𝑛−1+⋯+𝑝^(𝑚)_𝑛−1𝜆+𝑝_𝑛^(𝑚)𝑒^{−𝜆𝜏𝑚}, (31) where 𝜏_𝑖≥0(𝑖=1,2,…,𝑚) and 𝑝_𝑗^(𝑖)(𝑖=0,1,…,𝑚;𝑗=1,2,…,𝑛) are constants. As (𝜏₁,𝜏₂,…,𝜏_𝑚) vary, the sum of the order of the zeros of 𝑃(𝜆,𝑒^−𝜆𝜏1,…,𝑒^{−𝜆𝜏𝑚}) on the open right half-plane can change only if a zero appears on or crosses the imaginary axis.

From Lemma 2, it is easy to obtain the following theorem.

Theorem 3. Suppose that (H1), (H2), and (H3) hold. Then the following results hold true.(i)The positive equilibrium 𝐸^∗ of system (4) (𝑥^∗₁,𝑥^∗₂,𝑢^∗,𝑦^∗) is asymptotically stable for 𝜏∈[0,𝜏₀).(ii)The positive equilibrium 𝐸^∗ of system (4) undergoes a Hopf bifurcation when 𝜏=𝜏₀. That is system (4) has a periodic solution bifurcating from the positive equilibrium 𝐸^∗ near 𝜏=𝜏₀.

3. Direction and Stability of the Hopf Bifurcation

In this section, we will derive the explicit formulae for determining the properties of the Hopf bifurcation at critical value of 𝜏₀ by using the normal form and the center manifold theory [35]. Throughout this section, we always assume that system (4) undergoes Hopf bifurcation at the positive equilibrium 𝐸^∗ for 𝜏=𝜏₀, and then ±𝜔₀ is the corresponding purely imaginary roots of the characteristic equation at the positive equilibrium 𝐸^∗.

Let 𝑢₁=𝑥₁−𝑥^∗₁, 𝑢₂=𝑥₂−𝑥^∗₂, 𝑢₃=𝑢−𝑢^∗, 𝑢₄=𝑦−𝑦^∗, 𝑢_𝑖(𝑡)=𝑢_𝑖(𝜏𝑡), and 𝜏=𝜏₀+𝜇, and dropping the bars for simplification of notations, then system (4) is transformed into FDE defined in 𝐶=𝐶([−1,0],𝑅⁴) as 𝑢^′(𝑡)=𝐿_𝜇𝑢_𝑡+𝑔𝜇,𝑥_𝑡, (32) where 𝑢_𝑡=(𝑢₁(𝑡),𝑢₂(𝑡),𝑢₃(𝑡),𝑢₄(𝑡))^𝑇∈𝑅⁴, 𝐿_𝜇∶𝐶→𝑅,𝑔∶𝑅×𝐶→𝑅, and 𝐿_𝜇𝜙=𝜏₀+𝜇−𝑑₁𝑟000−𝑑₂−𝑏₁𝑦^∗0−𝑏₁𝑥^∗₂0𝛽−𝛽000𝑏₂𝑦^∗−𝛼𝑦^∗×𝜙₁(0)𝜙₂(0)𝜙₃(0)𝜙₄(0)+𝜏₀+𝜇0−𝑟𝑒^−𝑑1𝜏000𝑟𝑒^−𝑑1𝜏𝜙0000000000¹()𝜙−1₂()𝜙−1₃()𝜙−1₄(),𝑔=𝜏−1𝜇,𝜙₀0+𝜇−𝑏₁𝜙₂(0)𝜙₄(0)0𝑏₂𝜙₄(0)𝜙₃(0)−𝛼𝜙₄(0)², (33)(34) where 𝜙(𝜃)=(𝜙₁(𝜃),𝜙₂(𝜃),𝜙₃(𝜃),𝜙₄(𝜃))^𝑇∈𝐶. By the Riesz representation theorem, there exists a function 𝜂(𝜃,𝜇) of bounded variation for 𝜃∈[−1,0], such that 𝐿_𝜇𝜙=⁰₋₁d𝜂𝜃,𝜇𝜙(𝜃),for𝜙∈𝐶. (35)

In fact, we can choose 𝜂=𝜏𝜃,𝜇₀×+𝜇−𝑑₁𝑟000−𝑑₂−𝑏₁𝑦^∗0−𝑏₁𝑥^∗₂0𝛽−𝛽000𝑏₂𝑦^∗−𝛼𝑦^∗𝛿(𝜃)+𝜏₀+𝜇0−𝑟𝑒^{−𝑑1(𝜏0+𝜇)}000𝑟𝑒^{−𝑑1(𝜏0+𝜇)}𝛿(),0000000000𝜃+1 (36) where 𝛿 is the Dirac delta function. For 𝜙∈𝐶¹([−1,0],𝑅⁴), define 𝐴𝜇(𝜃)𝜙=d𝜙[),∫d𝜃,𝜃∈−1,0₀₋₁𝜙(𝑠)𝑅𝜇[),𝑔d𝜂𝑠,𝜇,𝜃=0,𝜙=0,𝜃∈−1,0𝜇,𝜙,𝜃=0. (37) Then system (32) is equivalent to 𝑢^′(𝑡)=𝐴𝜇𝑢_𝑡+𝑅𝜇𝑢_𝑡, (38) where 𝑢_𝑡=𝑢(𝑡+𝜃),𝜃∈[−1,0]. For 𝜓∈𝐶¹([0,1],𝑣(𝑅⁴)^∗), define 𝐴^∗𝜓(𝑠)=−d𝜓(𝑠)(],∫d𝑠,𝑠∈0,1₀₋₁d𝜂^𝑇(𝑡,0)𝜓(−𝑡),𝑠=0, (39) and a bilinear inner product 𝜓(𝑠),𝜙(𝜃)=𝜓(0)𝜙(0)−⁰₋₁^𝜃_𝜉=0𝜓(𝜉−𝜃)d𝜂(𝜃)𝜙(𝜉)d𝜉, (40) where 𝜂(𝜃)=𝜂(𝜃,0). Then 𝐴(0) and 𝐴^∗ are adjoint operators. By the discussion in Section 2, we know that ±𝑖𝜔₀𝜏₀ are eigenvalues of 𝐴(0). Thus, they are also eigenvalues of 𝐴^∗. We need to compute the eigenvector of 𝐴(0) and 𝐴^∗ corresponding to 𝑖𝜔₀𝜏₀ and −𝑖𝜔₀𝜏₀, respectively.

Suppose that 𝑞(𝜃)=(1,𝑞₁,𝑞₂,𝑞₃)^𝑇𝑒^{𝑖𝜃𝜔0𝜏0} is the eigenvectors of 𝐴(0) corresponding to 𝑖𝜔₀𝜏₀. Then 𝐴(0)𝑞(𝜃)=𝑖𝜔₀𝜏₀𝑞(𝜃). It follows from the definition of 𝐴(0) and 𝜂(𝜃,𝜇) that 𝜏₀−𝑑₁𝑟000−𝑑₂−𝑏₁𝑦^∗0−𝑏₁𝑥^∗₂0𝛽−𝛽000𝑏₂𝑦^∗−𝛼𝑦^∗𝑞(0)+𝜏₀0−𝑟𝑒^{−𝑑1𝜏0}000𝑟𝑒^{−𝑑1𝜏0}𝑞()0000000000−1=𝑖𝜔₀𝜏₀𝑞(0). (41) Because of 𝑞(−1)=𝑞(0)𝑒^{−𝑖𝜔0𝜏0}, then we get 𝑖𝜔₀+𝑑₁−𝑟+𝑟𝑒^{−𝑑1𝜏0}𝑒^{−𝑖𝜔0𝜏0}000𝑖𝜔₀+𝑑₂+𝑏₁𝑦^∗−𝑟𝑒^{−𝑑1𝜏0}𝑒^{−𝑖𝜔0𝜏0}0𝑏₁𝑥^∗₂0−𝛽𝑖𝜔₀+𝛽000−𝑏₂𝑦^∗𝑖𝜔₀+𝛼𝑦^∗×1𝑞₁𝑞₂𝑞₃=0000, (42) from which we obtain 𝑞₁=𝑖𝜔₀−𝑑₁𝑟1−𝑒^{−𝑑1𝜏0}𝑒^{−𝑖𝜔0𝜏0},𝑞₂=𝛽𝑑₁−𝑖𝜔₀𝑟𝑖𝜔₀+𝛽1−𝑒^{−𝑑1𝜏0}𝑒^{−𝑖𝜔0𝜏0},𝑞₃=𝑏₂𝑦^∗𝛽𝑖𝜔₀+𝑑₁𝑟𝑖𝜔₀+𝛽𝑖𝜔₀+𝛼𝑦^∗1−𝑒^{−𝑑1𝜏0}𝑒^{−𝑖𝜔0𝜏0}. (43)

Similarly, let 𝑞^∗(𝜃)=𝐷(1,𝑞^∗₁,𝑞^∗₂,𝑞^∗₃)𝑒^{𝑖𝜃𝜔0𝜏0} be the eigenvectors of 𝐴^∗ corresponding to −𝑖𝜔₀𝜏₀, and according to the definition of 𝐴^∗ we have 𝜏₀−𝑑₁000𝑟−𝑑₂−𝑏₁𝑦^∗𝛽000−𝛽𝑏₂𝑦^∗0−𝑏₁𝑥^∗₂0−𝛼𝑦^∗𝑞^∗(0)+𝜏₀0000−𝑟𝑒^{−𝑑1𝜏0}𝑟𝑒^{−𝑑1𝜏0}𝑞0000000000^∗(−1)=−𝑖𝜔₀𝜏₀𝑞^∗(0). (44) Note that 𝑞^∗(−1)=𝑞^∗(0)𝑒^{𝑖𝜔0𝜏0}. Then we get −𝑖𝜔₀+𝑑₁000𝑟𝑒^{−𝑑1𝜏0}𝑒^{𝑖𝜔0𝜏0}−𝑟𝑑₂−𝑖𝜔₀+𝑏₁𝑦^∗−𝑟𝑒^{−𝑑1𝜏0}𝑒^{𝑖𝜔0𝜏0}−𝛽000𝛽−𝑖𝜔₀−𝑏₂𝑦^∗0𝑏₁𝑥^∗₂0𝛼𝑦^∗−𝑖𝜔₀×1𝑞^∗₁𝑞^∗₂𝑞^∗₃=0000, (45) from which we can obtain 𝑞^∗₁=𝑖𝜔₀−𝛼𝑦^∗𝑖𝜔₀+𝛽1−𝑒^{−𝑑1𝜏0}𝑒^{𝑖𝜔0𝜏0}×𝑖𝜔₀+𝛽𝑖𝜔₀−𝛼𝑦^∗×−𝑖𝜔₀+𝑑₂+𝑏₁𝑦^∗−𝑟𝑒^{−𝑑1𝜏0}𝑒^{−𝑖𝜔0𝜏0}−𝛽𝑏₁𝑥^∗₂𝑏₂𝑦^∗−1,𝑞^∗₂=𝑟𝑏₁𝑥^∗₂𝑏₂𝑦^∗𝑖𝜔₀+𝛽1−𝑒^{−𝑑1𝜏0}𝑒^{𝑖𝜔0𝜏0}×𝜔²₀+𝛽²𝑖𝜔₀−𝛼𝑦^∗×−𝑖𝜔₀+𝑑₂+𝑏₁𝑦^∗−𝑟𝑒^{−𝑑1𝜏0}𝑒^{−𝑖𝜔0𝜏0}−𝛽−𝑖𝜔₀𝛽𝑏₁𝑥^∗₂𝑏₂𝑦^∗−1,𝑞^∗₃=𝑟𝑏₁𝑥^∗₂𝑖𝜔₀+𝛽1−𝑒^{−𝑑1𝜏0}𝑒^{𝑖𝜔0𝜏0}×𝑖𝜔₀+𝛽𝑖𝜔₀−𝛼𝑦^∗×−𝑖𝜔₀+𝑑₂+𝑏₁𝑦^∗−𝑟𝑒^{−𝑑1𝜏0}𝑒^{−𝑖𝜔0𝜏0}−𝛽𝑏₁𝑥^∗₂𝑏₂𝑦^∗−1. (46) By (40), we get 𝑞^∗(𝑠)(𝜃=,𝑞)𝐷1,𝑞^∗₁,𝑞^∗₂,𝑞^∗₃1,𝑞₁,𝑞₂,𝑞₃^𝑇−⁰₋₁^𝜃_𝜉=0𝐷1,𝑞^∗₁,𝑞^∗₂,𝑞^∗₃×𝑒^{−𝑖𝜔0𝜏0(𝜉−𝜃)}d𝜂(𝜃)1,𝑞₁,𝑞₂,𝑞₃^𝑇𝑒^{𝑖𝜔0𝜏0𝜉}=d𝜉𝐷1+𝑞₁𝑞^∗₁+𝑞₂𝑞^∗₂+𝑞₃𝑞^∗₃−⁰₋₁1,𝑞^∗₁,𝑞^∗₂,𝑞^∗₃𝜃𝑒^{𝑖𝜔0𝜏0𝜃}𝑑𝜂(𝜃)1,𝑞₁,𝑞₂,𝑞₃^𝑇=𝐷1+𝑞₁𝑞^∗₁+𝑞₂𝑞^∗₂+𝑞₃𝑞^∗₃+𝜏₀𝑟𝑞₁𝑞^∗₁−1𝑒^{−𝑑2𝜏0}𝑒^{−𝑖𝜔0𝜏0}. (47) Then we can choose 𝐷 such that ⟨𝑞^∗(𝑠),𝑞(𝜃)⟩=1, ⟨𝑞^∗(𝑠),𝑞(𝜃)⟩=0.

Nest, we will use the ideas in [35] to compute the coordinates describing center manifold 𝐶₀ at 𝜇=0. Define 𝑧(𝑡)=𝑞^∗,𝑢_𝑡,𝑊(𝑡,𝜃)=𝑢_𝑡(𝜃)−2Re𝑧(𝑡),𝑞(𝜃). (48)

On the center manifold 𝐶₀, we have 𝑊(𝑡,𝜃)=𝑊(𝑧(𝑡),𝑧(𝑡),𝜃), and 𝑊(𝑧(𝑡),𝑧(𝑡),𝜃)=𝑊₂₀(𝜃)𝑧²2+𝑊₁₁(𝜃)𝑧𝑧+𝑊₀₂𝑧²2+…, (49) where 𝑧 and 𝑧 are local coordinates for center manifold 𝐶₀ in the direction of 𝑞^∗(𝑠) and 𝑞(𝜃). Note that 𝑊 is real if 𝑢_𝑡 is real, and we only consider real solutions. For the solution 𝑢_𝑡∈𝐶₀ of (38), since 𝜇=0 and (38), we have 𝑧^′(𝑡)=𝑖𝜔₀𝜏₀𝑧+𝑞^∗(0)𝑔(0,𝑊(𝑧,𝑧,0))+2Re𝑧𝑞(0)=𝑖𝜔₀𝜏₀𝑧+𝑞^∗(0)𝑔₀(𝑧,𝑧). (50) Then, the above equation can be denoted as 𝑧^′(𝑡)=𝑖𝜔₀𝜏₀𝑧(𝑡)+𝑓(𝑧,𝑧), (51) where 𝑓(𝑧,𝑧)=𝑞^∗(0)𝑔₀(𝑧,𝑧)=𝑓₂₀𝑧²2+𝑓₁₁𝑧𝑧+𝑓₀₂𝑧²2+𝑓₂₁𝑧²𝑧2+⋯. (52) From (48) and (49), we have 𝑢_𝑡=𝑢_1𝑡(𝜃),𝑢_2𝑡(𝜃),𝑢_3𝑡(𝜃),𝑢_4𝑡(𝜃)()𝑧(𝑡)(𝜃()(𝜃)+=𝑊𝑡,𝜃+2Re,𝑞)=𝑊𝑡,𝜃+𝑧𝑞𝑧𝑞(𝜃), (53) and then we can obtain 𝑢_2𝑡(0)=𝑊⁽²⁾₂₀(0)𝑧²2+𝑊⁽²⁾₁₁(0)𝑧𝑧+𝑊⁽²⁾₀₂(0)𝑧²2+𝑞₁𝑧+𝑞₁𝑧+𝑂𝑧,𝑧³,𝑢_3𝑡(0)=𝑊⁽³⁾₂₀(0)𝑧²2+𝑊⁽³⁾₁₁(0)𝑧𝑧+𝑊⁽³⁾₀₂(0)𝑧²2+𝑞₂𝑧+𝑞₂𝑧+𝑂|𝑧,𝑧|³,𝑢_4𝑡(0)=𝑊⁽⁴⁾₂₀(0)𝑧²2+𝑊⁽⁴⁾₁₁(0)𝑧𝑧+𝑊⁽⁴⁾₀₂(0)𝑧²2+𝑞₃𝑧+𝑞₃𝑧+𝑂|(𝑧,𝑧)|³. (54)

From the definition of 𝑔(𝜇,𝑢_𝑡), we have 𝑓(𝑧,𝑧)=𝑞^∗(0)𝑔₀(𝑧,𝑧)=𝜏₀𝐷1,𝑞^∗₁,𝑞^∗₂,𝑞^∗₃×0−2𝑏₁𝑢_2𝑡(0)𝑢_4𝑡(0)0𝑏₂𝑢_3𝑡(0)𝑢_4𝑡(0)−𝛼𝑢²_4𝑡(0)=𝜏₀𝐷𝑧²−2𝑏₁𝑞^∗₁𝑞₁𝑞₃+𝑏₂𝑞^∗₃𝑞₂𝑞₃−𝛼𝑞^∗₃𝑞²₃+2𝑧𝑧−2𝑏₁𝑞^∗₁Re𝑞₁𝑞₃+𝑏₂𝑞^∗₃Re𝑞₂𝑞₃−𝛼𝑞^∗₃Re𝑞₃𝑞₃+𝑧²−𝑏₁𝑞^∗₁𝑞₁𝑞₃+𝑏₂𝑞^∗₃𝑞₂𝑞₃−𝛼𝑞^∗₃𝑞₃𝑞₃+12𝑧²𝑧−𝑏₁𝑞^∗₁𝑊⁽²⁾₂₀(0)𝑞₃+2𝑊⁽²⁾₁₁(0)𝑞₃+2𝑊⁽⁴⁾₁₁(0)𝑞₁+𝑊⁽⁴⁾₂₀(0)𝑞₁+𝑏₂𝑞^∗₃𝑊⁽³⁾₂₀(0)𝑞₃+2𝑊⁽³⁾₁₁(0)𝑞₃+2𝑊⁽⁴⁾₁₁(0)𝑞₂+𝑊⁽⁴⁾₂₀(0)𝑞₂−𝛼𝑞^∗₃2𝑊⁽⁴⁾₂₀(0)𝑞₃+4𝑊⁽⁴⁾₁₁(0)𝑞₃+⋯}. (55)(56)

Comparing the coefficients with those of (52), we obtain 𝑓₂₀=2𝜏₀𝐷−𝑏₁𝑞^∗₁𝑞₁𝑞₃+𝑏₂𝑞^∗₃𝑞₂𝑞₃−𝛼𝑞^∗₃𝑞²₃,𝑓₁₁=2𝜏₀𝐷−𝑏₁𝑞^∗₁Re𝑞₁𝑞₃+𝑏₂𝑞^∗₃Re𝑞₂𝑞₃−𝛼𝑞^∗₃Re𝑞₃𝑞₃,𝑓₀₂=2𝜏₀𝐷−𝑏₁𝑞^∗₁𝑞₁𝑞₃+𝑏₂𝑞^∗₃𝑞₂𝑞₃−𝛼𝑞^∗₃𝑞²₃,𝑓₂₁=𝜏₀𝐷−𝑏₁𝑞^∗₁𝑊⁽²⁾₂₀(0)𝑞₃+2𝑊⁽²⁾₁₁(0)𝑞₃+2𝑊⁽⁴⁾₁₁(0)𝑞₁+𝑊⁽⁴⁾₂₀(0)𝑞₁+𝑏₂𝑞^∗₃𝑊⁽³⁾₂₀(0)𝑞₃+2𝑊⁽³⁾₁₁(0)𝑞₃+2𝑊⁽⁴⁾₁₁(0)𝑞₂+𝑊⁽⁴⁾₂₀(0)𝑞₂−𝛼𝑞^∗₃2𝑊⁽⁴⁾₂₀(0)𝑞₃+4𝑊⁽⁴⁾₁₁(0)𝑞₃. (57)

In order to determine 𝑓₂₁ we need to compute 𝑊₂₀(𝜃) and 𝑊₁₁(𝜃). From (38) and (48), we have 𝑊̇=𝑢̇_𝑡−𝑧𝑞̇−̇=𝐴(0)𝑧𝑞𝑊−2Re𝑞^∗(0)𝑔₀𝑞(𝜃[),𝐴(0)),𝜃∈−1,0𝑊−2Re𝑞^∗(0)𝑔₀𝑞(𝜃)+𝑔₀(0)(,𝜃=0,=𝐴𝑊−𝐻𝑧,𝑧,𝜃), (58) where 𝐻(𝑧,𝑧,𝜃)=𝐻₂₀(𝜃)𝑧²2+𝐻₁₁(𝜃)𝑧𝑧+𝐻₀₂𝑧²2+⋯. (59)

Note that on the center manifold 𝐶₀ near the origin 𝑊̇=𝑊_𝑧𝑧̇+𝑊_𝑧𝑧̇‌, (60) and thus we obtain 𝐴(0)−2𝑖𝜔₀𝜏₀𝑊₂₀(𝜃)=−𝐻₂₀(𝜃),𝐴(0)𝑊₁₁(𝜃)=−𝐻₁₁(𝜃). (61)

By (58), we know that, for 𝜃∈[−1,0), 𝐻(𝑧,𝑧,𝜃)=−𝑞^∗(0)𝑔₀𝑞(𝜃)−𝑞^∗₀𝑔₀𝑞(𝜃)=−𝑓𝑞(𝜃)−𝑓𝑞(𝜃). (62)

Comparing the coefficients with those in (59), we get 𝐻₂₀(𝜃)=−𝑓₂₀𝑞(𝜃)−𝑓₀₂𝑞(𝜃),𝐻₁₁(𝜃)=−𝑓₁₁𝑞(𝜃)−𝑓₁₁𝑞(𝜃). (63)

From (61), (63), and the definition of 𝐴, we have 𝑊̇₂₀(𝜃)=2𝑖𝜔₀𝜏₀𝑊₂₀(𝜃)+𝑓₂₀𝑞(𝜃)+𝑓₀₂𝑞(𝜃). (64)

Noting that 𝑞(𝜃)=𝑞(0)𝑒^{𝑖𝜔0𝜏0𝜃}, we get 𝑊₂₀(𝜃)=𝑖𝑓₂₀𝜔₀𝜏₀𝑞(0)𝑒^{𝑖𝜔0𝜏0𝜃}+𝑖𝑓₀₂3𝜔₀𝜏₀𝑞(0)𝑒^{−𝑖𝜔0𝜏0𝜃}+𝐸₁𝑒^{2𝑖𝜔0𝜏0𝜃}, (65) where 𝐸₁=(𝐸₁⁽¹⁾,𝐸₁⁽²⁾,𝐸₁⁽³⁾,𝐸₁⁽⁴⁾)∈𝑅⁴ is a constant vector.

Similarly, from (61) and (62), we can get 𝑊₁₁(𝜃)=𝑖𝑓₁₁𝜔₀𝜏₀𝑞(0)𝑒^{−𝑖𝜔0𝜏0𝜃}+𝐸₂, (66) where 𝐸₂=(𝐸₂⁽¹⁾,𝐸₂⁽²⁾,𝐸₂⁽³⁾,𝐸₂⁽⁴⁾)∈𝑅⁴ is a constant vector.

Next, we will find out 𝐸₁ and 𝐸₂. In fact, from the definition of 𝐴 and (61), we can obtain d𝜂(⁰₋₁𝜃)𝑊₂₀(𝜃)=2𝑖𝜔₀𝜏₀𝑊₂₀(0)−𝐻₂₀(0),₀₋₁d𝜂(𝜃)𝑊₁₁(𝜃)=−𝐻₁₁(0), (67)(68) where 𝜂(𝜃)=𝜂(0,𝜃). From (59) and (61), when 𝜃=0 we have 𝐻(𝑧,𝑧,0)=−2Re𝑞^∗(0)𝑔₀𝑞(0)+𝑔₀=−𝑞^∗(0)𝑔₀𝑞(0)−𝑞^∗(0)𝑔₀𝑞(0)+𝑔₀=−𝑔(𝑧,𝑧)𝑞(0)−𝑔(𝑧,𝑧)𝑞(0)+𝑔₀. (69) That is, 𝐻₂₀(0)𝑧²2+𝐻₁₁(0)𝑧𝑧+𝐻₀₂(0)𝑧²2=−𝑞(0)𝑓₂₀𝑧²2+𝑓₁₁𝑧𝑧+𝑓₀₂𝑧²2−+⋯𝑞(0)𝑓₂₀𝑧²2+𝑓₁₁𝑧𝑧+𝑓₀₂𝑧²2+⋯+𝑔₀. (70)

By (34), we have 𝑔₀=𝜏₀0−𝑏₁𝑢_2𝑡(0)𝑢_4𝑡(0)0𝑏₂𝑢_3𝑡(0)𝑢_4𝑡(0)−𝛼𝑢²_4𝑡(0). (71)

By (48) and (49), we have 𝑢_𝑡(𝜃)()𝑧(𝑡)𝑞(𝜃()(𝑡)𝑞(𝜃)+=𝑊𝑡,𝜃+2Re)=𝑊𝑡,𝜃+𝑧𝑧(𝑡)𝑞(𝜃)=𝑊₂₀(𝜃)𝑧²2+𝑊₁₁(𝜃)𝑧(𝑡)𝑞(𝜃)+𝑧+𝑧𝑧(𝑡)𝑞(𝜃)+⋯. (72)

Thus, we obtain 𝐻₂₀(0)=−𝑓₂₀𝑞(0)−𝑓₀₂𝑞(0)+2𝜏₀0−𝑏₁𝑞₁𝑞₃0𝑏₂𝑞₂𝑞₃−𝛼𝑞²₃𝐻₁₁,(0)=−𝑓₁₁𝑞(0)−𝑓₁₁𝑞(0)+2𝜏₀×0−𝑏₁Re𝑞₁𝑞₃0𝑏₂Re𝑞₂𝑞₃−𝛼𝑞^∗₃𝑞₃². (73)(74)

Note that 𝑖𝜔₀𝜏₀𝐼−⁰₋₁𝑒^{𝑖𝜔0𝜏0𝜃}𝑑𝜂(𝜃)𝑞(0)=0,−𝑖𝜔₀𝜏₀𝐼−⁰₋₁𝑒^{−𝑖𝜔0𝜏0𝜃}𝑑𝜂(𝜃)𝑞(0)=0. (75) Then, by substituting (65) and (73) into (67), we will get 2𝑖𝜔₀𝜏₀𝐼−⁰₋₁𝑒^{2𝑖𝜔0𝜏0𝜃}𝑑𝜂(𝜃)𝐸₁=2𝜏₀0−𝑏₁𝑞₁𝑞₃0𝑏₂𝑞₂𝑞₃−𝛼𝑞²₃, (76) or equivalently, 2𝑖𝜔₀+𝑑₁−𝑟+𝑟𝑒^{−𝑑1𝜏0}𝑒^{−2𝑖𝜔0𝜏0}0002𝑖𝜔₀+𝑑₂+𝑏₁𝑦^∗−𝑟𝑒^{−𝑑1𝜏0}𝑒^{−2𝑖𝜔0𝜏0}0𝑏₁𝑥^∗₂0−𝛽2𝑖𝜔₀+𝛽000−𝑏₂𝑦^∗2𝑖𝜔₀−𝛼𝑦^∗𝐸₁0=2−𝑏₁𝑞₁𝑞₃0𝑏₂𝑞₂𝑞₃−𝛼𝑞²₃, (77) from which we can get 𝐸₁=22𝑖𝜔₀+𝑑₁−𝑟+𝑟𝑒^{−𝑑1𝜏0}𝑒^{−2𝑖𝜔0𝜏0}0002𝑖𝜔₀+𝑑₂+𝑏₁𝑦^∗−𝑟𝑒^{−𝑑1𝜏0}𝑒^{−2𝑖𝜔0𝜏0}0𝑏₁𝑥^∗₂0−𝛽2𝑖𝜔₀+𝛽000−𝑏₂𝑦^∗2𝑖𝜔₀−𝛼𝑦^∗−1×0−𝑏₁𝑞₁𝑞₃0𝑏₂𝑞₂𝑞₃−𝛼𝑞²₃. (78)

Similarly, by substituting (66) and (74) into (68), we will get 𝑑₁−𝑟+𝑟𝑒^{−𝑑1𝜏0}000𝑑₂+𝑏₁𝑦^∗−𝑟𝑒^{−𝑑1𝜏0}0𝑏₁𝑥^∗₂0−𝛽𝛽000−𝑏₂𝑦^∗−𝛼𝑦^∗𝐸₂0=2−𝑏₁Re𝑞₁𝑞₃0𝑏₂Re𝑞₂𝑞₃−𝛼Re𝑞₃𝑞₃,𝐸₂𝑑=2¹−𝑟+𝑟𝑒^{−𝑑1𝜏0}000𝑑₂+𝑏₁𝑦^∗−𝑟𝑒^{−𝑑1𝜏0}0𝑏₁𝑥^∗₂0−𝛽𝛽000−𝑏₂𝑦^∗−𝛼𝑦^∗−1×0−𝑏₁Re𝑞₁𝑞₃0𝑏₂Re𝑞₂𝑞₃−𝛼Re𝑞₃𝑞₃ (79)

Thus, we can determine 𝑊₂₀(𝜃) and 𝑊₁₁(𝜃) from (65) and (66). Furthermore, 𝑓₂₁ can be expressed by the parameters and delay. Thus, we can compute the following values: 𝑐₁(0)=𝑖2𝜔₀𝜏₀𝑓₂₀𝑓₁₁−2𝑓₁₁²−𝑓₀₂²3+𝑓₂₁2,𝜇₂𝑐=−Re₁(0)Re𝜆^′𝜏₀,𝛽₂=2Re𝑐₁(0),𝑇₂𝑐=−Im₁(0)+𝜇₂Im𝜆^′𝜏₀𝜔𝜏₀, (80) which determine the qualities of bifurcating periodic solution in the center manifold at critical value 𝜏₀; that is, 𝜇₂ determines the direction of the Hopf bifurcation: if 𝜇₂>0(𝜇₂<0), then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solution exists for 𝜏>𝜏₀(𝜏<𝜏₀);𝛽₂ determines the stability of the bifurcating periodic solution: the bifurcating periodic solution is stable (unstable) if 𝛽₂<0(𝛽₂>0); and 𝑇₂ determines period of the bifurcating periodic solution: the period increases (decreases) if 𝑇₂>0(<0).

4. Discussion

In this paper, we propose a three-dimensional stage-structured predatory-prey model with discrete and distributed delay. Then, by introducing a new variable, the original system is transformed into an equivalent four-dimensional system (4). In Section 2, we analyze the existence and local stability of the positive equilibrium of the system (4). In fact, by use of the Routh-Hurwitz criterion, we know that the positive equilibrium 𝐸^∗ is locally asymptotically stable for 𝜏=0. Then, by some computation, we get a threshold 𝜏=𝜏₀ in Theorem 3. The results indicate that if 𝜏<𝜏₀, the positive equilibrium 𝐸^∗ of system (4) is asymptotically stable; if 𝜏=𝜏₀, the positive equilibrium 𝐸^∗ of system (4) will undergo a Hopf bifurcation. In Section 3, by use of normal form theory and central manifold argument, we establish the formulae for the direction and the stability of the Hopf bifurcation. Our theoretical results show that stage structure and delay play an important role in the dynamics of the system, and delay may lead to complicated dynamic behaviors, such as Hopf bifurcation.

Acknowledgments

The first author is supported by Postdoctoral Science Foundation of China (no. 2011M501428) and Young Science Funds of Shanxi (no. 2013021002-2). The third author is supported by National Natural Science Foundation of China (no. 11171199). The authors would like to thank the referees of this paper for their helpful comments.