Stability and Hopf Bifurcation Analysis for a Gause-Type Predator-Prey System with Multiple Delays

Abstract

This paper is concerned with a Gause-type predator-prey system with two delays. Firstly, we study the stability and the existence of Hopf bifurcation at the coexistence equilibrium by analyzing the distribution of the roots of the associated characteristic equation. A group of sufficient conditions for the existence of Hopf bifurcation is obtained. Secondly, an explicit formula for determining the stability and the direction of periodic solutions that bifurcate from Hopf bifurcation is derived by using the normal form theory and center manifold argument. Finally, some numerical simulations are carried out to illustrate the main theoretical results.

1. Introduction

Multispecies predator-prey models have been studied by many scholars [1–7]. Guo and Jiang [7] studied the following three-species food-chain system: 𝑑𝑥(𝑡)(𝑡)𝑥(𝑡)𝑑𝑡=𝛼𝑥1−𝑘−𝛽𝑥(𝑡−𝜏)𝑦(𝑡)(),(𝑡)1+𝑝𝑥𝑡−𝜏𝑑𝑦(𝑡)(𝑡)𝑑𝑡=𝑦−ℎ+𝑒𝛽𝑥(𝑡)(𝑡)𝑧(𝑡),(𝑡)1+𝑝𝑥−𝑟𝑦𝑑𝑧(𝑡)(𝑡),𝑑𝑡=𝑧−𝑠+𝑚𝑟𝑦 (1) where 𝑥(𝑡), 𝑦(𝑡), and 𝑧(𝑡) are the population densities of the prey, the predator and the top predator at time 𝑡. The prey grows with intrinsic growth rate 𝛼 and carrying capacity 𝑘 in the absence of predation. The predator captures the prey with capture rate 𝛽 and Holling type II functional response 𝑥/(1+𝑝𝑥). The top predator captures its prey (the predator) with capture rate 𝑟 and Holling type I functional response 𝑟𝑦. The predator and the top predator contribute to their growth with the conversion rates 𝑒 and 𝑚, respectively. The parameters ℎ and 𝑠 are the death rates of the predator and the top predator, respectively. All the parameters 𝛼, 𝛽, 𝑒, ℎ, 𝑘, 𝑚, 𝑝, 𝑟, 𝑠 and in system (1) are assumed to be positive. The constant 𝜏≥0 represents the time delay due to the gestation of the prey. Guo and Jiang [7] investigated the bifurcation phenomenon and the properties of periodic solutions of system (1).

Predator-prey systems with single delay as system (1) have been investigated extensively [8–12]. However, there are some papers on the bifurcations of a population dynamics with multiple delays [13–16]. Gakkhar and Singh [15] studied the effects of two delays on a delayed predator-prey system with modified Leslie-Gower and Holling type II functional response and established the existence of periodic solutions via Hopf bifurcation with respect to both delays. Motivated by the work of Guo and Jiang [7] and Gakkhar and Singh [15], we consider the following predator-prey system with two delays: 𝑑𝑥(𝑡)(𝑡)𝑥(𝑡)𝑑𝑡=𝛼𝑥1−𝑘−𝛽𝑥(𝑡)𝑦(𝑡)(𝑡),(𝑡)1+𝑝𝑥𝑑𝑦(𝑡)𝑑𝑡=𝑦−ℎ+𝑒𝛽𝑥𝑡−𝜏₁1+𝑝𝑥𝑡−𝜏₁(𝑡)𝑧(𝑡),(𝑡)−𝑟𝑦𝑑𝑧(𝑡)𝑑𝑡=𝑧−𝑠+𝑚𝑟𝑦𝑡−𝜏₂, (2) where 𝜏₁ denotes the time delay due to the gestation of the predator and 𝜏₂ denotes the time delay due to the gestation of the top predator.

This paper is organized as follows. In the next section, we will consider the stability of the positive equilibrium of system (2) and the existence of local Hopf bifurcation at the positive equilibrium. In Section 3, we can determine the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions from the Hopf bifurcation. Some numerical simulations are also given to illustrate the theoretical prediction in Section 4.

2. Local Stability and Hopf Bifurcation

Because we are only interested in the case in which the species can coexist, then we only consider the positive equilibrium of system (2). It is not difficult to know that if conditions (𝐻₁)∶𝛽𝑠<𝛼𝑚𝑟 and (𝐻₂)∶𝑒𝛽𝑥_∗>ℎ(1+𝑝𝑥_∗) hold, then system (2) has a unique positive equilibrium 𝐸_∗(𝑥_∗,𝑦_∗,𝑧_∗), where 𝑥_∗=−𝐵+𝐵²−4𝐴𝐶2,𝑦_∗=𝑠𝑚𝑟,𝑧_∗=1𝑟𝑒𝛽𝑥_∗1+𝑝𝑥_∗−ℎ, (3) where 𝐴=𝛼𝑚𝑝𝑟,𝐵=𝛼𝑚𝑟1−𝑝𝑘,𝐶=𝑘𝛽𝑠−𝛼𝑚𝑟. (4)

Let 𝑥(𝑡)=𝑧₁(𝑡)+𝑥_∗, 𝑦(𝑡)=𝑧₂(𝑡)+𝑦_∗, and 𝑧(𝑡)=𝑧₃(𝑡)+𝑧_∗ and still denote 𝑧₁(𝑡), 𝑧₂(𝑡), and 𝑧₃(𝑡) by 𝑥(𝑡), 𝑦(𝑡), and 𝑧(𝑡) respectively. Then system (2) can be transformed to the following form: 𝑑𝑥(𝑡)𝑑𝑡=𝑎₁₁𝑥(𝑡)+𝑎₁₂𝑦(𝑡)+𝑓₁,𝑑𝑦(𝑡)𝑑𝑡=𝑎₂₃𝑧(𝑡)+𝑏₁₁𝑥𝑡−𝜏₁+𝑓₂,𝑑𝑧(𝑡)𝑑𝑡=𝑐₃₂𝑦𝑡−𝜏₂+𝑓₃, (5) where 𝑎₁₁=𝛼−2𝛼𝑘𝑥_∗−𝛽𝑦_∗1+𝑝𝑥_∗²,𝑎₁₂=−𝛽𝑥_∗1+𝑝𝑥_∗,𝑎₂₃=−𝑟𝑦_∗,𝑏₁₁=𝛽𝑒𝑦_∗1+𝑝𝑥_∗²,𝑐₃₂=𝑚𝑟𝑧_∗,𝑓₁=𝑎₁₃𝑥²(𝑡)+𝑎₁₄𝑥(𝑡)𝑦(𝑡)+𝑎₁₅𝑥³(𝑡)+𝑎₁₆𝑥²(𝑡)𝑦(𝑡)𝑓+⋯,₂=𝑎₂₄𝑦(𝑡)𝑧(𝑡)+𝑎₂₅𝑥𝑡−𝜏₁𝑦(𝑡)+𝑎₂₆𝑥²𝑡−𝜏₁+𝑎₂₇𝑥²𝑡−𝜏₁𝑦(𝑡)+𝑎₂₈𝑥³𝑡−𝜏₁𝑓+⋯,₃=𝑐₃₃𝑦𝑡−𝜏₂𝑧(𝑡),𝑎₁₃𝛼=−𝑘+𝛽𝑝𝑦_∗1+𝑝𝑥_∗³,𝑎₁₄𝛽=−1+𝑝𝑥_∗²,𝑎₁₅=−𝛽𝑝²𝑦_∗1+𝑝𝑥_∗⁴,𝑎₁₆=𝛽𝑝1+𝑝𝑥_∗³,𝑎₂₄=−𝑟,𝑎₂₅=𝛽𝑒1+𝑝𝑥_∗²,𝑎₂₆=−𝛽𝑒𝑝𝑦_∗1+𝑝𝑥_∗³,𝑎₂₇=−𝛽𝑒𝑝1+𝑝𝑥_∗³,𝑎₂₈=𝛽𝑒𝑝²𝑦_∗1+𝑝𝑥_∗⁴,𝑐₃₃=𝑚𝑟. (6) The linearized system of system (5) is 𝑑𝑥(𝑡)𝑑𝑡=𝑎₁₁𝑥(𝑡)+𝑎₁₂𝑦(𝑡),𝑑𝑦(𝑡)𝑑𝑡=𝑎₂₃𝑧(𝑡)+𝑏₁₁𝑥𝑡−𝜏₁,𝑑𝑧(𝑡)𝑑𝑡=𝑐₃₂𝑦𝑡−𝜏₂. (7) Then the associated characteristic equation of system (7) at the origin is of the form 𝜆³+𝑝₂𝜆²+𝑚₁𝜆𝑒^−𝜆𝜏1+𝑛_𝜆+𝑛₀𝑒^−𝜆𝜏2=0, (8) where 𝑝₂=−𝑎₁₁,𝑚₁=−𝑎₁₂𝑏₁₁,𝑛₁=−𝑎₂₃𝑐₃₂,𝑛₀=𝑎₁₁𝑎₂₃𝑐₃₂. (9)

Case 1. One has 𝜏₁=𝜏₂=𝜏=0.
Equation (8) becomes 𝜆³+𝑝₂𝜆²+𝑚₁+𝑚₁𝜆+𝑛₀=0. (10) Obviously, if conditions (𝐻₁₁)∶𝑝₂>0 and (𝐻₁₂)∶𝑝₂(𝑚₁+𝑛₁)>𝑛₀ hold, then all the roots of (10) must have negative real parts. Then, we can conclude that the positive equilibrium 𝐸_∗(𝑥_∗,𝑦_∗,𝑧_∗) is locally asymptotically stable in the absence of delay.

Case 2. One has 𝜏₁>0,𝜏₂=0.
Equation (8) becomes 𝜆³+𝑝₂𝜆²+𝑛₁𝜆+𝑛₀+𝑚₁𝜆𝑒^−𝜆𝜏1=0. (11) Letting 𝜆=𝑖𝜔₁(𝜔₁>0) be a root of (11), then we have 𝑚₁𝜔₁sin𝜏₁𝜔₁=𝑝₂𝜔²₁−𝑛₀,𝑚₁𝜔₁cos𝜏₁𝜔₁=𝜔³₁−𝑛₁𝜔₁. (12) It follows that 𝜔⁶₁+𝑐₂₂𝜔⁴₁+𝑐₂₁𝜔²₁+𝑐₂₀=0, (13) where 𝑐₂₂=𝑝²₂−2𝑛₁,𝑐₂₁=𝑛²₁−𝑚²₁−2𝑛₀𝑝₂,𝑐₂₀=𝑛²₀. (14) Letting 𝜔²₁=𝑣₁, then (13) becomes 𝑣³₁+𝑐₂₂𝑣²₁+𝑐₂₁𝑣₁+𝑐₂₀=0. (15) Obviously, 𝑐₂₀≥0. Thus, we assume that (15) has at least one positive solution. Without loss of generality, we assume that it has three positive roots, which are denoted as 𝑣₁₁, 𝑣₁₂, and 𝑣₁₃. Then (13) has three positive roots 𝜔_1𝑘=√𝑣_1𝑘, 𝑘=1,2,3.

From (12), we can get 𝜏^(𝑗)_1𝑘=1𝜔_1𝑘𝜔arccos^21𝑘−𝑛₁𝑚₁+2𝑗𝜋𝜔_1𝑘,𝑗=0,1,…,𝑘=1,2,3. (16) Then, we denote 𝜏₁₀=min𝜏⁽⁰⁾_1𝑘,𝑘=1,2,3,𝜔₁₀=𝜔_1𝑘0. (17) Next, we verify the transversality condition. Differentiating the two sides of (11) with respect to 𝜏₁ and noticing that 𝜆 is a function of 𝜏₁, we can get 𝑑𝜆𝑑𝜏₁⁻¹=−3𝜆²+2𝑝₂𝜆+𝑛₁𝜆𝜆³+𝑝₂𝜆²+𝑛₁𝜆+𝑛₀+1𝜆²−𝜏₁𝜆. (18) Therefore Re𝑑𝜆𝑑𝜏₁⁻¹_{𝜆=𝑖𝜔10}=3𝜔⁴₁₀+2𝑝²₂−2𝑛₁𝜔²₁₀+𝑛²₁−2𝑛₀𝑝₂𝜔³₁₀−𝑛₁𝜔₁₀²+𝑛₀−𝑝₂𝜔²₁₀²−𝑚²₁𝑚²₁𝜔²₁₀. (19) From (13), we have 𝜔³₁₀−𝑛₁𝜔₁₀²+𝑛₀−𝑝₂𝜔²₁₀²=𝑚²₁𝜔²₁₀. (20) Therefore Re𝑑𝜆𝑑𝜏₁⁻¹_{𝜆=𝑖𝜔10}=𝑓^′₁𝑣_1∗𝜔³₁₀−𝑛₁𝜔₁₀²+𝑛₀−𝑝₂𝜔²₁₀², (21) with 𝑓₁𝑣₁=𝑣³₁+𝑐₂₂𝑣²₁+𝑐₂₁𝑣₁+𝑐₂₀,𝑣_1∗=𝜔²₁₀. (22) Obviously, if (𝐻₂₂)∶𝑓^′₁(𝑣_1∗)≠0, then Re[𝑑𝜆/𝑑𝜏₁]⁻¹_{𝜆=𝑖𝜔10}≠0. Thus, if condition (𝐻₂₂) holds, the transversality condition is satisfied. In conclusion, we have the following results.

Theorem 1. Suppose that conditions (𝐻₂₁)-(𝐻₂₂) hold. The positive equilibrium 𝐸_∗(𝑥_∗,𝑦_∗,𝑧_∗) of system (2) is asymptotically stable for 𝜏₁∈[0,𝜏₁₀) and unstable when 𝜏₁>𝜏₁₀. And system (2) undergoes a Hopf bifurcation at 𝐸_∗(𝑥_∗,𝑦_∗,𝑧_∗) when 𝜏₁=𝜏₁₀.

Case 3. One has 𝜏₁=0,𝜏₂>0.
Equation (8) becomes 𝜆³+𝑝₂𝜆²+𝑚₁𝜆+𝑛₁𝜆+𝑛₀𝑒^−𝜆𝜏2=0. (23) Let 𝜆=𝑖𝜔₂(𝜔₂>0) be a root of (23), then we have 𝑛₁𝜔₂sin𝜏₂𝜔₂+𝑛₀cos𝜏₂𝜔₂=𝑝₂𝜔²₂,𝑛₁𝜔₂cos𝜏₂𝜔₂−𝑛₀sin𝜏₂𝜔₂=𝜔³₂−𝑚₁𝜔₂, (24) which follows that 𝜔⁶₂+𝑐₃₂𝜔⁴₂+𝑐₃₁𝜔²₂+𝑐₃₀=0, (25) with 𝑐₃₂=𝑝²₂−2𝑚₁,𝑐₃₁=𝑚²₁−𝑛²₁,𝑐₃₀=−𝑛²₀. (26) Let omega²₂=𝑣₂, then (24) becomes 𝑣³₂+𝑐₃₂𝑣²₂+𝑐₃₁𝑣₂+𝑐₃₀=0. (27) Similar as in Case 2, we assume that (𝐻₃₁)∶ (27) has at least one positive solution. Without loss of generality, we assume that it has three positive roots, which are denoted by 𝑣₂₁, 𝑣₂₂ and 𝑣₂₃. Then (25) has three positive roots 𝜔_2𝑘=√𝑣_2𝑘, 𝑘=1,2,3.
From (24), we get 𝜏^(𝑗)_2𝑘=1𝜔_2𝑘𝑛arccos¹𝜔⁴_2𝑘+𝑛₀𝑝₂−𝑚₁𝑛₁𝜔²_2𝑘𝑛²₁𝜔²_2𝑘+𝑛²₀+2𝑗𝜋𝜔_2𝑘,𝑗=0,1,2,…,𝑘=1,2,3. (28) Then, we denote 𝜏₂₀=min𝜏⁽⁰⁾_2𝑘,𝑘=1,2,3,𝜔₂₀=𝜔_2𝑘0. (29) Similar as in Case 2, we know that if condition (𝐻₃₂)∶𝑓^′₂(𝑣_2∗)≠0 holds, where 𝑓₂𝑣₂=𝑣³₂+𝑐₃₂𝑣²₂+𝑐₃₁𝑣₂+𝑐₃₀,𝑣_2∗=𝜔²₂₀, (30) then, Re[𝑑𝜆/𝑑𝜏₂]⁻¹_{𝜆=𝑖𝜔20}≠0. Namely, if condition (𝐻₃₂) holds, the transversality condition is satisfied. Therefore, we have the following results. Therefore, we have the following theorem.

Theorem 2. Suppose that conditions (𝐻₃₁)-(𝐻₃₂) hold. The positive equilibrium 𝐸_∗(𝑥_∗,𝑦_∗,𝑧_∗) of system (2) is asymptotically stable for 𝜏₂∈[0,𝜏₂₀) and unstable when 𝜏₂>𝜏₂₀. And system (2) undergoes a Hopf bifurcation at 𝐸_∗(𝑥_∗,𝑦_∗,𝑧_∗) when 𝜏₂=𝜏₂₀.

Case 4. One has 𝜏₁>0,𝜏₂∈[0,𝜏₂₀).
It is considered that with (8), 𝜏₂ in its stable interval and 𝜏₁ is considered as a parameter.
Let 𝜆=𝑖𝜔(𝜔>0) be the root of (8). Separating real and imaginary parts leads to 𝑚₁𝜔sin𝜏₁𝜔=𝑝₂𝜔²−𝑛₀cos𝜏₂𝜔−𝑛₁𝜔sin𝜏₂𝑚𝜔,₁𝜔cos𝜏₁𝜔=𝜔³+𝑛₀sin𝜏₂𝜔−𝑛₁𝜔cos𝜏₂𝜔. (31) Eliminating 𝜏₁ leads to 𝑐₄₀(𝜔)+𝑐₄₁cos𝜏₂𝜔+𝑐₄₂sin𝜏₂𝜔=0, (32) where 𝑐₄₀(𝜔)=𝜔⁶+𝑝²₂𝜔⁴+𝑛²₁+𝑚²₁𝜔²+𝑛²₀,𝑐₄₁(𝜔)=−2𝑛₁𝜔⁴+𝑛₀𝑝₂𝜔²,𝑐₄₂(𝜔)=2𝑛₀−𝑛₁𝑝₂𝜔³. (33) Suppose that (𝐻₄₁)∶ (32) has finite positive roots. If condition (𝐻₄₁) holds, we denote the roots of (32) by 𝜔₁,𝜔₂,…,𝜔_𝑛. For every fixed 𝜔_𝑖(𝑖=1,2,…,𝑛), there exists a sequence {𝜏^(𝑗)_1𝑘∣𝑘=1,2,…,𝑛,𝑗=0,1,2,…} satisfying (32).

From (31), we can get 𝜏^(𝑗)_1𝑘=1𝜔_𝑘𝜔arccos^3𝑘+𝑛₀sin𝜏₂𝜔_𝑘−𝑛₁𝜔_𝑘cos𝜏₂𝜔_𝑘𝑚₁𝜔_𝑘+2𝑗𝜋𝜔_𝑘,𝑘=1,2,…,𝑛,𝑗=0,1,2,…. (34) Let 𝜏^∗₁=𝜏⁽⁰⁾_1𝑘=min𝜏⁽⁰⁾_1𝑘,𝑘=1,2,…,𝑛, (35) when 𝜏₁=𝜏^∗₁ (8) has a pair of purely imaginary roots ±𝑖𝜔^∗₁ for 𝜏₂∈[0,𝜏₂₀).

To verify the transversality condition of Hopf bifurcation, differentiating (8) with respect to 𝜏₁ and substituting 𝜏₁=𝜏^∗₁, we can get Re𝑑𝜆𝑑𝜏₁⁻¹_{𝜆=𝑖𝜔∗1}=𝑃_𝑅𝑄_𝑅+𝑃_𝐼𝑄_𝐼𝑄²_𝑅+𝑄²_𝐼, (36) where 𝑃_𝑅=𝑛₁−𝑛₀𝜏₂cos𝜏₂𝜔^∗₁−𝑛₁𝜔^∗₁𝜏₂sin𝜏₂𝜔^∗₁−3𝜔^∗₁²+𝑚₁cos𝜏^∗₁𝜔^∗₁,𝑃_𝐼=𝑛₁−𝑛₀𝜏₂sin𝜏₂𝜔^∗₁−𝑛₁𝜔^∗₁𝜏₂cos𝜏₂𝜔^∗₁+2𝑝₂𝜔^∗₁−𝑚₁sin𝜏^∗₁𝜔^∗₁,𝑄_𝑅=−𝑚₁𝜔^∗₁²cos𝜏^∗₁𝜔^∗₁,𝑄_𝐼=𝑚₁𝜔^∗₁²sin𝜏^∗₁𝜔^∗₁. (37) Clearly, if condition (𝐻₄₂)∶𝑃_𝑅𝑄_𝑅+𝑃_𝐼𝑄_𝐼≠0 holds, then [Re(𝑑𝜆/𝑑𝜏)]⁻¹_{𝜆=𝑖𝜔}^∗₁≠0. Namely, if condition (𝐻₄₂) holds, the transversality condition is satisfied. Therefore, we have the following results. Thus, we have the following theorem.

Theorem 3. Suppose that conditions (𝐻₄₁)-(𝐻₄₂) hold and 𝜏₂∈[0,𝜏₂₀). The positive equilibrium 𝐸_∗(𝑥_∗,𝑦_∗,𝑧_∗) of system (2) is asymptotically stable for 𝜏₁∈[0,𝜏_1∗) and unstable when 𝜏₁>𝜏_1∗. And system (2) undergoes a Hopf bifurcation at 𝐸_∗(𝑥_∗,𝑦_∗,𝑧_∗) when 𝜏₁=𝜏_1∗.

Case 5. One has 𝜏₂>0,𝜏₁∈[0,𝜏₁₀).
We consider (8) with 𝜏₁ in its stable interval, regarding 𝜏₂ as a parameter.
Let 𝜆=𝑖𝜔(𝜔>0) be a root of (8). Then we get 𝑛₁𝜔sin𝜏₂𝜔+𝑛₀cos𝜏₂𝜔=𝑝₂𝜔²−𝑚₁𝜔sin𝜏₁𝑛𝜔,₁𝜔cos𝜏₂𝜔−𝑛₀sin𝜏₂𝜔=𝜔³−𝑚₁𝜔cos𝜏₁𝜔. (38) It follows that 𝑐₅₀(𝜔)+𝑐₅₁(𝜔)cos𝜏₁𝜔+𝑐₅₂(𝜔)sin𝜏₁𝜔=0, (39) where 𝑐₅₀(𝜔)=𝜔⁶+𝑝²₂𝜔⁴+𝑚²₁−𝑛²₁𝜔²−𝑛²₀,𝑐₅₁(𝜔)=−2𝑚₁𝜔⁴,𝑐₅₂(𝜔)=−2𝑚₁𝑝₂𝜔³. (40) Similar as in Case 4, we suppose that (𝐻₅₁)∶ (39) has finite positive roots. And we denote the roots of (39) by 𝜔₁,𝜔₂,…,𝜔_𝑛. The corresponding critical value of 𝜏₂ is 𝜏^(𝑗)_2𝑘=𝑛₁𝜔²_𝑘+𝑛₀𝑝₂𝜔²_𝑘−𝑚₁𝑛₁𝜔²_𝑘cos𝜏₁𝜔_𝑘−𝑚₁𝑛₀𝜔_𝑘sin𝜏₁𝜔_𝑘𝑛₂𝜔²_𝑘+𝑛²₀,𝑘=1,2,…,𝑛,𝑗=0,1,2,…. (41) Let 𝜏^∗₂=𝜏⁽⁰⁾_2𝑘=min𝜏⁽⁰⁾_2𝑘,𝑘=1,2,…,𝑛, (42) when 𝜏₂=𝜏^∗₂ (8) has a pair of purely imaginary roots ±𝑖𝜔^∗₂ for 𝜏₁∈[0,𝜏₁₀).

Similar as in Case 4, we give the following assumption (𝐻₅₂)∶𝑃^′_𝑅𝑄^′_𝑅+𝑃^′_𝐼𝑄^′_𝐼≠0, where 𝑃^′_𝑅=𝑚₁cos𝜏₁𝜔^∗₂−𝜏₁𝑚₁𝜔^∗₂sin𝜏₁𝜔^∗₂−3𝜔^∗₂²+𝑛₁cos𝜏^∗₂𝜔^∗₂,𝑃^′_𝐼=−𝑚₁sin𝜏₁𝜔^∗₂−𝜏₁𝑚₁𝜔^∗₂cos𝜏₁𝜔^∗₂+2𝑝₂𝜔^∗₂−𝑛₁sin𝜏^∗₂𝜔^∗₂,𝑄^′_𝑅=𝑛₀𝜔^∗₂sin𝜏^∗₂𝜔^∗₂−𝑛₁𝜔^∗₂²cos𝜏^∗₂𝜔^∗₂,𝑄^′_𝐼=𝑛₀𝜔^∗₂cos𝜏^∗₂𝜔^∗₂+𝑛₁𝜔^∗₂²sin𝜏^∗₂𝜔^∗₂. (43)

Therefore, if condition (𝐻₅₂)∶𝑃^′_𝑅𝑄^′_𝑅+𝑃^′_𝐼𝑄^′_𝐼≠0 holds, then we can get Re[𝑑𝜆/𝑑𝜏₂]⁻¹_{𝜆=𝑖𝜔}^∗₂. That is, the transversality condition is satisfied. Hence, we have the following theorem.

Theorem 4. Suppose that conditions (𝐻₅₁)-(𝐻₅₂) hold and 𝜏₁∈[0,𝜏₁₀). The positive equilibrium 𝐸_∗(𝑥_∗,𝑦_∗,𝑧_∗) of system (2) is asymptotically stable for 𝜏₂∈[0,𝜏^∗₂) and unstable when 𝜏₂>𝜏^∗₂. And system (2) undergoes a Hopf bifurcation at 𝐸_∗(𝑥_∗,𝑦_∗,𝑧_∗) when 𝜏₂=𝜏^∗₂.

3. Direction and Stability of the Hopf Bifurcation

In this section, we will employ the normal form method and center manifold theorem introduced by Hassard et al. [17] to determine the direction of Hopf bifurcation and stability of the bifurcated periodic solutions of system (2) with respect to 𝜏₁ for 𝜏₂∈(0,𝜏₂₀). Without loss of generality, we assume that 𝜏^′₂<𝜏^∗₁, where 𝜏^′₂∈(0,𝜏₂₀).

Let 𝜏₁=𝜏^∗₁+𝜇,𝜇∈𝑅. Then 𝜇=0 is the Hopf bifurcation value of system (2). Rescaling the time delay 𝑡→(𝑡/𝜏₁), then system (2) can be rewritten as 𝑢̇(𝑡)=𝐿_𝜇𝑢_𝑡+𝐹𝜇,𝑢_𝑡, (44) where 𝐿_𝜇𝜙=𝜏_1∗(0)−𝜏+𝜇𝐴𝜙+𝐶𝜙^′2𝜏^∗₁(),𝐹=𝜏+𝐵𝜙−1𝜇,𝜙^∗₁+𝜇𝐹₁,𝐹₂,𝐹₃^𝑇 (45) With 𝑎𝐴=¹¹𝑎₁₂000𝑎₂₃𝑏000,𝐵=000₁₁00,000𝐶=0000000𝑐₃₂0,𝐹₁=𝑎₁₃𝜙²₁(0)+𝑎₁₄𝜙₁(0)𝜙₂(0)+𝑎₁₅𝜙³₁(0)𝜙₁(0)+𝑎₁₈𝜙²₁(0)𝜙₂(0)𝐹+⋯,₂=𝑎₂₄𝜙₂(0)𝜙₃(0)+𝑎₂₅𝜙₁(−1)𝜙₂(0)+𝑎₂₆𝜙²₁(−1)+𝑎₂₇𝜙²₁(−1)𝜙₂(0)+𝑎₂₈𝜙³₁()𝐹−1⋯,₃=𝑐₃₃𝜙₂−𝜏^′₁𝜏^∗₁𝜙₃(0). (46) By Riesz representation theorem, there exists a 3×3 matrix function 𝜂(𝜃,𝜇)∶[−1,0]→𝑅³ whose elements are of bounded variation, such that 𝐿_𝜇𝜙=⁰₋₁𝑑𝜂𝜃,𝜇𝜙(𝜃),𝜙∈𝐶[−1,0],𝑅³. (47) In fact, we can choose 𝜂=𝜏𝜃,𝜇^∗1()𝜏+𝜇𝐴+𝐵+𝐶,𝜃=0,_∗1()−𝜏+𝜇𝐵+𝐶,𝜃∈^′2𝜏^∗₁,𝜏,0_∗1𝜏+𝜇𝐶,𝜃∈−1,−^′2𝜏^∗₁,0,𝜃=−1. (48) For 𝜙∈𝐶([−1,0],𝑅³), we define 𝐴𝜇(𝜃)𝜙=𝑑𝜙∫𝑑𝜃,−1≤𝜃<0,₀₋₁𝜙(𝜃)𝑅𝜇𝐹𝑑𝜂𝜃,𝜇,𝜃=0,𝜙=0,−1≤𝜃<0,𝜇,𝜙,𝜃=0. (49) Then system (44) can be transformed into the following operator equation: 𝑢̇(𝑡)=𝐴𝜇𝑢_𝑡+𝑅𝜇𝑢_𝑡. (50)

The adjoint operator 𝐴^∗ of 𝐴 is defined by 𝐴^∗𝜇−(𝑠)𝜑=𝑑𝜑∫𝑑𝑠,0<𝑠≤1,₀₋₁𝑑𝜂^𝑇𝑠,𝜇𝜑(−𝑠),𝑠=0, (51) associated with a bilinear form ⟨𝜑,𝜙⟩=𝜑(0)𝜙(0)−⁰_𝜃=−1^𝜃_𝜉=0𝜑(𝜉−𝜃)𝑑𝜂(𝜃)𝜙(𝜉)𝑑𝜉, (52) where 𝜂(𝜃)=𝜂(𝜃,0).

From the above discussion, we know that ±𝑖𝜏^∗₁𝜔^∗₁ are the eigenvalues of 𝐴(0) and they are also eigenvalues of 𝐴^∗(0). We assume that 𝜌(𝜃)=(1,𝜌₂,𝜌₃)^𝑇𝑒^{𝑖𝜏∗1𝜔∗1𝜃} is the eigenvector belonging to the eigenvalue +𝑖𝜏^∗₁𝜔^∗₁ and 𝜌^∗(𝜃)=𝐷(1,𝜌^∗₂,𝜌^∗₃)𝑒^{𝑖𝜏∗1𝜔∗1𝑠} is the eigenvector belonging to the eigenvalue −𝑖𝜏^∗₁𝜔^∗₁. Then, by a simple computation, we can obtain 𝜌₂=𝑖𝜔^∗₁−𝑎₁₁𝑎₁₂,𝜌₃=𝑖𝜔^∗₁𝑖𝜔^∗₁−𝑎₁₁−𝑏₁₁𝑒^{−𝑖𝜔∗1𝜔∗1}𝑎₁₂𝑎₂₃,𝜌^∗₂=−𝑖𝜔^∗₁+𝑎₁₁𝑏₁₁𝑒^{𝑖𝜏∗1𝜔∗1},𝜌^∗₃=𝑖𝜔^∗₁𝑖𝜔^∗₁+𝑎₁₁𝑒^{−𝑖𝜏∗1𝜔∗1}−𝑎₁₂𝑏₁₁𝑏₁₁𝑐₃₂𝑒^{𝑖𝜏′2𝜔∗1},𝐷=1+𝜌₂𝜌^∗₂+𝜌₃𝜌^∗₃+𝑏₁₁𝜌^∗₂𝜏^∗₁𝑒^{−𝑖𝜏∗1𝜔∗1}+𝑐₃₂𝜌₂𝜌^∗₃𝜏^′₂𝑒^{−𝑖𝜏′2𝜔∗1−1}. (53) Then we have ⟨𝜌^∗,𝜌⟩=1.

Next, we get the coefficients used to determine the important quantities of the periodic solution by using a computation process similar to that in [18]: 𝑔₂₀=2𝜏^∗₁𝐷𝑎₁₃+𝑎₁₄𝜌⁽²⁾(0)+𝜌^∗₂𝑎₂₄𝜌⁽²⁾(0)𝜌⁽³⁾(0)+𝑎₂₅𝜌⁽¹⁾(−1)𝜌⁽²⁾(0)+𝑎₂₆𝜌⁽¹⁾(−1)²+𝜌^∗₃𝑐₃₃𝜌⁽²⁾−𝜏^′₂𝜏^∗₁𝜌⁽³⁾(0,𝑔₁₁)=𝜏^∗₁𝐷2𝑎₁₃+𝑎₁₄𝜌⁽²⁾(0)+𝜌⁽²⁾(0)+𝜌^∗₂𝑎₂₄𝜌⁽²⁾(0)𝜌⁽³⁾(0)+𝜌⁽²⁾(0)𝜌⁽³⁾(0)+𝑎₂₅𝜌⁽¹⁾(−1)𝜌⁽²⁾(0)+𝜌⁽¹⁾(−1)𝜌⁽²⁾(0)+2𝑎₂₆𝜌⁽¹⁾(−1)𝜌⁽¹⁾(−1)+𝑐₃₃𝜌^∗₃𝜌⁽²⁾−𝜏^′₂𝜏^∗₁𝜌⁽³⁾(0)+𝜌⁽²⁾−𝜏^′₂𝜏^∗₁𝜌⁽³⁾(0),𝑔₀₂=2𝜏^∗₁𝐷𝑎₁₃+𝑎₁₄𝜌⁽²⁾(0)+𝜌^∗₂𝑎₂₄𝜌⁽²⁾(0)𝜌⁽³⁾(0)+𝑎₂₅𝜌⁽¹⁾(−1)𝜌⁽²⁾(0)+𝑎₂₆𝜌⁽¹⁾(−1)²+𝜌^∗₃𝑐₃₃𝜌⁽²⁾−𝜏^′₂𝜏^∗₁𝜌⁽³⁾(0,𝑔₂₁)=2𝜏^∗₁𝐷𝑎₁₃2𝑊⁽¹⁾₁₁(0)+𝑊⁽¹⁾₂₀(0)+𝑎₁₄𝑊⁽¹⁾₁₁(0)𝜌⁽²⁾(0)+12𝑊⁽¹⁾₂₀(0)𝜌⁽²⁾(0)+𝑊⁽²⁾₁₁(0)+12𝑊⁽²⁾₂₀(0)+3𝑎₁₅+𝑎₁₆𝜌⁽²⁾(0)+2𝜌⁽²⁾(0+)𝜌^∗₂𝑎₂₄𝑊⁽²⁾₁₁(0)𝜌⁽³⁾(0)+12𝑊⁽²⁾₂₀(0)𝜌⁽³⁾(0)+𝑊⁽³⁾₁₁(0)𝜌⁽²⁾(0)+12𝑊⁽³⁾₂₀(0)𝜌⁽²⁾(0)+𝑎₂₅𝑊⁽¹⁾₁₁(−1)𝜌⁽²⁾(0)+12𝑊⁽¹⁾₂₀(−1)𝜌⁽²⁾(0)+𝑊⁽²⁾₁₁(0)𝜌⁽²⁾()+1−12𝑊⁽²⁾₂₀(0)𝜌⁽¹⁾(−1)+𝑎₂₆2𝑊⁽¹⁾₁₁(−1)𝜌⁽¹⁾(−1)+𝑊⁽¹⁾₂₀(−1)𝜌⁽¹⁾(−1)+𝑎₂₇𝜌⁽¹⁾(−1)²+2𝜌⁽¹⁾×(−1)𝜌⁽²⁾(0)𝜌⁽¹⁾(−1)+3𝑎₂₈𝜌⁽¹⁾(−1)²𝜌⁽¹⁾(−1)+𝑐₃₃𝜌^∗₃𝑊⁽²⁾₁₁−𝜏^′₂𝜏^∗₁𝜌⁽³⁾(0)+12𝑊⁽²⁾₂₀−𝜏^′₂𝜏^∗₁𝜌⁽³⁾(0)+𝑊⁽³⁾₁₁(0)𝜌⁽²⁾−𝜏^′₂𝜏^∗₁+12𝑊⁽³⁾₂₀(0)𝜌⁽²⁾−𝜏^′₂𝜏^∗₁, (54) with 𝑊₂₀(𝜃)=𝑖𝑔₂₀𝑞(0)𝜏^∗₁𝜔^∗₁𝑒^{𝑖𝜏∗1𝜔∗1𝜃}+𝑖𝑔₀₂𝑞(0)3𝜏^∗₁𝜔^∗₁𝑒^{−𝑖𝜏∗1𝜔∗1𝜃}+𝐸₁𝑒^{2𝑖𝜏∗1𝜔∗1𝜃},𝑊₁₁(𝜃)=−𝑖𝑔₁₁𝑞(0)𝜏^∗₁𝜔^∗₁𝑒^{𝑖𝜏∗1𝜔∗1𝜃}+𝑖𝑔₁₁𝑞(0)𝜏^∗₁𝜔^∗₁𝑒^{−𝑖𝜏∗1𝜔∗1𝜃}+𝐸₂, (55) where 𝐸₁ and 𝐸₂ can be computed as the following equations, respectively 2𝑖𝜔^∗₁−𝑎₁₁−𝑎₁₂0−𝑏₁₁𝑒^{−2𝑖𝜏∗1𝜔∗1}2𝑖𝜔^∗₁−𝑎₂₃0−𝑐₃₂𝑒^{−2𝑖𝜏∗1𝜔∗1}2𝑖𝜔^∗₁𝐸₁𝐸=2𝐸¹⁽¹⁾₁⁽²⁾𝐸₁⁽³⁾,𝑎₁₁𝑎₁₂0𝑏₁₁𝑎₂₂00𝑐₃₂0𝐸₂𝐸=−𝐸²⁽¹⁾₂⁽²⁾𝐸₂⁽³⁾, (56) with 𝐸₁⁽¹⁾=𝑎₁₃+𝑎₁₄𝜌⁽²⁾(0),𝐸₁⁽²⁾=𝑎₂₄𝜌⁽²⁾(0)𝜌⁽³⁾(0)+𝑎₂₅𝜌⁽¹⁾(−1)𝜌⁽²⁾(0)+𝑎₂₆𝜌⁽¹⁾(−1))²,𝐸₁⁽³⁾=𝑐₃₃𝜌⁽²⁾−𝜏^′₂𝜏^∗₁𝜌⁽³⁾(0),𝐸₂⁽¹⁾=2𝑎₁₃+𝑎₁₄𝜌⁽²⁾(0)+𝜌⁽²⁾(0),𝐸₂⁽²⁾=𝑎₂₄𝜌⁽²⁾(0)𝜌⁽³⁾(0)+𝜌⁽²⁾(0)𝜌⁽³⁾(0)+𝑎₂₅𝜌⁽¹⁾(−1)𝜌⁽²⁾(0)+𝜌⁽¹⁾(−1)𝜌⁽²⁾(0)+2𝑎₂₆𝜌⁽¹⁾(−1)𝜌⁽¹⁾(),𝐸−1₂₍₃₎=𝑐₃₃𝜌⁽²⁾−𝜏^′₂𝜏^∗₁𝜌⁽³⁾(0)+𝜌⁽²⁾−𝜏^′₂𝜏^∗₁𝜌⁽³⁾(0). (57) Thus, we can calculate the following values: 𝐶₁(0)=𝑖2𝜏^∗₁𝜔^∗₁𝑔₁₁𝑔₂₀−2𝑔₁₁²−𝑔₀₂²3+𝑔₂₁2,𝜇₂𝐶=−Re₁(0)Re𝜆^′𝜏^∗₁,𝛽₂=2Re𝐶₁(0,𝑇₂)𝐶=−Im₁(0)+𝜇₂Im𝜆^′𝜏^∗₁𝜏^∗₁𝜔^∗₁. (58) Based on the discussion above, we can obtain the following results.

Theorem 5. If 𝜇₂>0(𝜇₂<0), then the Hopf bifurcation is supercritical (subcritical); if 𝛽₂<0(𝛽₂>0), the bifurcating periodic solutions are stable (unstable); if 𝑇₂>0(𝑇₂<0), the period of the bifurcating periodic solutions increases (decreases).

4. Numerical Simulation and Discussion

In this section, we present some numerical simulations to illustrate the analytical results obtained in the previous sections. Let 𝛼=3, 𝑘=2, 𝛽=2, 𝑝=0.3, ℎ=0.2, 𝑒=0.4, 𝑟=1, 𝑠=0.5, and 𝑚=0.6. Then we have the following particular case of system (2): 𝑑𝑥(𝑡)(𝑡)𝑥(𝑡)𝑑𝑡=3𝑥1−2−2𝑥(𝑡)𝑦(𝑡)(𝑡),(𝑡)1+0.3𝑥𝑑𝑦(𝑡)𝑑𝑡=𝑦−0.2+0.8𝑥𝑡−𝜏₁1+0.3𝑥𝑡−𝜏₁(𝑡)𝑧(𝑡),(𝑡)−𝑦𝑑𝑧(𝑡)𝑑𝑡=𝑧−0.5+0.6𝑦𝑡−𝜏₂, (59) which has a positive equilibrium 𝐸_∗(1.1809,0.8333,0.4976).

For 𝜏₁>0,𝜏₂=0, we have 𝜔₁₀=0.7318, 𝜏₁₀=1.4887. From Theorem 1, we know that the positive equilibrium 𝐸_∗(1.1809,0.8333,0.4976) is asymptotically stable for 𝜏₁∈[0,1.4887). As can be seen from Figure 1, if 𝜏₁=1.38∈[0,1.4887), 𝐸_∗(1.1809,0.8333,0.4976) is asymptotically stable. However, if 𝜏₁=1.498>𝜏₁₀=1.4887, then 𝐸_∗(1.1809,0.8333,0.4976) is unstable and system (59) undergoes a Hopf bifurcation at 𝐸_∗(1.1809,0.8333,0.4976), and a family of periodic solutions bifurcate from the positive equilibrium 𝐸_∗(1.1809,0.8333,0.4976). This property can be illustrated by Figure 2. For 𝜏₂>0,𝜏₁=0, by a simple computation, we can easily get 𝜔₂₀=0.4737, 𝜏₂₀=1.7990. The corresponding waveform and the phase plots are shown in Figures 3 and 4.

Figure 1 

𝐸∗ is asymptotically stable for 𝜏1=1.38<𝜏10=1.4887 with initial values 0.975, 1.025, and 0.625.

Figure 2 

𝐸∗ is unstable for 𝜏1=1.498>𝜏10=1.4887 with initial values 0.975, 1.025, and 0.625.

Figure 3 

𝐸∗ is asymptotically stable for 𝜏2=1.60<𝜏20=1.7990 with initial values 0.975, 1.025, and 0.625.

Figure 4 

𝐸∗ is unstable for 𝜏2=1.96>𝜏20=1.7990 with initial values 0.975, 1.025, and 0.625.

For 𝜏₁>0 and 𝜏₂=0.8∈[0,𝜏₂₀), we get 𝜔^∗₁=0.6698, 𝜏^∗₁=0.8095. That is, when 𝜏₁ increases from zero to the critical value 𝜏^∗₁, the positive equilibrium 𝐸_∗(1.1809,0.8333,0.4976) is asymptotically stable; then it will lose stability, and a Hopf bifurcation occurs once 𝜏₁>𝜏^∗₁=0.8095. This property can be illustrated by Figures 5 and 6. Further, we get 𝐶₁(0)=−11.2213+16.3520𝑖, [𝑑𝜆/𝑑𝜏₁]_{𝜆=𝑖𝜔}^∗₁=4.1212−2.5252𝑖. Then we have 𝜇₂=2.7228, 𝛽₂=−22.4426, 𝑇₂=−17.4776. Therefore, from Theorem 5, we can know that the Hopf bifurcation is supercritical and the bifurcating periodic solutions are stable.

Figure 5 

𝐸∗ is asymptotically stable for 𝜏1=0.75<𝜏∗1=0.8095 and 𝜏2=0.8 with initial values 0.975, 1.025, and 0.625.

Figure 6 

𝐸∗ is unstable for 𝜏1=0.936>𝜏∗1=0.8095 and 𝜏2=0.8 with initial values 0.975, 1.025, and 0.625.

At last, for 𝜏₂>0 and 𝜏₁=0.5∈(0,𝜏₁₀), we obtain 𝜔^∗₂=0.5525, 𝜏^∗₂=1.4818. The corresponding waveform and the phase plots are shown in Figures 7 and 8.

Figure 7 

𝐸∗ is asymptotically stable for 𝜏2=1.25<𝜏∗2=1.4818 and 𝜏1=0.5 with initial values 0.975, 1.025, and 0.625.

Figure 8 

𝐸∗ is unstable for 𝜏2=1.50>𝜏∗2=1.4818 and 𝜏1=0.5 with initial values 0.975, 1.025, and 0.625.

Guo and Jiang [7] have obtained that the three species in system (2) with only one time delay can coexist, however, we get that the species could also coexist with some available time delays of the predator and the top predator. This is valuable from the view of ecology. As the future work, we shall consider the following more general and more complicated system with multiple delays: 𝑑𝑥(𝑡)(𝑡)𝑥𝑑𝑡=𝛼𝑥1−𝑡−𝜏₁𝑘−𝛽𝑥(𝑡)𝑦(𝑡)(𝑡),(𝑡)1+𝑝𝑥𝑑𝑦(𝑡)𝑑𝑡=𝑦−ℎ+𝑒𝛽𝑥𝑡−𝜏₂1+𝑝𝑥𝑡−1𝜏₂(𝑡)𝑧(𝑡),(𝑡)−𝑟𝑦𝑑𝑧(𝑡)𝑑𝑡=𝑧−𝑠+𝑚𝑟𝑦𝑡−𝜏₃, (60) where 𝜏₁ is feedback delay of the prey and 𝜏₂, 𝜏₃ are the time delays due to the gestation of the predator and the top predator, respectively.

Acknowledgments

This work was supported by National Natural Science Foundation of China (11072090), Natural Science Foundation of the Higher Education Institutions of Anhui Province (KJ2013B137) and Anhui Provincial Natural Science Foundation under Grant no. 1208085QA11.