Hopf Bifurcation of a Predator-Prey System with Delays and Stage Structure for the Prey

Abstract

This paper is concerned with a Holling type III predator-prey system with stage structure for the prey population and two time delays. The main result is given in terms of local stability and bifurcation. By choosing the time delay as a bifurcation parameter, sufficient conditions for the local stability of the positive equilibrium and the existence of periodic solutions via Hopf bifurcation with respect to both delays are obtained. In particular, explicit formulas that can determine the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions are established by using the normal form method and center manifold theorem. Finally, numerical simulations supporting the theoretical analysis are also included.

1. Introduction

Predator-prey dynamics continues to draw interest from both applied mathematicians and ecologists due to its universal existence and importance. Many kinds of predator-prey models have been studied extensively [1–6]. It is well known that there are many species whose individual members have a life history that takes them through immature stage and mature stage. To analyze the effect of a stage structure for the predator or the prey on the dynamics of a predator-prey system, many scholars have investigated predator-prey systems with stage structure in the last two decades [7–15]. In [7], Wang considered the following predator-prey system with stage structure for the predator and obtained the sufficient conditions for the global stability of a coexistence equilibrium of the system: 𝑑𝑥(𝑡(𝑡−𝑎𝑑𝑡=𝑥)(𝑟−𝑎𝑥))¹𝑦₂(𝑡)𝑥(𝑡)(𝑡),1+𝑚𝑥𝑑𝑦₁=𝑎𝑑𝑡²𝑥(𝑡)𝑦₂(𝑡)1+𝑚𝑥(𝑡)−𝑟₁𝑦₁(𝑡)−𝐷𝑦₁(𝑡),𝑑𝑦₂𝑑𝑡=𝐷𝑦₁(𝑡)−𝑟₂𝑦₂(𝑡), (1.1) where 𝑥(𝑡) represents the density of the prey at time 𝑡. 𝑦₁(𝑡) and 𝑦₂(𝑡) represent the densities of the immature predator and the mature predator at time 𝑡, respectively. For the meanings of all the parameters in system (1.1), one can refer to [7]. Considering the gestation time of the mature predator, Xu [8] incorporated the time delay due to the gestation of the mature predator into system (1.1) and considered the effect of the time delay on the dynamics of system (1.1).

There has also been a significant body of work on the predator-prey system with stage structure for the prey. In [12], Xu considered a delayed predator-prey system with a stage structure for the prey: 𝑑𝑥₁𝑑𝑡=𝑎𝑥₂(𝑡)−𝑟₁𝑥₁(𝑡)−𝑏𝑥₁(𝑡),𝑑𝑥₂𝑑𝑡=𝑏𝑥₁(𝑡)−𝑟₂𝑥₂(𝑡)−𝑏₁𝑥²₂(𝑡)−𝑎₁𝑥₂(𝑡)𝑦(𝑡)1+𝑚𝑥₂(𝑡),𝑑𝑦=𝑎𝑑𝑡²𝑥₂(𝑡−𝜏)𝑦(𝑡−𝜏)1+𝑚𝑥₂()(𝑡),𝑡−𝜏−𝑟𝑦 (1.2) where 𝑥₁(𝑡) and 𝑥₂(𝑡) denote the population densities of the immature prey and the mature prey at time 𝑡, respectively. 𝑦(𝑡) denotes the population density of the predator at time 𝑡. All the parameters in system (1.2) are assumed positive. 𝑎 is the birth rate of the immature prey. 𝑏 is the transformation rate from immature individual to mature individuals. 𝑏₁ is the intraspecific competition coefficient of the mature prey. 𝑟₁ and 𝑟₂ are the death rates of the immature and the mature prey, respectively. 𝑟 is the death rate of the predator. 𝑎₁ and 𝑎₂ are the interspecific interaction coefficients between the mature prey and the predator, respectively. 𝑎₁𝑥₂/(1+𝑚𝑥₂) is the response function of the predator. And 𝜏 is a constant delay due to the gestation of the predator. In [12], Xu investigated the persistence of system (1.2) by means of the persistence theory on infinite dimensional systems, and sufficient conditions are obtained for the global stability of nonnegative equilibrium of the model by constructing appropriate Lyapunov function. But studies on the predator-prey system not only involve the persistence and stability, but also involve many other behaviors such as periodic phenomenon, attractivity, and bifurcation [16–19]. In particular, the properties of periodic solutions are of great interest [20–24]. Therefore, F. Li and H. W. Li [14] considered the property of periodic solutions of the following system: 𝑑𝑥₁𝑑𝑡=𝑎𝑥₂(𝑡)−𝑟₁𝑥₁(𝑡)−𝑏𝑥₁(𝑡),𝑑𝑥₂𝑑𝑡=𝑏𝑥₁(𝑡)−𝑟₂𝑥₂(𝑡)−𝑏₁𝑥²₂(𝑡)−𝑎₁𝑥²₂(𝑡)𝑦(𝑡)1+𝑚𝑥²₂(𝑡),𝑑𝑦=𝑎𝑑𝑡²𝑥²₂(𝑡−𝜏)𝑦(𝑡−𝜏)1+𝑚𝑥²₂()(𝑡).𝑡−𝜏−𝑟𝑦 (1.3) Motivated by the work of Xu [12] and F. Li and H. W. Li [14] and considering the intraspecific competition of the immature prey population, we consider the following system: 𝑑𝑥₁𝑑𝑡=𝑎𝑥₂(𝑡)−𝑟₁𝑥₁(𝑡)−𝑏𝑥₁(𝑡)−𝑐𝑥²₁(𝑡),𝑑𝑥₂𝑑𝑡=𝑏𝑥₁(𝑡)−𝑟₂𝑥₂(𝑡)−𝑏₁𝑥₂(𝑡)𝑥₂𝑡−𝜏₁−𝑎₁𝑥²₂(𝑡)𝑦(𝑡)1+𝑚𝑥²₂(𝑡),𝑑𝑦=𝑎𝑑𝑡²𝑥²₂𝑡−𝜏₂𝑦𝑡−𝜏₂1+𝑚𝑥²₂𝑡−𝜏₂−𝑟𝑦(𝑡), (1.4) where 𝑥₁(𝑡) and 𝑥₂(𝑡) denote the population densities of the immature prey and the mature prey at time 𝑡, respectively. 𝑦(𝑡) denotes the population density of the predator at time 𝑡. The parameters 𝑎, 𝑎₁, 𝑎₂, 𝑏, 𝑏₁, 𝑟, 𝑟₁, 𝑟₂, and 𝑚 are defined as in system (1.3). 𝑐 is the intraspecific competition of the immature prey, 𝜏₁is the feedback delay of the mature prey, and 𝜏₂ is the time delay due to the gestation of the predator.

The organization of this paper is as follows. In Section 2, by analyzing the corresponding characteristic equations, the local stability of the positive equilibrium of system (1.4) is discussed, and the existence of Hopf bifurcation at the positive equilibrium is established. In Section 3, we determine the direction of Hopf bifurcation and the stability of bifurcating periodic solutions by using the normal form theory and center manifold theorem in [20]. And numerical simulations are carried out in Section 4 to illustrate the main theoretical results. Finally, main conclusions are included.

2. Local Stability and Hopf Bifurcation

From the viewpoint of biology, we are only interested in the positive equilibrium of system (1.4). It is not difficult to verify that system (1.4) has a positive equilibrium 𝐸⁰(𝑥⁰₁,𝑥⁰₂,𝑦⁰), where 𝑥⁰₁=−𝑏+𝑟₁+𝑏+𝑟₁²+4𝑎𝑐𝑥⁰₂,𝑥2𝑐₀₂=𝑟𝑎₂,𝑦−𝑚𝑟₀=𝑏𝑥⁰₁−𝑟₂𝑥⁰₂−𝑏₁𝑥⁰₂²1+𝑚𝑥⁰₂²𝑎₁𝑥⁰₂², (2.1) if the following conditions hold: 𝐻₁∶𝑎₂>𝑚𝑟, 𝐻₂∶𝑏𝑥⁰₁>(𝑟₂+𝑏₁𝑥⁰₂)𝑥⁰₂.

Let 𝑥₁(𝑡)=𝑧₁(𝑡)+𝑥⁰₁, 𝑥₂(𝑡)=𝑧₂(𝑡)+𝑥⁰₂, 𝑦(𝑡)=𝑧₃(𝑡)+𝑦⁰, and we still denote 𝑧₁(𝑡), 𝑧₂(𝑡), and 𝑧₃(𝑡) by 𝑥₁(𝑡), 𝑥₂(𝑡), and 𝑦(𝑡). Then system (1.4) can be transformed to the following form: 𝑑𝑥₁𝑑𝑡=𝛼₁₁𝑥₁(𝑡)+𝛼₁₂𝑥₂(𝑡)+_{𝑖+𝑗≥2}𝑓₁^(𝑖𝑗)𝑥^𝑖₁𝑥^𝑗₂,𝑑𝑥₂𝑑𝑡=𝛼₂₁𝑥₁(𝑡)+𝛼₂₂𝑥₂(𝑡)+𝛼₂₃𝑦(𝑡)+𝛽₂₂𝑥₂𝑡−𝜏₁+_{𝑖+𝑗+𝑘+𝑙≥2}𝑓₂^{(𝑖𝑗𝑘𝑙)}𝑥^𝑖₁𝑥^𝑗₂𝑦^𝑘𝑥^𝑙₂𝑡−𝜏₁,𝑑𝑦𝑑𝑡=𝛼₃₃𝑦(𝑡)+𝛾₃₂𝑥₂𝑡−𝜏₂+𝛾₃₃𝑦𝑡−𝜏₂+_{𝑖+𝑗+𝑘≥2}𝑓₃^{(𝑖𝑗𝑘)}𝑦^𝑖𝑥^𝑗₂𝑡−𝜏₂𝑦^𝑘𝑡−𝜏₂, (2.2) where 𝛼₁₁=−𝑟₁−𝑏−2𝑐𝑥⁰₁,𝛼₁₂=𝑎,𝛼₂₁𝛼=𝑏,₂₂=−𝑟₂−𝑏₁𝑥⁰₂−2𝑎₁𝑥⁰₂𝑦⁰1+𝑚𝑥⁰₂²²,𝛼₂₃𝑎=−¹𝑥⁰₂²1+𝑚𝑥⁰₂²,𝛼₃₃=−𝑟,𝛽₂₂=−𝑏₁𝑥⁰₂,𝛾₃₂=2𝑎₂𝑥⁰₂𝑦⁰1+𝑚𝑥⁰₂²²,𝛾₃₃𝑓=𝑟,_1(𝑖𝑗)=1𝜕𝑖!𝑗!^𝑖+𝑗𝑓₁𝜕𝑥^𝑖₁(𝑡)𝜕𝑥^𝑗₂(𝑡)∣𝑥⁰₁,𝑥⁰₂,𝑦⁰,𝑓₂^{(𝑖𝑗𝑘𝑙)}=1𝜕𝑖!𝑗!𝑘!𝑙!^{𝑖+𝑗+𝑘+𝑙}𝑓₂𝜕𝑥^𝑖₁(𝑡)𝜕𝑥^𝑗₂(𝑡)𝜕𝑦^𝑘(𝑡)𝜕𝑥^𝑙₂𝑡−𝜏₁∣𝑥⁰₁,𝑥⁰₂,𝑦⁰,𝑓₃^{(𝑖𝑗𝑘)}=1𝜕𝑖!𝑗!𝑘!^{𝑖+𝑗+𝑘}𝑓₃𝜕^𝑖(𝑡)𝜕𝑥^𝑗₂𝑡−𝜏₂𝜕𝑦^𝑘𝑡−𝜏₂∣𝑥⁰₁,𝑥⁰₂,𝑦⁰,𝑓₁=𝑎𝑥₂(𝑡)−𝑟₁𝑥₁(𝑡)−𝑏𝑥₁(𝑡)−𝑐𝑥²₁(𝑡),𝑓₂=𝑏𝑥₁(𝑡)−𝑟₂𝑥₂(𝑡)−𝑏₁𝑥₂(𝑡)𝑥₂𝑡−𝜏₁−𝑎₁𝑥²₂(𝑡)𝑦(𝑡)1+𝑚𝑥²₂(𝑡),𝑓₃=𝑎₂𝑥²₂𝑡−𝜏₂𝑦𝑡−𝜏₂1+𝑚𝑥²₂𝑡−𝜏₂−𝑟𝑦(𝑡). (2.3) Then we can get the linearized system of system (2.2) 𝑑𝑥₁𝑑𝑡=𝛼₁₁𝑥₁(𝑡)+𝛼₁₂𝑥₂(𝑡),𝑑𝑥₂𝑑𝑡=𝛼₂₁𝑥₁(𝑡)+𝛼₂₂𝑥₂(𝑡)+𝛼₂₃𝑦(𝑡)+𝛽₂₂𝑥₂𝑡−𝜏₁,𝑑𝑦𝑑𝑡=𝛼₃₃𝑦(𝑡)+𝛾₃₂𝑥₂𝑡−𝜏₂+𝛾₃₃𝑦𝑡−𝜏₂. (2.4) Therefore, the corresponding characteristic equation of system (2.4) is 𝜆³+𝑚₂𝜆²+𝑚₁𝜆+𝑚₀+𝑛₂𝜆²+𝑛₁𝜆+𝑛₀𝑒^−𝜆𝜏1+𝑝₂𝜆²+𝑝₁𝜆+𝑝₀𝑒^−𝜆𝜏2+𝑞₁𝜆+𝑞₀𝑒^{−𝜆(𝜏1+𝜏2)}=0, (2.5) where 𝑚₀=(𝛼₁₂𝛼₂₁−𝛼₁₁𝛼₂₂)𝛼₃₃, 𝑚₁=𝛼₁₁𝛼₂₂+𝛼₁₁𝛼₃₃+𝛼₂₂𝛼₃₃−𝛼₁₂𝛼₂₁, 𝑚₂=−(𝛼₁₁+𝛼₂₂+𝛼₃₃), 𝑛₀=−𝛼₁₁𝛼₃₃𝛽₂₂, 𝑛₁=(𝛼₁₁+𝛼₃₃)𝛽₂₂, 𝑛₂=−𝛽₂₂, 𝑝₀=𝛼₁₁𝛼₂₃𝛾₃₂+𝛼₁₂𝛼₂₁𝛾₃₃−𝛼₁₁𝛼₂₂𝛾₃₃, 𝑝₁=𝛼₁₁𝛾₃₃+𝛼₂₂𝛾₃₃−𝛼₂₃𝛾₃₂, 𝑝₂=−𝛾₃₃,𝑞₀=−𝛼₁₁𝛽₂₂𝛾₃₃, 𝑞₁=𝛽₂₂𝛾₃₃.

Next, we consider the local stability of the positive equilibrium 𝐸⁰(𝑥⁰₁,𝑥⁰₂,𝑦⁰) and the Hopf bifurcation of system (1.4) for the different combination of 𝜏₁ and 𝜏₂.

Case 1. (𝜏₁=𝜏₂=0). The characteristic equation (2.5) becomes 𝜆³+𝑚₁₂𝜆²+𝑚₁₁𝜆+𝑚₁₀=0, (2.6) where 𝑚₁₂=𝑚₂+𝑛₂+𝑝₂, 𝑚₁₁=𝑚₁+𝑛₁+𝑝₁+𝑞₁, 𝑚₁₀=𝑚₀+𝑛₀+𝑝₀+𝑞₀.
It is not difficult to verify that 𝑚₁₂>0 and 𝑚₁₀>0. Thus, all the roots of (2.6) must have negative real parts, if the following condition holds: 𝐻₁₁∶𝑚₁₂𝑚₁₁>𝑚₁₀. Namely, the positive equilibrium 𝐸⁰(𝑥⁰₁,𝑥⁰₂,𝑦⁰) is locally stable in the absence of time delay, if 𝐻₁₁ holds.

Case 2. (𝜏₁>0,𝜏₂=0). On substituting 𝜏₂=0, (2.5) becomes 𝜆³+𝑚₂₂𝜆²+𝑚₂₁𝜆+𝑚₂₀+𝑛₂₂𝜆²+𝑛₂₁𝜆+𝑛₂₀𝑒^−𝜆𝜏1=0, (2.7) where 𝑚₂₂=𝑚₂+𝑝₂, 𝑚₂₁=𝑚₁+𝑝₁, 𝑚₂₀=𝑚₀+𝑝₀, 𝑛₂₂=𝑛₂, 𝑛₂₁=𝑛₁+𝑞₁, 𝑛₂₀=𝑛₀+𝑞₀.
Let 𝜆=𝑖𝜔₁(𝜔₁>0) be a root of (2.7). Then, we have 𝑛₂₁𝜔₁sin𝜏₁𝜔₁+𝑛₂₀−𝑛₂₂𝜔²₁cos𝜏₁𝜔₁=𝑚₂₂𝜔²₁−𝑚₂₀,𝑛₂₁𝜔₁cos𝜏₁𝜔₁−𝑛₂₀−𝑛₂₂𝜔²₁sin𝜏₁𝜔₁=𝜔³₁−𝑚₂₁𝜔₁. (2.8) Squaring both sides and adding them up, we get the following sixth-degree polynomial equation: 𝜔⁶₁+𝑒₂₂𝜔⁴₁+𝑒₂₁𝜔²₁+𝑒₂₀=0, (2.9) where 𝑒₂₂=𝑚²₂₂−𝑛²₂₂−2𝑚₂₁, 𝑒₂₁=𝑚²₂₁−2𝑚₂₀𝑚₂₂+2𝑛₂₀𝑛₂₂−𝑛²₂₁, 𝑒₂₀=𝑚²₂₀−𝑛²₂₀.
Let 𝜔²₁=𝑣₁, then (2.9) becomes 𝑣³₁+𝑒₂₂𝑣²₁+𝑒₂₁𝑣₁+𝑒₂₀=0. (2.10) Define 𝑓₁𝑣₁=𝑣³₁+𝑒₂₂𝑣²₁+𝑒₂₁𝑣₁+𝑒₂₀. (2.11) If 𝑒₂₀<0, it is easy to know that (2.10) has at least one positive root. On the other hand, if 𝑒₂₀≥0, according to Lemma 2.2 in [25], (2.10) has positive roots if 𝑒²₂₂−3𝑒₂₁>0 and 𝑣^∗₁=(−𝑒₂₂+𝑒²₂₂−3𝑒₂₁)/3>0, 𝑓₁(𝑣^∗₁)≤0 hold. Therefore, we give the following assumption.

𝐻₂₁: equation (2.10) has at least one positive root.

Without loss of generality, we assume that it has three positive roots which are denoted as 𝑣₁₁, 𝑣₁₂, and 𝑣₁₃. Thus, (2.9) has three positive roots 𝜔_1𝑘=√𝑣_1𝑘, 𝑘=1,2,3. The corresponding critical value of time delay 𝜏^(𝑗)_1𝑘 is 𝜏^(𝑗)_1𝑘=1𝜔_1𝑘𝐴arccos²⁴𝜔⁴_1𝑘+𝐴₂₂𝜔²_1𝑘+𝐴₂₀𝐵₂₄𝜔⁴_1𝑘+𝐵₂₂𝜔²_1𝑘+𝐵₂₀+2𝑗𝜋𝜔_1𝑘,𝑘=1,2,3,𝑗=0,1,2,…, (2.12) where 𝐴₂₄=𝑛₂₁−𝑚₂₂𝑛₂₂, 𝐴₂₂=𝑚₂₀𝑛₂₂+𝑚₂₂𝑛₂₀−𝑚₂₁𝑛₂₁, 𝐴₂₀=−𝑚₂₀𝑛₂₀, 𝐵₂₄=𝑛²₂₂, 𝐵₂₂=𝑛²₂₁−2𝑛₂₀𝑛₂₂, 𝐵₂₀=𝑛²₂₀.

Let 𝜏₁₀=min{𝜏⁽⁰⁾_1𝑘}, 𝑘∈{1,2,3}, 𝜔₁₀=𝜔_1𝑘0.

To verify the transversality condition of Hopf bifurcation, differentiating the two sides of (2.7) with respect to 𝜏₁, and noticing that 𝜆 is a function of 𝜏₁, we can obtain 𝑑𝜆𝑑𝜏₁⁻¹=−3𝜆²+2𝑚₂₂𝜆+𝑚₂₁𝜆𝜆³+𝑚₂₂𝜆²+𝑚₂₁𝜆+𝑚₂₀+2𝑛₂₂𝜆+𝑛₂₁𝜆𝑛₂₂𝜆²+𝑛₂₁𝜆+𝑛₂₀−𝜏₁𝜆. (2.13) Thus, Re𝑑𝜆𝑑𝜏₁⁻¹=Re−3𝜆²+2𝑚₂₂𝜆+𝑚₂₁𝜆𝜆³+𝑚₂₂𝜆²+𝑚₂₁𝜆+𝑚₂₀+Re2𝑛₂₂𝜆+𝑛₂₁𝜆𝑛₂₂𝜆²+𝑛₂₁𝜆+𝑛₂₀. (2.14) Therefore, Re𝑑𝜆𝑑𝜏₁⁻¹_{𝜆=𝑖𝜔10}=3𝜔⁴₁₀+2𝑚²₂₂−𝑛²₂₂−2𝑚₂₁𝜔²₁₀+𝑚²₂₁−2𝑚₂₀𝑚₂₂𝜔³₁₀−𝑚₂₁𝜔₁₀²+𝑚₂₀−𝑚₂₂𝜔²₁₀²−2𝑛²₂₂𝜔²₁₀+𝑛²₂₁−2𝑛₂₀𝑛₂₂𝑛₂₂𝜔²₁₀−𝑛₂₀²+𝑛²₂₁𝜔²₁₀. (2.15) From (2.9), we can get 𝜔³₁₀−𝑚₂₁𝜔₁₀²+𝑚₂₀−𝑚₂₂𝜔²₁₀²=𝑛₂₂𝜔²₁₀−𝑛₂₀²+𝑛²₂₁𝜔²₁₀. (2.16) Then, we have Re𝑑𝜆𝑑𝜏₁⁻¹_{𝜆=𝑖𝜔10}=3𝑣²_1∗+2𝑒₂₂𝑣_1∗+𝑒₂₁𝑛₂₂𝜔²₁₀−𝑛₂₀²+𝑛²₂₁𝜔²₁₀=𝑓^′₁𝑣_1∗𝑛₂₂𝜔²₁₀−𝑛₂₀²+𝑛²₂₁𝜔²₁₀, (2.17) where 𝑣_1∗=𝜔²₁₀∈{𝑣₁₁,𝑣₁₂,𝑣₁₃}.

Therefore, Re[𝑑𝜆/𝑑𝜏₁]⁻¹_{𝜆=𝑖𝜔10}≠0 if 𝐻₂₂∶𝑓^′₁(𝑣_1∗)≠0 holds. Notice that Re[𝑑𝜆/𝑑𝜏₁]⁻¹_{𝜆=𝑖𝜔10} and [𝑑Re(𝜆)/𝑑𝜏₁]_{𝜆=𝑖𝜔10} have the same sign. Then we have [𝑑Re(𝜆)/𝑑𝜏₁]_{𝜆=𝑖𝜔10}≠0 if 𝐻₂₂ holds. In conclusion, we have the following results.

Theorem 2.1. Suppose that the conditions 𝐻₂₁ and 𝐻₂₂ hold. The positive equilibrium 𝐸⁰ of system (1.4) is asymptotically stable for 𝜏₁∈[0,𝜏₁₀) and unstable when 𝜏₁>𝜏₁₀. Further, system (1.4) undergoes a Hopf bifurcation when 𝜏₁=𝜏₁₀.

Case 3. (𝜏₂>0,𝜏₁=0). On substituting 𝜏₂=0, (2.5) becomes 𝜆³+𝑚₃₂𝜆²+𝑚₃₁𝜆+𝑚₃₀+𝑝₃₂𝜆²+𝑝₃₁𝜆+𝑝₃₀𝑒^−𝜆𝜏2=0, (2.18) where 𝑚₃₂=𝑚₂+𝑛₂, 𝑚₃₁=𝑚₁+𝑛₁, 𝑚₃₀=𝑚₀+𝑛₀, 𝑝₃₂=𝑝₂, 𝑝₃₁=𝑝₁+𝑞₁, 𝑝₃₀=𝑝₀+𝑞₀.
Let 𝜆=𝑖𝜔₂(𝜔₂>0) be a root of (2.18). Then, we get 𝑝₃₁𝜔₂sin𝜏₂𝜔₂+𝑝₃₀−𝑝₃₂𝜔²₂cos𝜏₂𝜔₂=𝑚₃₂𝜔²₂−𝑚₃₀,𝑝₃₁𝜔₂cos𝜏₂𝜔₂−𝑝₃₀−𝑝₃₂𝜔²₂sin𝜏₂𝜔₂=𝜔³₂−𝑚₃₁𝜔₂, (2.19) which follows that 𝜔⁶₂+𝑒₃₂𝜔⁴₂+𝑒₃₁𝜔²₂+𝑒₃₀=0, (2.20) where 𝑒₃₂=𝑚²₃₂−𝑝²₃₂−2𝑚₃₁, 𝑒₃₁=𝑚²₃₁−2𝑚₃₀𝑚₃₂+2𝑝₃₀𝑝₃₂−𝑝²₃₁, 𝑒₃₀=𝑚²₃₀−𝑝²₃₀.
Let 𝜔²₂=𝑣₂, then (2.20) becomes 𝑣³₂+𝑒₃₂𝑣²₂+𝑒₃₁𝑣₂+𝑒₃₀=0. (2.21) Define 𝑓₂𝑣₂=𝑣³₂+𝑒₃₂𝑣²₂+𝑒₃₁𝑣₂+𝑒₃₀. (2.22) Similar as in case (2), we give the following assumption.
𝐻₃₁: equation (2.21) has at least one positive root.

Without loss of generality, we assume that it has three positive roots denoted by 𝑣₂₁, 𝑣₂₂, and 𝑣₂₃. Thus, (2.20) has three positive roots 𝜔_2𝑘=√𝑣_2𝑘,𝑘=1,2,3.

The corresponding critical value of time delay 𝜏^(𝑗)_2𝑘 is 𝜏^(𝑗)_2𝑘=1𝜔_2𝑘𝐴arccos³⁴𝜔⁴_2𝑘+𝐴₃₂𝜔²_2𝑘+𝐴₃₀𝐵₃₄𝜔⁴_2𝑘+𝐵₃₂𝜔²_2𝑘+𝐵₃₀+2𝑗𝜋𝜔_2𝑘,𝑘=1,2,3,𝑗=0,1,2,…, (2.23) where 𝐴₃₄=𝑝₃₁−𝑚₃₂𝑝₃₂, 𝐴₃₂=𝑚₃₀𝑝₃₂+𝑚₃₂𝑝₃₀−𝑚₃₁𝑝₃₁, 𝐴₃₀=−𝑚₃₀𝑝₃₀, 𝐵₃₄=𝑝²₃₂, 𝐵₃₂=𝑝²₃₁−2𝑝₃₀𝑝₃₂, 𝐵₃₀=𝑝²₃₀.

Let 𝜏₂₀=min{𝜏⁽⁰⁾_2𝑘}, 𝑘∈{1,2,3}, 𝜔₂₀=𝜔_2𝑘0.

Similar as in case (1), next, we suppose that the condition 𝐻₃₂∶𝑓^′₂(𝑣_2∗)≠0 holds, where 𝑣_2∗=𝜔²₂₀∈{𝑣₂₁,𝑣₂₂,𝑣₂₃}. Then we have [𝑑Re(𝜆)/𝑑𝜏₂]_{𝜆=𝑖𝜔20}≠0. By the above analysis, we have the following results.

Theorem 2.2. Suppose that the conditions 𝐻₃₁ and 𝐻₃₂ hold. The positive equilibrium 𝐸⁰ of system (1.4) is asymptotically stable for 𝜏₂∈[0,𝜏₂₀) and unstable when 𝜏₂>𝜏₂₀. Further, system (1.4) undergoes a Hopf bifurcation when 𝜏₂=𝜏₂₀.

Case 4. (𝜏₁=𝜏₂=𝜏>0).
For 𝜏₁=𝜏₂=𝜏>0, (2.5) can be rewritten in the following form: 𝜆³+𝑚₄₂𝜆²+𝑚₄₁𝜆+𝑚₄₀+𝑛₄₂𝜆²+𝑛₄₁𝜆+𝑛₄₀𝑒^−𝜆𝜏+𝑞₄₁+𝑞₄₀𝑒^−2𝜆𝜏=0, (2.24) where 𝑚₄₂=𝑚₂, 𝑚₄₁=𝑚₁, 𝑚₄₀=𝑚₀, 𝑛₄₂=𝑛₂+𝑝₂, 𝑛₄₁=𝑛₁+𝑝₁, 𝑛₄₀=𝑛₀+𝑝₀, 𝑞₄₁=𝑞₀, 𝑞₄₀=𝑞₀.
Multiplying 𝑒^𝜆𝜏 on both sides of (2.24), it is obvious to get 𝑛₄₂𝜆²+𝑛₄₁𝜆+𝑛₄₀+𝜆³+𝑚₄₂𝜆²+𝑚₄₁𝜆+𝑚₄₀𝑒^𝜆𝜏+𝑞₄₁𝜆+𝑞₄₀𝑒^−𝜆𝜏=0. (2.25) Let 𝜆=𝑖𝜔(𝜔>0) be the root of (2.25). Then, we have Δ₄₁sin𝜏𝜔+Δ₄₂cos𝜏𝜔=Δ₄₅,Δ₄₃cos𝜏𝜔+Δ₄₄sin𝜏𝜔=Δ₄₆, (2.26) where Δ₄₁=𝑞₄₁𝜔−𝑚₄₁𝜔+𝜔³, Δ₄₂=𝑚₄₀−𝑚₄₂𝜔²+𝑞₄₀, Δ₄₃=𝑞₄₁𝜔+𝑚₄₁𝜔−𝜔³, Δ₄₄=𝑚₄₀−𝑚₄₂𝜔²−𝑞₄₀, Δ₄₅=𝑛₄₂𝜔²−𝑛₄₀, Δ₄₆=−𝑛₄₁𝜔.
It follows that 𝐴sin𝜏𝜔=⁵𝜔⁵+𝐴₃𝜔³+𝐴₁𝜔𝜔⁶+𝐵₄𝜔⁴+𝐵₂𝜔²+𝐵₀,𝐴cos𝜏𝜔=⁴𝜔⁴+𝐴₂𝜔²+𝐴₀𝜔⁶+𝐵₄𝜔⁴+𝐵₂𝜔²+𝐵₀, (2.27) where 𝐴₀=(𝑞₄₀−𝑚₄₀)𝑛₄₀, 𝐴₁=(𝑚₄₁+𝑞₄₁)𝑛₄₀−(𝑚₄₀+𝑞₄₀)𝑛₄₁, 𝐴₂=(𝑚₄₀−𝑞₄₀)𝑛₄₂+(𝑞₄₁−𝑚₄₁)𝑛₄₁−𝑚₄₂𝑛₄₀, 𝐴₃=𝑚₄₂𝑛₄₁−𝑛₄₀−(𝑚₄₁+𝑞₄₁)𝑛₄₂, 𝐴₄=𝑛₄₁−𝑚₄₂𝑛₄₂, 𝐴₅=𝑛₄₂, 𝐵₀=𝑚²₄₀−𝑞²₄₀, 𝐵₂=𝑚²₄₁−𝑞²₄₁−2𝑚₄₀𝑚₄₂, 𝐵₄=𝑚²₄₂−2𝑚₄₁.
From (2.27), we can get 𝜔¹²+𝑒₄₅𝜔¹⁰+𝑒₄₄𝜔⁸+𝑒₄₃𝜔⁶+𝑒₄₂𝜔⁴+𝑒₄₁𝜔²+𝑒₄₀=0, (2.28) where 𝑒₄₅=2𝐵₄−𝐴²₅, 𝑒₄₄=𝐵²₄+2𝐵₂−𝐴²₄−2𝐴₃𝐴₅, 𝑒₄₃=2𝐵₀+2𝐵₂𝐵₄−𝐴²₃−2𝐴₁𝐴₅−2𝐴₂𝐴₄, 𝑒₄₂=𝐵²₂+2𝐵₀𝐵₄−𝐴²₂−2𝐴₀𝐴₄−2𝐴₁𝐴₃, 𝑒₄₁=2𝐵₀𝐵₂−𝐴²₁−2𝐴₀𝐴₂, 𝑒₄₀=𝐵²₀−𝐴²₀.
Let 𝜔²=𝑣₃, then (2.28) becomes 𝑣⁶₃+𝑒₄₅𝑣⁵₃+𝑒₄₄𝑣⁴₃+𝑒₄₃𝑣³₃+𝑒₄₂𝑣²₃+𝑒₄₁𝑣₃+𝑒₄₀=0. (2.29) Suppose that (2.29) has at least one positive root, and, without loss of generality, we assume that it has six positive roots which are denoted as 𝑣₃₁, 𝑣₃₂, 𝑣₃₃, 𝑣₃₄, 𝑣₃₅, and 𝑣₃₆. Then, (2.28) has six positive roots 𝜔_𝑘=√𝑣_3𝑘,𝑘=1,2,3,4,5,6.
The corresponding critical value of time delay 𝜏_𝑘^(𝑗) is 𝜏_𝑘^(𝑗)=1𝜔_𝑘𝐴arccos⁴𝜔⁴_𝑘+𝐴₂𝜔²_𝑘+𝐴₀𝜔⁶_𝑘+𝐵₄𝜔⁴_𝑘+𝐵₂𝜔²_2𝑘+𝐵₀+2𝑗𝜋𝜔_𝑘,𝑘=1,2,3,4,5,6,𝑗=0,1,2,…. (2.30) Let 𝜏₀=min{𝜏_𝑘⁽⁰⁾}, 𝑘∈{1,2,3,4,5,6}, 𝜔₀=𝜔_𝑘0.
Next, we verify the transversality condition. Differentiating (2.25) regarding 𝜏 and substituting 𝜏=𝜏₀, we get Re𝑑𝜆𝑑𝜏⁻¹_𝜏=𝜏0=Re𝐴+𝐵𝑖=𝐶+𝐷𝑖𝐴𝐶+𝐵𝐷𝐶²+𝐷², (2.31) where 𝐴=𝑚₄₁−3𝜔²₀cos𝜏₀𝜔₀−2𝑚₄₂𝜔₀sin𝜏₀𝜔₀+𝑞₄₁cos𝜏₀𝜔₀+𝑛₄₁,𝐵=𝑚₄₁−3𝜔²₀sin𝜏₀𝜔₀+2𝑚₄₂𝜔₀cos𝜏₀𝜔₀−𝑞₄₁sin𝜏₀𝜔₀+2𝑛₄₂𝜔₀,𝐶=𝑚₄₁−𝑞₄₁−𝜔²₀𝜔²₀cos𝜏₀𝜔₀+𝑞₄₀+𝑚₄₀−𝑚₄₂𝜔²₀𝜔₀sin𝜏₀𝜔₀,𝐷=𝑚₄₁+𝑞₄₁−𝜔²₀𝜔²₀sin𝜏₀𝜔₀+𝑞₄₀−𝑚₄₀+𝑚₄₂𝜔²₀𝜔₀cos𝜏₀𝜔₀. (2.32) Thus, if the condition (𝐻₄₂)∶𝐴𝐶+𝐵𝐷≠0 holds, the transversality condition is satisfied.

Theorem 2.3. Suppose that the conditions 𝐻₄₁ and 𝐻₄₂ hold. The positive equilibrium 𝐸⁰ of system (1.4) is asymptotically stable for 𝜏∈[0,𝜏₀) and unstable when 𝜏>𝜏₀. Further, system (1.4) undergoes a Hopf bifurcation when 𝜏=𝜏₀.

Case 5. (𝜏₁>0and𝜏₂∈[0,𝜏₂₀)).
We consider (2.5) with 𝜏₂ in its stable interval, and 𝜏₁ is considered as a parameter.
Let 𝜆=𝑖𝜔_1∗(𝜔_1∗>0) be the root of (2.5). Then, we have Δ₅₁sin𝜏₁𝜔_1∗+Δ₅₂cos𝜏₁𝜔_1∗=Δ₅₃,Δ₅₁cos𝜏₁𝜔_1∗−Δ₅₂sin𝜏₁𝜔_1∗=Δ₅₄, (2.33) where Δ₅₁=𝑛₁𝜔_1∗−𝑞₀sin𝜏₂𝜔_1∗+𝑞₁𝜔_1∗cos𝜏₂𝜔_1∗,Δ₅₂=𝑛₀−𝑛₂𝜔²_1∗+𝑞₀cos𝜏₂𝜔_1∗+𝑞₁𝜔_1∗sin𝜏₂𝜔_1∗,Δ₅₃=𝑚₂𝜔²_1∗−𝑚₀−𝑝₁𝜔_1∗sin𝜏₂𝜔_1∗+𝑝₂𝜔²_1∗−𝑝₀cos𝜏₂𝜔_1∗,Δ₅₄=𝜔³_1∗−𝑚₁𝜔_1∗−𝑝₁𝜔_1∗cos𝜏₂𝜔_1∗−𝑝₂𝜔²_1∗−𝑝₀sin𝜏₂𝜔_1∗. (2.34) From (2.33), we can get the following transcendental equation: 𝜔⁶_1∗+𝑒₅₂𝜔⁴_1∗+𝑒₅₁𝜔²_1∗+𝑒₅₀+𝑐₅₄𝜔⁴_1∗+𝑐₅₂𝜔²_1∗+𝑐₅₀cos𝜏₂𝜔_1∗+𝑐₅₅𝜔⁵_1∗+𝑐₅₃𝜔³_1∗+𝑐₅₁𝜔_1∗sin𝜏₂𝜔_1∗=0, (2.35) where 𝑒₅₀=𝑚²₀+𝑝²₀−𝑛²₀−𝑞²₀, 𝑒₅₁=𝑚²₁+𝑝²₁−𝑛²₁−𝑞²₁+2𝑛₀𝑛₂−2𝑚₀𝑚₂−2𝑝₀𝑝₂, 𝑒₅₂=𝑚²₂+𝑝²₂−𝑛²₂−2𝑚₁, 𝑐₅₀=2𝑚₀𝑝₀−2𝑛₀𝑞₀, 𝑐₅₁=2𝑚₀𝑝₁+2𝑛₁𝑞₀−2𝑚₁𝑝₀−2𝑛₀𝑞₁, 𝑐₅₂=2𝑚₁𝑝₁+2𝑛₂𝑞₀−2𝑚₂𝑝₀−2𝑚₀𝑝₂−2𝑛₁𝑞₁, 𝑐₅₃=2𝑚₁𝑝₂+2𝑛₂𝑞₁+2𝑝₀−2𝑚₂𝑝₁, 𝑐₅₄=2𝑚₂𝑝₂−2𝑝₁, 𝑐₅₅=−2𝑝₂.
In order to give the main results, we suppose that (2.35) has finite positive root. We denote the positive roots of (2.35) as 𝜔₅₁,𝜔₅₂⋯𝜔_5𝑘. For every 𝜔_5𝑖(𝑖=1,2,…,𝑘), the corresponding critical value of time delay 𝜏₁^(𝑗)_𝑖|𝑗=1,2,… is 𝜏₁^(𝑗)_𝑖=1𝜔_1∗Δarccos⁵¹Δ₅₄+Δ₅₂Δ₅₃Δ²₅₁+Δ²₅₂+2𝑗𝜋_{𝜔1∗=𝜔5𝑖},𝑖=1,2,…,𝑘,𝑗=0,1,2,…. (2.36) Let 𝜏^′_1∗=min{𝜏₁⁽⁰⁾_𝑖|𝑖=1,2,…,𝑘}, and 𝜔^′_1∗ is the corresponding root of (2.35) with 𝜏^′_1∗.
In the following, we differentiate the two sides of (2.5) with respect to 𝜏₁ to verify the transversality condition.
Taking the derivative of 𝜆 with respect to 𝜏₁ in (2.5) and substituting 𝜏₁=𝜏^′_1∗, we get Re𝑑𝜆𝑑𝜏_𝜏⁻¹_1=𝜏1∗𝑃=Re^𝑅+𝑃_𝐼𝑖𝑄_𝑅+𝑄_𝐼𝑖=𝑃_𝑅𝑄_𝑅+𝑃_𝐼𝑄_𝐼𝑄²_𝑅+𝑄²_𝐼, (2.37) where 𝑃_𝑅=𝑚₁−3𝜔^′_1∗²+2𝑛₂𝜔^′_1∗sin𝜏^′_1∗𝜔^′_1∗+𝑛₁cos𝜏^′_1∗𝜔^′_1∗+sin𝜏₂𝜔^′_1∗2𝑝₂𝜔^′_1∗−𝑝₁𝜏₂𝜔^′_1∗−𝑞₁sin𝜏^′_1∗𝜔^′_1∗+cos𝜏₂𝜔^′_1∗𝑝₂𝜏₂𝜔^′_1∗²+𝑝₁−𝑝₀𝜏₂+𝑞₁cos𝜏^′_1∗𝜔^′_1∗,𝑃_𝐼=2𝑚₂𝜔^′_1∗−𝑛₁sin𝜏^′_1∗𝜔^′_1∗+2𝑛₂𝜔^′_1∗cos𝜏^′_1∗𝜔^′_1∗+sin𝜏₂𝜔^′_1∗𝑝₀𝜏₂−𝑝₁−𝑝₂𝜏₂𝜔^′_1∗²−𝑞₁cos𝜏^′_1∗𝜔^′_1∗+cos𝜏₂𝜔^′_1∗2𝑝₂𝜔^′_1∗−𝑝₁𝜏₂𝜔^′_1∗−𝑞₁sin𝜏^′_1∗𝜔^′_1∗,𝑄_𝑅=𝑛₀𝜔^′_1∗−𝑛₂𝜔^′_1∗³sin𝜏^′_1∗𝜔^′_1∗−𝑛₁𝜔^′_1∗²cos𝜏^′_1∗𝜔^′_1∗+sin𝜏₂𝜔^′_1∗𝑞₀𝜔^′_1∗cos𝜏^′_1∗𝜔^′_1∗+𝑞₁𝜔^′_1∗²sin𝜏^′_1∗𝜔^′_1∗+cos𝜏₂𝜔^′_1∗𝑞₀𝜔^′_1∗sin𝜏^′_1∗𝜔^′_1∗−𝑞₁𝜔^′_1∗²cos𝜏^′_1∗𝜔^′_1∗,𝑄_𝐼=𝑛₀𝜔^′_1∗−𝑛₂𝜔^′_1∗³cos𝜏^′_1∗𝜔^′_1∗+𝑛₁𝜔^′_1∗²sin𝜏^′_1∗𝜔^′_1∗+sin𝜏₂𝜔^′_1∗𝑞₁𝜔^′_1∗²cos𝜏^′_1∗𝜔^′_1∗−𝑞₀𝜔^′_1∗sin𝜏^′_1∗𝜔^′_1∗+cos𝜏₂𝜔^′_1∗𝑞₁𝜔^′_1∗²sin𝜏^′_1∗𝜔^′_1∗+𝑞₀𝜔^′_1∗cos𝜏^′_1∗𝜔^′_1∗. (2.38) Obviously, if the condition 𝐻₅₂∶𝑃_𝑅𝑄_𝑅+𝑃_𝐼𝑄_𝐼≠0 holds, the transversality condition is satisfied. Through the above analysis, we have the following results.

Theorem 2.4. Suppose that the conditions 𝐻₅₁ and 𝐻₅₂ hold and 𝜏₂∈[0,𝜏₂₀). The positive equilibrium 𝐸⁰ of system (1.4) is asymptotically stable for 𝜏₁∈[0,𝜏^′_1∗) and unstable when 𝜏₁>𝜏^′_1∗. Further, system (1.4) undergoes a Hopf bifurcation when 𝜏₁=𝜏^′_1∗.

Case 6. (𝜏₂>0and𝜏₁∈[0,𝜏₁₀)).
We consider (2.5) with 𝜏₁ in its stable interval, and 𝜏₂ is considered as a parameter.
Substitute 𝜆=𝑖𝜔_2∗(𝜔_2∗>0) into (2.5). Then, we get 𝜔⁶_2∗+𝑒₆₂𝜔⁴_2∗+𝑒₆₁𝜔²_1∗+𝑒₆₀+𝑐₆₄𝜔⁴_2∗+𝑐₆₂𝜔²_2∗+𝑐₆₀cos𝜏₁𝜔_2∗+𝑐₆₅𝜔⁵_2∗+𝑐₆₃𝜔³_2∗+𝑐₆₁𝜔_2∗sin𝜏₁𝜔_2∗, (2.39) where 𝑒₆₀=𝑚²₀+𝑛²₀−𝑝²₀−𝑞²₀, 𝑒₆₁=𝑚²₁+𝑛²₁−𝑝²₁−𝑞²₁−2𝑛₀𝑛₂−2𝑚₀𝑚₂+2𝑝₀𝑝₂, 𝑒₆₂=𝑚²₂+𝑛²₂−𝑝²₂−2𝑚₁, 𝑐₆₀=2𝑚₀𝑛₀−2𝑝₀𝑞₀, 𝑐₆₁=2𝑚₀𝑛₁+2𝑝₁𝑞₀−2𝑚₁𝑛₀−2𝑝₀𝑞₁, 𝑐₆₂=2𝑚₁𝑛₁+2𝑝₂𝑞₀−2𝑚₂𝑛₀−2𝑚₀𝑛₂−2𝑝₁𝑞₁, 𝑐₆₃=2𝑚₁𝑛₂+2𝑝₂𝑞₁+2𝑛₀−2𝑚₂𝑛₁, 𝑐₆₄=2𝑚₂𝑛₂−2𝑛₁, 𝑐₆₅=−2𝑛₂.
Similar as in case (5), we give the following assumption. 𝐻₆₁: (2.39) has finite positive root.
The positive roots of (2.39) are denoted as 𝜔₆₁, 𝜔₆₂,…,𝜔_6𝑘. For every 𝜔_6𝑖(𝑖=1,2,…,𝑘), the corresponding critical value of time delay 𝜏₂^(𝑗)_𝑖∣𝑗=1,2,… is 𝜏₂^(𝑗)_𝑖=1𝜔_2∗Δarccos⁶¹Δ₆₄+Δ₆₂Δ₆₃Δ²₆₁+Δ²₆₂+2𝑗𝜋_{𝜔2∗=𝜔6𝑖},𝑘=1,2,3,4,5,6,𝑗=0,1,2,…, (2.40) where Δ₆₁=𝑝₁𝜔_2∗−𝑞₀sin𝜏₁𝜔_2∗+𝑞₁𝜔_2∗cos𝜏₁𝜔_2∗,Δ₆₂=𝑝₀−𝑝₂𝜔²_2∗+𝑞₀cos𝜏₁𝜔_2∗+𝑞₁𝜔_2∗sin𝜏₁𝜔_2∗,Δ₆₃=𝑚₂𝜔²_1∗−𝑚₀−𝑛₁𝜔_1∗sin𝜏₁𝜔_2∗+𝑛₂𝜔²_2∗−𝑛₀cos𝜏₁𝜔_2∗,Δ₆₄=𝜔³_2∗−𝑚₁𝜔_2∗−𝑛₁𝜔_2∗cos𝜏₁𝜔_2∗−𝑛₂𝜔²_2∗−𝑛₀sin𝜏₁𝜔_2∗. (2.41)(2.42)(2.43)(2.44)
Let 𝜏^′_2∗=min{𝜏₂⁽⁰⁾_𝑖∣𝑖=1,2,…,𝑘}, and 𝜔^′_2∗ is the corresponding root of (2.39) with 𝜏^′_2∗.
Then, we suppose that 𝐻₆₂∶[𝑑Re(𝜆)/𝑑𝜏₂]_{𝜆=𝑖𝜔′2∗} holds. By the general Hopf bifurcation theorem for FDEs in Hale [26], we have the following results.

Theorem 2.5. Suppose that the conditions 𝐻₆₁ and 𝐻₆₂ hold and 𝜏₁∈[0,𝜏₁₀). The positive equilibrium 𝐸⁰ of system (1.4) is asymptotically stable for 𝜏₂∈[0,𝜏^′_2∗) and unstable when 𝜏₂>𝜏^′_2∗. Further, system (1.4) undergoes a Hopf bifurcation at when 𝜏₂=𝜏^′_2∗.

3. Direction and Stability of Bifurcated Periodic Solutions

In Section 2, we have obtained the conditions under which a family of periodic solutions bifurcate from the positive equilibrium of system (1.4) when the delay crosses through the critical value. In this section, we will determine the direction of Hopf bifurcation and stability of bifurcating periodic solutions of system (1.4) with respect to 𝜏₁ for 𝜏₂∈(0,𝜏₂₀) by using the normal form method and center manifold theorem introduced by Hassard et al. [20]. It is considered that system (1.4) undergoes Hopf bifurcation at 𝜏₁=𝜏^′_1∗,𝜏₂∈(0,𝜏₂₀). Without loss of generality, we assume that 𝜏^′_1∗>𝜏^′_2∗, where 𝜏^′_1∗∈(0,𝜏₂₀).

Let 𝜏₁=𝜏^′_1∗+𝜇,𝜇∈𝑅, 𝑡=𝑠𝜏₁, 𝑥₁(𝑠𝜏₁)=𝑧₁(𝑠), 𝑥₂(𝑠𝜏₁)=𝑧₂(𝑠), 𝑦(𝑠𝜏₁)=𝑧₃(𝑠). We still denote 𝑠 by 𝑡. Then, system (1.4) can be transformed into the following system: 𝑢̇(𝑡)=𝐿_𝜇𝑢_𝑡+𝐹𝜇,𝑢_𝑡, (3.1) where 𝑢(𝑡)=(𝑢₁(𝑡),𝑢₂(𝑡),𝑢₃(𝑡))^𝑇∈𝐶=𝐶([−1,0],𝑅³) and 𝐿_𝜇∶𝐶→𝑅³, 𝐹∶𝑅×𝐶→𝑅³ are given, respectively, by 𝐿_𝜇𝜙=𝜏^′_1∗𝐴^′𝜙(0)+𝐵^′𝜙−𝜏^′_2∗𝜏^′_1∗+𝐶^′𝜙(),𝐹=𝜏−1𝜇,𝜙^′_1∗+𝜇𝐹₁,𝐹₂,𝐹₃^𝑇, (3.2) where 𝜙(𝜃)=𝜙₁(𝜃),𝜙₂(𝜃),𝜙₃(𝜃)^𝑇𝐴∈𝐶,_′=𝛼₁₁𝛼₁₂0𝛼₂₁𝛼₂₂𝛼₂₃00𝛼₃₃,𝐵^′=0000000𝛾₃₂𝛾₃₃,𝐶^′=0000𝛽₂₂0,𝐹000₁=𝑎₁₁𝜙²₁(0)𝐹+⋯,₂=𝑎₂₁𝜙²₂(0)+𝑎₂₂𝜙₂(0)𝜙₂(−1)+𝑎₂₃𝜙₂(0)𝜙₃(0)+𝑎₂₄𝜙³₂(0)+𝑎₂₅𝜙²₂(0)𝜙₃(0)𝐹+⋯,₃=𝑎₃₁𝜙²₂−𝜏^′_2∗𝜏^′_1∗+𝑎₃₂𝜙₂−𝜏^′_2∗𝜏^′_1∗𝜙₃−𝜏^′_2∗𝜏^′_1∗+𝑎₃₃𝜙³₂−𝜏^′_2∗𝜏^′_1∗+𝑎₃₄𝜙²₂−𝜏^′_2∗𝜏^′_1∗𝜙₃−𝜏^′_2∗𝜏^′_1∗𝑎+⋯,₁₁=−𝑐,𝑎₂₁=3𝑚𝑎₁𝑥⁰₂²𝑦⁰−𝑎₁𝑦⁰1+𝑚𝑥⁰₂²³,𝑎₂₂=−𝑏₁,𝑎₂₃=−2𝑎₁𝑥⁰₂1+𝑚𝑥⁰₂²²,𝑎₂₄=4𝑚𝑎₁𝑥⁰₂𝑦⁰1−𝑚𝑥⁰₂²1+𝑚𝑥⁰₂²⁴,𝑎₂₅=3𝑚𝑎₁𝑥⁰₂²−𝑎₁1+𝑚𝑥⁰₂²³,𝑎₃₁=𝑎₂𝑦⁰−3𝑚𝑎₂𝑥⁰₂²𝑦⁰1+𝑚𝑥⁰₂²³,𝑎₃₂=2𝑎₂𝑥⁰₂1+𝑚𝑥⁰₂²²,𝑎₃₃=4𝑚𝑎₂𝑥⁰₂𝑦⁰𝑚𝑥⁰₂²−11+𝑚𝑥⁰₂²⁴,𝑎₃₄=𝑎₂−3𝑚𝑎₂𝑥⁰₂²1+𝑚𝑥⁰₂²³. (3.3) Hence, by the Riesz representation theorem, there exists a 3×3 matrix function 𝜂(𝜃,𝜇)∶[−1,0]→𝑅³ whose elements are of bounded variation such that 𝐿_𝜇𝜙=⁰₋₁𝑑𝜂𝜃,𝜇𝜙(𝜃),𝜙∈𝐶[−1,0],𝑅³. (3.4) In fact, we choose 𝜂=𝜏𝜃,𝜇^′1∗+𝜇𝐴^′+𝐵^′+𝐶^′𝜏,𝜃=0,_′1∗+𝜇𝐵^′+𝐶^′−𝜏,𝜃∈^2∗𝜏₁,𝜏,0_′1∗+𝜇𝐶^′𝜏,𝜃∈−1,−^2∗𝜏₁,0,𝜃=−1. (3.5) For 𝜙∈𝐶([−1,0],𝑅³), we define 𝐴𝜇(𝜃)𝜙=𝑑𝜙∫𝑑𝜃,−1≤𝜃<0,₀₋₁𝑑𝜂𝜃,𝜇𝜙(𝜃),𝜃=0, (3.6) and 𝑅𝜇𝐹𝜙=0,−1≤𝜃<0,𝜇,𝜙,𝜃=0. (3.7) Then, system (3.1) can be transformed into the following operator equation: 𝑢̇(𝑡)=𝐴𝜇𝑢_𝑡+𝑅𝜇𝑢_𝑡. (3.8) For 𝜑∈𝐶^′([0,1],(𝑅³)^∗), where (𝑅³)^∗ is the 3-dimensional space of row vectors, we define the adjoint operator 𝐴^∗ of 𝐴(0): 𝐴^∗𝜑=−𝑑𝜑(𝑠)∫𝑑𝑠,0<𝑠≤1,₀₋₁𝑑^𝑇𝜂(𝜉,0)𝜑(−𝜉),𝑠=0. (3.9) For 𝜙∈𝐶([−1,0],𝑅³) and 𝜙∈𝐶^′([0,1],(𝑅³)^∗), we define a bilinear inner product: 𝜑(𝑠),𝜙(𝜃)=𝜑(0)𝜙(0)−⁰_𝜃=−1^𝜃_𝜉=0𝜑(𝜉−𝜃)𝑑𝜂(𝜃)𝜙(𝜉)𝑑𝜉, (3.10) where 𝜂(𝜃)=𝜂(𝜃,0).

By the discussion in Section 2, we know that ±𝑖𝜔^′_1∗𝜏^′_1∗ are eigenvalues of 𝐴(0). Thus, they are also eigenvalues of 𝐴^∗.

Suppose that 𝑞(𝜃)=(1,𝑞₂,𝑞₃)^𝑇𝑒^{𝑖𝜔′1∗𝜏′1∗𝜃} is the eigenvector of 𝐴(0) corresponding to 𝑖𝜔^′_1∗𝜏^′_1∗ and 𝑞^∗(𝑠)=(1/𝜌)(1,𝑞^∗₂,𝑞^∗₃)𝑒^{𝑖𝜔′1∗𝜏′1∗𝑠} is the eigenvector of corresponding to −𝑖𝜔^′_1∗𝜏^′_1∗. By direction computation, we can get 𝑞₂=𝑖𝜔^′_1∗−𝛼₁₁𝛼₁₂,𝑞₃=𝛾₃₂𝑖𝜔^′_1∗−𝛼₁₁𝑒^{−𝑖𝜏′2∗𝜔′1∗}𝛼₁₂𝑖𝜔^′_1∗−𝛼₃₃−𝛾₃₃𝑒^{−𝑖𝜏′2∗𝜔′1∗},𝑞^∗₂𝛼=−¹¹+𝑖𝜔^′_1∗𝛼₂₁,𝑞^∗₃=𝛼₂₃𝛼₁₁−𝑖𝜔^′_1∗𝛼₁₂𝛼₃₃+𝛾₃₃𝑒^{−𝑖𝜏′2∗𝜔′1∗}+𝑖𝜔^′_1∗. (3.11) Then, from (3.10), we can get 𝑞^∗(𝑠),𝑞(𝜃)=𝑞^∗(0)𝑞(0)−⁰_𝜃=−1^𝜃_𝜉=0𝑞^∗()(𝜃)𝑞(𝜉)=1𝜉−𝜃𝑑𝜂𝑑𝜉𝜌1+𝑞₂𝑞^∗₂+𝑞₃𝑞^∗₃−⁰₋₁1,𝑞^∗₂,𝑞^∗₃𝜃𝑒^{𝑖𝜏′1∗𝜔′1∗𝜃}𝑑𝜂(𝜃)1,𝑞₂,𝑞₃^𝑇=1𝜌1+𝑞₂𝑞^∗₂+𝑞₃𝑞^∗₃+𝜏^′_1∗𝛽₂₂𝑞₂𝑞^∗₂𝑒^{−𝑖𝜏′1∗𝜔′1∗}+𝜏^′_2∗𝑒^{−𝑖𝜏′2∗𝜔′1∗}𝛾₃₂𝑞₂+𝛾₃₃𝑞₃𝑞^∗₃. (3.12) Therefore, we can choose 𝜌=1+𝑞₂𝑞^∗₂+𝑞₃𝑞^∗₃+𝜏^′_1∗𝛽₂₂𝑞₂𝑞^∗₂𝑒^{−𝑖𝜏′1∗𝜔′1∗}+𝜏^′_2∗𝑒^{−𝑖𝜏′2∗𝜔′1∗}𝛾₃₂𝑞₂+𝛾₃₃𝑞₃𝑞^∗₃, (3.13) such that ⟨𝑞^∗,𝑞⟩=1 and ⟨𝑞^∗,𝑞⟩=0.

In the remainder of this section, following the algorithms given in [20] and using similar computation process in [27], we can get the coefficients that can be used to determine the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions:𝑔₂₀=2𝜏^′_1∗𝜌𝑎₁₁𝑞⁽¹⁾(0)²+𝑞^∗₂𝑎₂₁𝑞⁽²⁾(0)²+𝑎₂₂𝑞⁽²⁾(0)𝑞⁽²⁾(−1)+𝑎₂₃𝑞⁽²⁾(0)𝑞⁽³⁾(0)+𝑞^∗₃𝑎₃₁𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗²+𝑎₃₂𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗𝑞⁽³⁾−𝜏^′_2∗𝜏^′_1∗,𝑔₁₁=𝜏^′_1∗𝜌2𝑎₁₁𝑞⁽¹⁾(0)𝑞⁽¹⁾(0)+𝑞^∗₂2𝑎₂₁𝑞⁽²⁾(0)𝑞⁽²⁾(0)+𝑎₂₂𝑞⁽²⁾(0)𝑞⁽²⁾(−1)+𝑞⁽²⁾(0)𝑞⁽²⁾(−1)+𝑎₂₃𝑞⁽²⁾(0)𝑞⁽³⁾(0)+𝑞⁽²⁾(0)𝑞⁽³⁾(0)+𝑞^∗₃2𝑎₃₁𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗+𝑎₃₂𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗𝑞⁽³⁾−𝜏^′_2∗𝜏^′_1∗+𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗𝑞⁽³⁾−𝜏^′_2∗𝜏^′_1∗,𝑔₀₂=2𝜏^′_1∗𝜌𝑎₁₁𝑞⁽¹⁾(0)²+𝑞^∗₂𝑎₂₁𝑞⁽²⁾(0)²+𝑎₂₂𝑞⁽²⁾(0)𝑞⁽²⁾(−1)+𝑎₂₃𝑞⁽²⁾(0)𝑞⁽³⁾(0)+𝑞^∗₃𝑎₃₁𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗²+𝑎₃₂𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗𝑞⁽³⁾−𝜏^′_2∗𝜏^′_1∗,𝑔₂₁=2𝜏^′_1∗𝜌𝑎₁₁2𝑊⁽¹⁾₁₁(0)𝑞⁽¹⁾(0)+𝑊⁽¹⁾₂₀(0)𝑞⁽¹⁾(0)+𝑞^∗₂𝑎₂₁2𝑊⁽²⁾₁₁(0)𝑞⁽²⁾(0)+𝑊⁽²⁾₂₀(0)𝑞⁽²⁾(0)+𝑎₂₂𝑊⁽²⁾₁₁(0)𝑞⁽²⁾()+1−12𝑊⁽²⁾₂₀(0)𝑞⁽²⁾(−1)+𝑊⁽²⁾₁₁(−1)𝑞⁽²⁾(0)+12𝑊⁽²⁾₂₀(−1)𝑞⁽²⁾(0)+𝑎₂₃𝑊⁽²⁾₁₁(0)𝑞⁽³⁾(0)+12𝑊⁽²⁾₂₀(0)𝑞⁽³⁾(0)+𝑊⁽³⁾₁₁(0)𝑞⁽²⁾(0)+12𝑊⁽³⁾₂₀(−1)𝑞⁽²⁾(0)+3𝑎₂₄𝑞⁽²⁾(0)²𝑞⁽²⁾(0)+𝑎₂₅𝑞⁽²⁾(0)²𝑞⁽³⁾(0)+2𝑞⁽²⁾(0)𝑞⁽²⁾(0)𝑞⁽³⁾(0)+𝑞^∗₃𝑎₃₁2𝑊⁽²⁾₁₁−𝜏^′_2∗𝜏^′_1∗𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗+𝑊⁽²⁾₂₀−𝜏^′_2∗𝜏^′_1∗𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗+𝑎₃₂𝑊⁽²⁾₁₁−𝜏^′_2∗𝜏^′_1∗𝑞⁽³⁾−𝜏^′_2∗𝜏^′_1∗+12𝑊⁽²⁾₂₀−𝜏^′_2∗𝜏^′_1∗𝑞⁽³⁾−𝜏^′_2∗𝜏^′_1∗+𝑊⁽³⁾₁₁−𝜏^′_2∗𝜏^′_1∗𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗+12𝑊⁽³⁾₂₀−𝜏^′_2∗𝜏^′_1∗𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗+3𝑎₃₃𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗²𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗+𝑎₃₄𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗²𝑞⁽³⁾−𝜏^′_2∗𝜏^′_1∗+2𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗𝑞⁽³⁾−𝜏^′_2∗𝜏^′_1∗, (3.14)(3.15) with 𝑊₂₀(𝜃)=𝑖𝑔₂₀𝑞(0)𝜏^′_1∗𝜔^′_1∗𝑒^{𝑖𝜏′1∗𝜔′1∗𝜃}+𝑖𝑔₀₂𝑞(0)3𝜏^′_1∗𝜔^′_1∗𝑒^{−𝑖𝜏′1∗𝜔′1∗𝜃}+𝐸₁𝑒^{2𝑖𝜏′1∗𝜔′1∗𝜃},𝑊₁₁(𝜃)=𝑔₁₁𝑞(0)𝑖𝜏^′_1∗𝜔^′_1∗𝑒^{𝑖𝜏′1∗𝜔′1∗𝜃}−𝑔₁₁𝑞(0)𝑖𝜏^′_1∗𝜔^′_1∗𝑒^{−𝑖𝜏′1∗𝜔′1∗𝜃}+𝐸₂, (3.16) where 𝐸₁, 𝐸₂ can be determined by the following equations, respectively, 𝛼^′₁₁−𝛼₁₂0−𝛼₂₁𝛼^′₂₂−𝛼₂₃0−𝛾₃₂𝑒^{−2𝑖𝜏′2∗𝜔′1∗}𝛼^′₃₃𝐸₁Δ=2Δ¹⁽¹⁾₁⁽²⁾Δ₁⁽³⁾,𝛼₁₁𝛼₁₂0𝛼₂₁𝛼₂₂+𝛽₂₂𝛼₂₃0𝛾₃₂𝛼₃₃+𝛾₃₃𝐸₂Δ=−Δ²⁽¹⁾₂⁽²⁾Δ₂⁽³⁾, (3.17) with Δ₁⁽¹⁾=𝑎₁₁𝑞⁽¹⁾(0)²,Δ₁⁽²⁾=𝑎₂₁𝑞⁽²⁾(0)²+𝑎₂₂𝑞⁽²⁾(0)𝑞⁽²⁾(−1)+𝑎₂₃𝑞⁽²⁾(0)𝑞⁽³⁾(0),Δ₁⁽³⁾=𝑎₃₁𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗²+𝑎₃₂𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗𝑞⁽³⁾−𝜏^′_2∗𝜏^′_1∗,Δ₂⁽¹⁾=2𝑎₁₁𝑞⁽¹⁾(0)𝑞⁽¹⁾(0),Δ₂⁽²⁾=2𝑎₂₁𝑞⁽²⁾(0)𝑞⁽²⁾(0)+𝑎₂₂𝑞⁽²⁾(0)𝑞⁽²⁾(−1)+𝑞⁽²⁾(0)𝑞⁽²⁾(−1)+𝑎₂₃𝑞⁽²⁾(0)𝑞⁽³⁾(0)+𝑞⁽²⁾(0)𝑞⁽³⁾(0),Δ₂⁽³⁾=2𝑎₃₁𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗+𝑎₃₂𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗𝑞⁽³⁾−𝜏^′_2∗𝜏^′_1∗+𝑞⁽²⁾−𝜏^′_2∗𝜏^′_1∗𝑞⁽³⁾−𝜏^′_2∗𝜏^′_1∗,𝛼^′₁₁=2𝑖𝜔^′_1∗−𝛼₁₁,𝛼^′₂₂=2𝑖𝜔^′_1∗−𝛼₂₂−𝛽₂₂𝑒^{−2𝑖𝜏′1∗𝜔′1∗},𝛼^′₃₃=2𝑖𝜔^′_1∗−𝛼₃₃−𝛾₃₃𝑒^{−2𝑖𝜏′2∗𝜔′1∗}. (3.18)(3.19) Then, we can calculate the following values: 𝐶₁(0)=𝑖2𝜏^′_1∗𝜔^′_1∗𝑔₁₁𝑔₂₀−2𝑔₁₁²−𝑔₀₂²3+𝑔₂₁2,𝐶𝛿=−Re₁(0)Re𝜆^′𝜏^′_1∗,𝜎=2Re𝐶₁(0,𝐶)𝑇=−Im₁(0)+𝛿Im𝜆^′𝜏^′_1∗𝜏^′_1∗𝜔^′_1∗. (3.20)

Based on the above discussion, we can obtain the following results.

Theorem 3.1. From (3.20) one has(i)the direction of the Hopf bifurcation is determined by the sign of 𝛿: if 𝛿>0(𝛿<0), then the Hopf bifurcation is supercritical (subcritical);(ii)the stability of bifurcating periodic solutions is determined by the sign of 𝜎: if 𝜎<0(𝜎>0), the bifurcating periodic solutions are stable (unstable);(iii)the period of the bifurcating periodic solution is determined by the sign of 𝑇: if 𝑇>0(𝑇<0), the bifurcating periodic solution increases (decreases).

4. Numerical Example

In this section, we give some numerical simulations to verify the theoretical analysis in Sections 2 and 3. Let 𝑎=8, 𝑎₁=4.25, 𝑎₂=3, 𝑏=5, 𝑏₁=1, 𝑐=0.5, 𝑚=2, 𝑟=1, 𝑟₁=1, and 𝑟₂=2. Then, we have the following particular case of system (1.4): 𝑑𝑥₁𝑑𝑡=8𝑥₂(𝑡)−𝑥₁(𝑡)−5𝑥₁(𝑡)−0.5𝑥²₁(𝑡),𝑑𝑥₂𝑑𝑡=5𝑥₁(𝑡)−2𝑥₂(𝑡)−𝑥₂(𝑡)𝑥₂𝑡−𝜏₁−4.25𝑥²₂(𝑡)𝑦(𝑡)1+2𝑥²₂(𝑡),𝑑𝑦=𝑑𝑡3𝑥²₂𝑡−𝜏₂𝑦𝑡−𝜏₂1+2𝑥²₂𝑡−𝜏₂−𝑦(𝑡). (4.1) It is not difficult to verify that 𝑎₂>𝑚𝑟, 𝑏𝑥⁰₁>(𝑟₂−𝑏₁𝑥⁰₂)𝑥⁰₂, namely, the conditions 𝐻₁and 𝐻₂ hold. Therefore, system (4.1) has at least a positive equilibrium. By means of Matlab, we can get that the positive equilibrium of (4.1) is 𝐸⁰_∗(1.2111,1.0000,2.1568).

For 𝜏₁>0, 𝜏₂=0, we can get 𝜔₁₀=1.3881, 𝜏₁₀=0.9032. From Theorem 2.2, we know that the positive equilibrium 𝐸⁰_∗ is asymptotically stable when 𝜏₁∈[0,𝜏₁₀). The corresponding waveform and the phase plot are illustrated by Figure 1. When the time delay 𝜏₁ passes through the critical value 𝜏₁₀, the positive equilibrium 𝐸⁰_∗ will lose its stability and a Hopf bifurcation occurs, and a family of periodic solutions bifurcate from the positive equilibrium 𝐸⁰_∗. This property is illustrated by the numerical simulation in Figure 2. Similarly, we have 𝜔₂₀=0.8497, 𝜏₂₀=0.5124, when 𝜏₂>0, 𝜏₁=0. The corresponding waveform and the phase plots are shown in Figures 3 and 4.

Figure 1 

𝐸0 is locally asymptotically stable for 𝜏1=0.8500<𝜏10=0.9032 with initial value “2.31; 1.9; 2.34.”

Figure 2 

𝐸0 is unstable for 𝜏1=0.9200>𝜏10=0.9032 with initial value “2.31; 1.9; 2.34.”

Figure 3 

𝐸0 is locally asymptotically stable for 𝜏2=0.4400<𝜏20=0.5124 with initial value “2.31; 1.9; 2.34.”

Figure 4 

𝐸0 is unstable for 𝜏2=0.5500>𝜏20=0.5124 with initial value “2.31; 1.9; 2.34.”

For 𝜏₁=𝜏₂=𝜏>0, we can obtain 𝜔₀=1.0000, 𝜏₀=0.4178. From Theorem 2.2, we know that, when the time delay 𝜏 increases from zero to 𝜏₀, the positive equilibrium 𝐸⁰_∗ is asymptotically stable. Once the time delay 𝜏 passes through the critical value 𝜏₀, the positive equilibrium 𝐸⁰_∗ will lose its stability and a Hopf bifurcation occurs. This property is illustrated by the numerical simulation in Figures 5 and 6.

Figure 5 

𝐸0 is locally asymptotically stable for 𝜏=0.3900<𝜏0=0.4178 with initial value “2.31; 1.9; 2.34.”

Figure 6 

𝐸0 is unstable for 𝜏=0.4600>𝜏0=0.4178 with initial value “2.31; 1.9; 2.34.”

For 𝜏₁>0 and 𝜏^′_2∗=0.25∈[0,𝜏₂₀), we have 𝜔^′_1∗=0.8020, 𝜏^′_1∗=0.6125. According to Theorem 2.2, the positive equilibrium 𝐸⁰_∗ is asymptotically stable when 𝜏₁∈[0,𝜏^′_1∗) and unstable when 𝜏₁>𝜏^′_1∗, which can be depicted by the numerical simulation in Figures 7 and 8. In addition, from (3.20), we can obtain 𝐶₁(0)=−1.1949+3.2464𝑖, 𝛿=−23.4294, 𝜎=−2.3898, 𝑇=2.4949. Thus, from Theorem 2.3, we know that the Hopf bifurcation with respect to 𝜏₁ with 𝜏^′_2∗=0.25∈[0,𝜏₂₀) is subcritical, the bifurcating periodic solutions are stable and increase. Similarly, we have 𝜔^′_2∗=0.8872, 𝜏^′_2∗=0.5094, for 𝜏₂>0 and 𝜏^′_1∗=0.15∈[0,𝜏₁₀). The corresponding waveform and the phase plots are shown in Figures 9 and 10.

Figure 7 

𝐸0 is locally asymptotically stable for 𝜏1=0.6000<𝜏′1∗=0.6125 and 𝜏′2∗=0.25 with initial value “2.31; 1.9; 2.34.”

Figure 8 

𝐸0 is unstable for 𝜏1=0.6800>𝜏′1∗=0.6125 and 𝜏′2∗=0.25 with initial value “2.31; 1.9; 2.34.”

Figure 9 

𝐸0 is locally asymptotically stable for 𝜏2=0.4500<𝜏′2∗=0.5094 and 𝜏′1∗=0.15 with initial value “2.31; 1.9; 2.34.”

Figure 10 

𝐸0 is unstable for 𝜏2=0.5700>𝜏′2∗=0.5094 and 𝜏′1∗=0.15 with initial value “2.31; 1.9; 2.34.”

5. Conclusions

In this paper, a delayed predator-prey system with Holling type III functional response and stage structure for the prey population is investigated. Compared with literature [14], we consider not only the time delay due to the gestation of the predator but also the negative feedback of the mature prey density and the intraspecific competition of the immature prey population. F. Li and H. W. Li [14] has obtained that the species in system (4.1) with only the time delay due to the gestation of the predator could coexist. However, we get that the species could also coexist with some available time delays of the mature prey and the predator. This is valuable from the view of ecology.

The sufficient conditions for the local stability of the positive equilibrium and the existence of local Hopf bifurcation for the possible combinations of two delays are obtained. The main results are given in Theorems 2.1–2.5. By computation, we find that the time delay due to the gestation of the predator is marked because the critical value of 𝜏₂ is smaller than that of 𝜏₁ when we only consider them, respectively. Furthermore, the explicit formulae which determines the direction of the bifurcation and the stability of the bifurcating periodic solutions is established when 𝜏>0 and 𝜏₂∈[0,𝜏₂₀) by using the normal form theory and center manifold theorem. The main results are given in Theorem 2.3. Finally, numerical simulations are carried out to support the obtained theoretical results.

Acknowledgments

The authors are grateful to the anonymous reviewers for their helpful comments and valuable suggestions on improving the paper. This work was supported by the National Natural Science Foundation of China (61273070), Doctor Candidate Foundation of Jiangnan University (JUDCF12030) and Anhui Provincial Natural Science Foundation under Grant no. 1208085QA11.