Stability and Hopf Bifurcation in a Modified Holling-Tanner Predator-Prey System with Multiple Delays

Abstract

A modified Holling-Tanner predator-prey system with multiple delays is investigated. By analyzing the associated characteristic equation, the local stability and the existence of periodic solutions via Hopf bifurcation with respect to both delays are established. Direction and stability of the periodic solutions are obtained by using normal form and center manifold theory. Finally, numerical simulations are carried out to substantiate the analytical results.

1. Introduction

Predator-prey dynamics has long been and will continue to be of interest to both applied mathematicians and ecologists due to its universal existence and importance [1, 2]. Many population models investigating the dynamic relationship between predators and their preys have been proposed and studied. For example, Lotka-Volterra model [3–5], Leslie-Gower model [6–10], and Holling-Tanner model [11–16]. Among these widely used models, Holling-Tanner model plays a special role in view of the interesting dynamics it possesses. Holling-Tanner model for predator-prey interaction is governed by the following nonlinear coupled ordinary differential equations: 𝑑𝑋𝑋𝑑𝑇=𝑟𝑋1−𝐾−𝑚𝑋𝑌,𝑎+𝑋𝑑𝑌𝑠𝑌𝑑𝑇=𝑌1−ℎ𝑋, (1.1) where 𝑋 and 𝑌 denote the population densities of prey species and predator species at time 𝑇, respectively. The first equation in system (1.1) shows that the prey grows logistically with the carrying capacity 𝐾 and the intrinsic growth rate 𝑟 in the absence of the predator. And the growth of the prey is hampered by the predator at a rate proportional to the functional response 𝑚𝑋/(𝑎+𝑋) in the presence of the predator. The second equation shows that the predator consumes the prey according to the functional response 𝑚𝑋/(𝑎+𝑋) and grows logistically with the intrinsic growth rate 𝑠 and carrying capacity 𝑋/ℎ proportional to the number of the prey. The parameter 𝑚 denotes the maximal predator per capita consumption rate. 𝑎 is a saturation value; it corresponds to the number of prey necessary to achieve one half the maximum rate 𝑚. The parameter ℎ denotes the number of prey required to support one predator at equilibrium when 𝑦 equals 𝑋/ℎ.

Recently, there has been considerable interest in predator-prey systems with the Beddington-DeAngelis functional response. And it has been shown that the predator-prey systems with the Beddington-DeAngelis functional response have rich but biologically reasonable dynamics. For more details about this functional response one can refer to [17–21]. Zhang [16], Lu and Liu [22] considered the following modified Holling-Tanner delayed predator-prey system: 𝑑𝑋𝑋𝑑𝑇=𝑟𝑋1−𝐾−𝛼𝑋𝑌,𝑎+𝑏𝑋+𝑐𝑌𝑑𝑌𝑠𝑌(𝑑𝑇=𝑌1−ℎ𝑡−𝜏)𝑋(𝑡−𝜏), (1.2) where 𝜏 is incorporated in the negative feedback of the predator density. 𝛼𝑋𝑌/(𝑎+𝑏𝑋+𝑐𝑌) is the Beddington-DeAngelis functional response. The parameters 𝛼, 𝑎, 𝑏, and 𝑐 are assumed to be positive. 𝛼 is the maximum value at which per capita reduction rate of the prey can attain. 𝑎 measures the extent to which environment provides protection to the prey. 𝑏 describes the effect of handling time on the feeding rate. 𝑐 describes the magnitude of interference among predators. Zhang [16] investigated the local Hopf bifurcation of system (1.2). Lu and Liu [22] proved that system (1.2) is permanent under some conditions and obtained the sufficient conditions of local and global stability of system (1.2). Since both the species are growing logistically, it is reasonable to assume delay in prey species as well. Based on this consideration, we incorporate the negative feedback of the prey density into system (1.2) and obtain the following system: 𝑑𝑋𝑋𝑑𝑇=𝑟𝑋1−𝑇−𝑇₁𝐾−𝛼𝑋𝑌,𝑎+𝑏𝑋+𝑐𝑌𝑑𝑌𝑠𝑌𝑑𝑇=𝑌1−ℎ𝑇−𝑇₂𝑋𝑇−𝑇₂, (1.3) where 𝑇₁ and 𝑇₂ are the feedback time delays of the prey density and the predator density respectively. Let 𝑋=𝐾𝑥, 𝑌=(𝑟𝐾/𝛼)𝑦, 𝑡=𝑟𝑇, 𝜏₁=𝑟𝑇₁, 𝜏₂=𝑟𝑇₂, system (1.3) can be transformed into the following nondimensional form: 𝑑𝑥𝑑𝑡=𝑥1−𝑥𝑡−𝜏₁−𝑥𝑦𝑎₁+𝑏𝑥+𝑐₁𝑦,𝑑𝑦𝑦𝑑𝑡=𝑦𝛿−𝛽𝑡−𝜏₂𝑥𝑡−𝜏₂, (1.4) where 𝑎₁=𝑎/𝐾,𝑐₁=𝑐𝑟/𝛼,𝛿=𝑠/𝑟,𝛽=𝑠ℎ/𝛼 are the non-dimensional parameters and they are positive.

The main purpose of this paper is to consider the effect of multiple delays on system (1.4). The local stability of the positive equilibrium and the existence of Hopf bifurcation are investigated. By employing normal form and center manifold theory, the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions are determined. Finally, some numerical simulations are also included to illustrate the theoretical analysis.

2. Local Stability and the Existence of Hopf Bifurcation

In this section, we study the local stability of each of feasible equilibria and the existence of Hopf bifurcation at the positive equilibrium. Obviously, system (1.4) has a unique boundary equilibrium 𝐸₁(1,0) and a unique positive equilibrium 𝐸_∗(𝑥_∗,𝑦_∗), where 𝑥_∗=−𝑎₁−𝑏𝛽+1−𝑐₁𝛿+𝑎₁−𝑏𝛽+1−𝑐₁𝛿²+4𝑎₁𝛽𝑏𝛽+𝑐₁𝛿2𝑏𝛽+𝑐₁𝛿,𝑦_∗=𝛿𝛽𝑥_∗. (2.1) The Jacobian matrix of system (1.4) at 𝐸₁ takes the form 𝐽𝐸₁=−𝑒^−𝜆𝜏1−1𝑎₁+𝑏0𝛿. (2.2) The characteristic equation of system (1.4) at 𝐸₁ is of the form (𝜆−𝛿)𝜆+𝑒^−𝜆𝜏1=0. (2.3) Clearly, the boundary equilibrium 𝐸₁(1,0) is unable.

Next, we discuss the existence of Hopf bifurcation at the positive equilibrium 𝐸(𝑥_∗,𝑦_∗). Let 𝑥(𝑡)=𝑧₁(𝑡)+𝑥_∗, 𝑦(𝑡)=𝑧₂(𝑡)+𝑦_∗, and still denote 𝑧₁(𝑡) and 𝑧₂(𝑡) by 𝑥(𝑡) and 𝑦(𝑡), respectively, then system (1.4) becomes 𝑑𝑥𝑑𝑡=𝑎₁₁𝑥(𝑡)+𝑎₁₂𝑦(𝑡)+𝑏₁₁𝑥𝑡−𝜏₁+_{𝑖+𝑗+𝑘≥2}𝑓₁^𝑖𝑗𝑘𝑥^𝑖𝑦^𝑗𝑥^𝑘𝑡−𝜏₁,𝑑𝑦𝑑𝑡=𝑐₂₁𝑥𝑡−𝜏₂+𝑐₂₂𝑦𝑡−𝜏₂+_{𝑖+𝑗+𝑘≥2}𝑓₂^𝑖𝑗𝑘𝑦^𝑖𝑥^𝑗𝑡−𝜏₂𝑦^𝑘𝑡−𝜏₂, (2.4) where 𝑎₁₁=𝑏𝑥_∗𝑦_∗𝑎₁+𝑏𝑥_∗+𝑐₁𝑦_∗²,𝑎₁₂𝑎=−¹+𝑏𝑥_∗𝑥_∗𝑎₁+𝑏𝑥_∗+𝑐₁𝑦_∗²,𝑏₁₁=−𝑥_∗,𝑐₂₁=𝛿²𝛽,𝑐₂₂𝑓=−𝛿,₁=𝑥1−𝑥𝑡−𝜏₁−𝑥𝑦𝑎₁+𝑏𝑥+𝑐₁𝑦,𝑓₂𝑦=𝑦𝛿−𝛽𝑡−𝜏₂𝑥𝑡−𝜏₂,𝑓₁^𝑖𝑗𝑘𝑥^𝑖𝑦^𝑗𝑥^𝑘𝑡−𝜏₁=1𝜕𝑖!𝑗!𝑘!^{𝑖+𝑗+𝑘}𝑓₁𝜕𝑥^𝑖𝜕𝑦^𝑗𝜕𝑥^𝑘𝑡−𝜏₁|_{(𝑥∗,𝑦∗)},𝑓₂^𝑖𝑗𝑘𝑦^𝑖𝑥^𝑗𝑡−𝜏₂𝑦^𝑘𝑡−𝜏₂=1𝜕𝑖!𝑗!𝑘!^{𝑖+𝑗+𝑘}𝑓₂𝜕𝑦^𝑖𝜕𝑥^𝑗𝑡−𝜏₂𝜕𝑦^𝑘𝑡−𝜏₂∣_{(𝑥∗,𝑦∗)}. (2.5) Then we can obtain the linearized system of system (1.4) 𝑑𝑥𝑑𝑡=𝑎₁₁𝑥(𝑡)+𝑎₁₂𝑦(𝑡)+𝑏₁₁𝑥𝑡−𝜏₁,𝑑𝑦𝑑𝑡=𝑐₂₁𝑥𝑡−𝜏₂+𝑐₂₂𝑦𝑡−𝜏₂. (2.6) The characteristic equation of system (2.6) is 𝜆²−𝑎₁₁𝜆−𝑏₁₁𝜆𝑒^−𝜆𝜏1+𝑎₁₁𝑐₂₂−𝑎₁₂𝑐₂₁−𝑐₂₂𝜆𝑒^−𝜆𝜏2+𝑏₁₁𝑐₂₂𝑒^{−𝜆(𝜏1+𝜏2)}=0. (2.7)

Case 1. 𝜏₁=𝜏₂=𝜏=0.
The characteristic equation of system (1.4) becomes 𝜆²+(𝐴+𝐵+𝐷)𝜆+𝐶+𝐸=0, (2.8) where 𝐴=−𝑎₁₁,𝐵=−𝑏₁₁,𝐶=𝑎₁₁𝑐₂₂−𝑎₁₂𝑐₂₁,𝐷=−𝑐₂₂,𝐸=𝑏₁₁𝑐₂₂. (2.9) It is easy to verify that 𝐶+𝐸=𝛿𝑥_∗+𝑎₁𝛿²𝑥_∗𝛽𝑎₁+𝑏𝑥_∗+𝑐₁𝑦_∗²>0. (2.10) Therefore, if (𝐻₁)∶𝐴+𝐵+𝐷>0, the roots of (2.8) must have negative real parts. Then, we know that the positive equilibrium 𝐸_∗(𝑥_∗,𝑦_∗) of system (1.4) is locally stable in the absence of delay, if (𝐻₁) holds.

Case 2. 𝜏₁=𝜏₂=𝜏>0.
The associated characteristic equation of the system is 𝜆²+𝐴₁𝜆+𝐵₁+𝐶₁𝜆𝑒^−𝜆𝜏+𝐷₁𝑒^−2𝜆𝜏=0, (2.11) where 𝐴₁=−𝑎₁₁,𝐵₁=𝑎₁₁𝑐₂₂−𝑎₁₂𝑐₂₁,𝐶₁=−𝑏₁₁+𝑐₂₂,𝐷₁=𝑏₁₁𝑐₂₂. (2.12) Multiplying 𝑒^𝜆𝜏 on both sides of (2.11), we can obtain 𝜆²+𝐴₁𝜆𝑒^𝜆𝜏+𝐵₁+𝐶₁𝜆+𝐷₁𝑒^−𝜆𝜏=0. (2.13) Now, for 𝜏>0, if 𝜆=𝑖𝜔(𝜔>0) be a root of (2.13). Then, we have 𝐷₁−𝜔²cos𝜏𝜔−𝐴₁𝜔sin𝜏𝜔=−𝐵₁,𝐷₁+𝜔²sin𝜏𝜔−𝐴₁𝜔cos𝜏𝜔=𝐶₁𝜔. (2.14) It follows from (2.14) that 𝐶sin𝜏𝜔=¹𝜔²+𝐴₁𝐵₁−𝐶₁𝐷₁𝜔𝜔⁴+𝐴²₁𝜔²−𝐷²₁𝐵,cos𝜏𝜔=¹−𝐴₁𝐶₁𝜔²+𝐵₁𝐷₁𝜔⁴+𝐴²₁𝜔²−𝐷²₁. (2.15) Then we have 𝜔⁸+𝑒₃𝜔⁶+𝑒₂𝜔⁴+𝑒₁𝜔²+𝑒₀=0, (2.16) where 𝑒₃=2𝐴²₁−𝐶²₁,𝑒₂=𝐴⁴₁+2𝐶²₁𝐷₁−𝐴²₁𝐶²₁−𝐵²₁−2𝐷²₁,𝑒₁=4𝐴₁𝐵₁𝐶₁𝐷₁−𝐴²₁𝐵²₁−𝐶²₁𝐷²₁−2𝐴²₁𝐷²₁−2𝐵²₁𝐷₁,𝑒₀=𝐷⁴₁−𝐵²₁𝐷²₁. (2.17) Let 𝑣=𝜔², then (2.16) becomes 𝑣⁴+𝑒₃𝑣³+𝑒₂𝑣²+𝑒₁𝑣+𝑒₀=0. (2.18) Next, we give the following assumption. (𝐻₂): (2.18) has at least one positive real root.

Suppose that (𝐻₂) holds. Without loss of generality, we assume that (2.18) has four real positive roots, which are defined by 𝑣₁, 𝑣₂, 𝑣₃, and 𝑣₄, respectively. Then (2.16) has four positive roots 𝜔_𝑘=√𝑣_𝑘,𝑘=1,2,3,4. Therefore, 𝜏_𝑘^(𝑗)=1𝜔_𝑘𝐵arccos¹−𝐴₁𝐶₁𝜔²_𝑘+𝐵₁𝐷₁𝜔⁴_𝑘+𝐴²₁𝜔²_𝑘−𝐷²₁+2𝑗𝜋,𝑘=1,2,3,4;𝑗=0,1,2… (2.19) Then we can know that ± 𝑖𝜔_𝑘 are a pair of purely imaginary roots of (2.11) with 𝜏=𝜏_𝑘^(𝑗). Define 𝜏₀=𝜏_𝑘⁽⁰⁾=min𝜏_𝑘⁽⁰⁾,𝜔₀=𝜔_𝑘0,𝑘=1,2,3,4. (2.20) Let 𝜆(𝜏)=𝛼(𝜏)+𝑖𝜔(𝜏) be the root of (2.11) near 𝜏=𝜏₀ which satisfies 𝛼(𝜏₀)=0,𝜔(𝜏₀)=𝜔₀. Taking the derivative of 𝜆 with respect to 𝜏 in (2.13), we obtain 𝐶₁𝑑𝜆+𝑑𝜏2𝜆+𝐴₁𝑒^𝜆𝜏𝑑𝜆+𝜆𝑑𝜏²+𝐴₁𝜆𝑒^𝜆𝜏𝜆+𝜏𝑑𝜆𝑑𝜏−𝐷₁𝑒^−𝜆𝜏𝜆+𝜏𝑑𝜆𝑑𝜏=0. (2.21) it follows that 𝑑𝜆=𝜆𝐷𝑑𝜏¹𝑒^−𝜆𝜏−𝜆²+𝐴₁𝜆𝑒^𝜆𝜏2𝜆+𝐴₁𝑒^𝜆𝜏+𝐶₁−𝜏𝐷₁𝑒^−𝜆𝜏−𝜆²+𝐴₁𝜆𝑒^𝜆𝜏. (2.22) Thus 𝑑𝜆𝑑𝜏⁻¹=2𝜆+𝐴₁𝑒^𝜆𝜏+𝐶₁𝐷₁𝜆𝑒^−𝜆𝜏−𝜆³+𝐴₁𝜆²𝑒^𝜆𝜏−𝜏𝜆. (2.23) Let Λ₁=𝐷₁𝜔₀−𝜔³₀sin𝜏₀𝜔₀+𝐴₁𝜔²₀cos𝜏₀𝜔₀,Λ₂=𝐷₁𝜔₀+𝜔³₀cos𝜏₀𝜔₀+𝐴₁𝜔²₀sin𝜏₀𝜔₀,Λ₃=𝐴₁cos𝜏₀𝜔₀−2𝜔₀sin𝜏₀𝜔₀+𝐶₁,Λ₄=𝐴₁sin𝜏₀𝜔₀−2𝜔₀cos𝜏₀𝜔₀. (2.24) Substitute 𝜆=𝑖𝜔₀(𝜔₀>0) into (2.23), we can get 𝜏Re𝑑𝜆₀𝑑𝜏⁻¹=Re2𝜆+𝐴₁𝑒^𝜆𝜏+𝐶₁𝐷₁𝜆𝑒^−𝜆𝜏−𝜆³+𝐴₁𝜆²𝑒^𝜆𝜏_{𝜆=𝑖𝜔0}=Λ₁×Λ₃+Λ₂×Λ₄Λ²₁+Λ²₂. (2.25) Noting that 𝜏signRe𝑑𝜆₀𝜏𝑑𝜏=sign𝑑Re𝜆₀𝑑𝜏⁻¹. (2.26) Therefore, we make the following assumption in order to give the main results: (𝐻₃)∶Λ₁×Λ₃+Λ₂×Λ₄≠0. Then, by Corollary 2.4 in [23, 24], we have the following theorem.

Theorem 2.1. For system (1.4), if the conditions (𝐻₁)-(𝐻₃) hold, then the equilibrium 𝐸_∗(𝑥_∗,𝑦_∗) of system (1.4) is asymptotically stable for 𝜏∈[0,𝜏₀) and unstable when 𝜏>𝜏₀. And system (1.4) has a branch of periodic solution bifurcation from the zero solution near 𝜏=𝜏₀.

Case 3. 𝜏₁≠𝜏₂, 𝜏₁>0 and 𝜏₂>0.
The associated characteristic equation of the system is 𝜆²+𝐴₂𝜆+𝐵₂𝜆𝑒^−𝜆𝜏1+𝐶₂+𝐷₂𝜆𝑒^−𝜆𝜏2+𝐸₂𝑒^{−𝜆(𝜏1+𝜏2)}, (2.27) where 𝐴₂=−𝑎₁₁,𝐵₂=−𝑏₁₁,𝐶₂=𝑎₁₁𝑐₂₂−𝑎₁₂𝑐₂₁,𝐷₂=−𝑐₂₂,𝐸₂=𝑏₁₁𝑐₂₂. (2.28) We consider (2.27) with 𝜏₂ in its stable interval, regarding 𝜏₁ as a parameter. Without loss of generality, we consider system (1.4) under the case considered in [16], and 𝜏₂∈[0,𝜏₂₀). 𝜏₂₀ is defined as in [16] and can be obtained by 𝜏₂₀=1𝜔₊𝑎arccos¹¹𝑐₂₂+𝑏₁₁𝑐₂₂−𝑎₁₂𝑐₂₁𝜔²₊−𝑎₁₁+𝑏₁₁𝑐₂₂𝜔²₊𝑎₁₁𝑐₂₂+𝑏₁₁𝑐₂₂−𝑎₁₂𝑐₂₁²+𝑐₂₂𝜔₊², (2.29) with 𝜔₊=−𝑎₁₁+𝑏₁₁²−𝑐²₂₂+𝑎₁₁+𝑏₁₁²−𝑐²₂₂²+4𝑎₁₁𝑐₂₂+𝑏₁₁𝑐₂₂−𝑎₁₂𝑐₂₁²2. (2.30) Let 𝜆=𝑖𝜔(𝜔>0) be a root of (2.27). Then we obtain 𝐵₂𝜔−𝐸₂sin𝜏₂𝜔sin𝜏₁𝜔+𝐸₂cos𝜏₂𝜔cos𝜏₁𝜔=𝜔²−𝐶₂cos𝜏₂𝜔−𝐷₂𝜔sin𝜏₂𝐵𝜔,₂𝜔−𝐸₂sin𝜏₂𝜔cos𝜏₁𝜔−𝐸₂cos𝜏₂𝜔sin𝜏₁𝜔=𝐶₂sin𝜏₂𝜔−𝐷₂𝜔cos𝜏₂𝜔−𝐴₂𝜔. (2.31) It follows from (2.31) that sin𝜏₁𝑀𝜔=¹𝑁₁−𝑀₂𝑁₂𝑀²₁+𝑀²₂,cos𝜏₁𝑀𝜔=¹𝑁₂+𝑀₂𝑁₁𝑀²₁+𝑀²₂, (2.32) With 𝑀₁=𝐵₂𝜔−𝐸₂sin𝜏₂𝜔,𝑀₂=𝐸₂cos𝜏₂𝑁𝜔,₁=𝜔²−𝐶₂cos𝜏₂𝜔−𝐷₂𝜔sin𝜏₂𝜔,𝑁₂=−𝐴₂𝜔+𝐶₂sin𝜏₂𝜔−𝐷₂𝜔cos𝜏₂𝜔. (2.33) Then we have 𝑃₁(𝜔)+𝑃₂(𝜔)sin𝜏₂𝜔+𝑃₃(𝜔)cos𝜏₂𝜔=0, (2.34) where 𝑃₁(𝜔)=𝜔⁴+𝐴²₂+𝐷²₂−𝐵²₂𝜔²+𝐶²₂−𝐸²₂,𝑃₂(𝜔)=−2𝐷₂𝜔³−2𝐴₂𝐶₂𝜔+2𝐵₂𝐸₂𝜔,𝑃₃(𝜔)=2𝐴₂𝐷₂−𝐶₂𝜔². (2.35) Suppose that (𝐻₄): (2.34) has at least finite positive roots. If (𝐻₄) holds, we define the roots of (2.34) as 𝜔₁,𝜔₂,…,𝜔_𝑘. Then, for every fixed 𝜔_𝑖(𝑖=1,2,…,𝑘), there exists a sequence {𝜏₁^(𝑗)_𝑖∣𝑗=1,2,…} which satisfies (2.34). Let 𝜏_1∗=min𝜏₁^(𝑗)_𝑖∣𝑖=1,2,…,𝑘,𝑗=0,1,…. (2.36) When 𝜏₁=𝜏_1∗, (2.27) has a pair of purely imaginary roots ±𝑖𝜔_∗ for 𝜏₂∈[0,𝜏₂₀).

To verify the transversality condition of Hopf bifurcation, we take the derivative of 𝜆 with respect to 𝜏₁ in (2.27), we can obtain 𝑑𝜆𝑑𝜏₁=𝜆𝑒^−𝜆𝜏1𝐵₂𝜆+𝐸₂𝑒^−𝜆𝜏22𝜆+𝐴₂+𝐵₂𝑒^−𝜆𝜏1+𝐷₂𝑒^−𝜆𝜏2−𝜏₁𝑒^−𝜆𝜏1𝐵₂𝜆+𝐸₂𝑒^−𝜆𝜏2−𝜏₂𝑒^−𝜆𝜏2𝐶₂+𝐷₂𝜆+𝐸₂𝑒^−𝜆𝜏1. (2.37) Thus 𝑑𝜆𝑑𝜏₁⁻¹=2𝜆+𝐴₂+𝐵₂𝑒^−𝜆𝜏1+𝐷₂𝑒^−𝜆𝜏2−𝜏₂𝑒^−𝜆𝜏2𝐶₂+𝐷₂𝜆+𝐸₂𝑒^−𝜆𝜏1𝜆𝑒^−𝜆𝜏1𝐵₂𝜆+𝐸₂𝑒^−𝜆𝜏2−𝜏₁𝜆. (2.38) Substitute 𝜆=𝑖𝜔_∗(𝜔_∗>0) into (2.38), we can get 𝜏Re𝑑𝜆_1∗𝑑𝜏₁⁻¹=Δ₁×Δ₃+Δ₂×Δ₄Δ²₁+Δ²₂, (2.39) where Δ₁=𝐸₂𝜔_∗cos𝜏₂𝜔_∗sin𝜏_1∗𝜔_∗−𝜔_∗cos𝜏_1∗𝜔_∗𝐵₂𝜔_∗−𝐸₂𝜔_∗sin𝜏₂𝜔_∗,Δ₂=𝐸₂𝜔_∗cos𝜏₂𝜔_∗cos𝜏_1∗𝜔_∗+𝜔_∗sin𝜏_1∗𝜔_∗𝐵₂𝜔_∗−𝐸₂𝜔_∗sin𝜏₂𝜔_∗,Δ₃=𝐴₂+𝐷₂−𝜏₂𝐶₂cos𝜏₂𝜔_∗−𝜏₂𝐷₂𝜔_∗sin𝜏₂𝜔_∗+𝜏₂𝐸₂sin𝜏₂𝜔_∗sin𝜏_1∗𝜔_∗+𝐵₂−𝜏₂𝐸₂cos𝜏₂𝜔_∗cos𝜏_1∗𝜔_∗,Δ₄=2𝜔_∗+𝜏₂𝐶₂+𝜏₂𝐷₂𝜔_∗−𝐷₂sin𝜏₂𝜔_∗+𝜏₂𝐸₂sin𝜏₂𝜔_∗cos𝜏_1∗𝜔_∗−𝐵₂−𝜏₂𝐸₂cos𝜏₂𝜔_∗sin𝜏_1∗𝜔_∗. (2.40) Next, we make the following assumption: (𝐻₅): Δ₁×Δ₃+Δ₂×Δ₄≠0.

Thus, by the discussion above and by the general Hopf bifurcation theorem for FDEs in Hale [25], we have the following results.

Theorem 2.2. For 𝜏₂∈(0,𝜏₂₀), 𝜏₂₀ is defined by (2.29). If the conditions (𝐻₄)-(𝐻₅) hold, then the equilibrium 𝐸_∗(𝑥_∗,𝑦_∗) of system (1.4) is asymptotically stable for 𝜏₁∈[0,𝜏_1∗) and unstable when 𝜏>𝜏_1∗. System (1.4) has a branch of periodic solution bifurcation from the zero solution near 𝜏=𝜏_1∗.

3. Direction and Stability of Bifurcated Periodic Solutions

In this section, we shall investigate the direction of the Hopf bifurcation and the stability of bifurcating periodic solution of system (1.4) w.r. to 𝜏₁ for 𝜏₂∈(0,𝜏₂₀), and 𝜏₂₀ is defined by (2.29). The idea employed here is the normal form and center manifold theory described in Hassard et al. [26]. Throughout this section, it is considered that system (1.4) undergoes the Hopf bifurcation at 𝜏₁=𝜏_1∗, 𝜏₂∈(0,𝜏₂₀) at 𝐸_∗(𝑥_∗,𝑦_∗). Let 𝜏₁=𝜏_1∗+𝜇,𝜇∈𝑅 so that the Hopf bifurcation occurs at 𝜇=0. Without loss of generality, we assume that 𝜏_2∗<𝜏_1∗, where 𝜏_2∗∈(0,𝜏₂₀).

Let 𝑢₁(𝑡)=𝑥(𝑡)−𝑥_∗, 𝑢₂(𝑡)=𝑦(𝑡)−𝑦_∗, and rescaling the time delay 𝑡→(𝑡/𝜏₁), Then system (1.4) can be transformed into an FDE in 𝐶=𝐶([−1,0],𝑅²) as: 𝑢̇(𝑡)=𝐿_𝜇𝑢_𝑡+𝐹𝜇,𝑢_𝑡, (3.1) where 𝑢(𝑡)=(𝑢₁(𝑡),𝑢₂(𝑡))^𝑇∈𝑅² and 𝐿_𝜇∶𝐶→𝑅², 𝐹∶𝑅×𝐶→𝑅² are given, respectively, by 𝐿_𝜇𝜙=𝜏_1∗+𝜇𝐴^′𝜙(0)+𝐶^′𝜙−𝜏_2∗𝜏_1∗+𝐵^′𝜙(),𝐹=𝜏−1𝜇,𝜙_1∗+𝜇𝑓₁,𝑓₂^𝑇, (3.2) with 𝐴^′=𝑎₁₁𝑎₁₂00,𝐵^′=𝑏₁₁000,𝐶^′=𝑐00₂₁𝑐₂₂,𝑓₁=𝑔₁𝜙²₁(0)+𝑔₂𝜙₁(0)𝜙₂(0)+𝑔₃𝜙²₂(0)+𝑔₄𝜙₁(0)𝜙₁(−1)+ℎ₁𝜙³₁(0)+ℎ₂𝜙²₁(0)𝜙₂(0)+ℎ₃𝜙₁(0)𝜙²₂(0)+ℎ₄𝜙³₂(0)𝑓+⋯,₂=𝑔^′₁𝜙²₁−𝜏_2∗𝜏_1∗+𝑔^′₂𝜙₁−𝜏_2∗𝜏_1∗𝜙₂(0)+𝑔^′₃𝜙₁−𝜏_2∗𝜏_1∗𝜙₂−𝜏_2∗𝜏_1∗,+𝑔^′₄𝜙₂(0)𝜙₂−𝜏_2∗𝜏_1∗+ℎ^′₁𝜙³₁−𝜏_2∗𝜏_1∗+ℎ^′₂𝜙²₁−𝜏_2∗𝜏_1∗𝜙₂(0),+ℎ^′₃𝜙²₁−𝜏_2∗𝜏_1∗𝜙₂−𝜏_2∗𝜏_1∗+ℎ₄𝜙³₂(0)+⋯, (3.3) where 𝑔₁=𝑏𝑦_∗𝑎₁+𝑐₁𝑦_∗𝑎₁+𝑏𝑥_∗+𝑐₁𝑦_∗³,𝑔₂𝑎=−²¹+𝑎₁𝑏𝑥_∗+𝑎₁𝑐₁𝑦_∗+2𝑏𝑐₁𝑥_∗𝑦_∗𝑎₁+𝑏𝑥_∗+𝑐₁𝑦_∗³,𝑔₃=𝑐₁𝑥_∗𝑎₁+𝑏𝑥_∗𝑎₁+𝑏𝑥_∗+𝑐₁𝑦_∗³,𝑔₄ℎ=−1,₁𝑏=−²𝑦_∗𝑎₁+𝑐₁𝑦_∗𝑎₁+𝑏𝑥_∗+𝑐₁𝑦_∗⁴,ℎ₂=𝑎²₁𝑏+𝑎₁𝑏²𝑥_∗+2𝑏²𝑐₁𝑥_∗𝑦_∗−𝑏𝑐²₁𝑦²_∗𝑎₁+𝑏𝑥_∗+𝑐₁𝑦_∗⁴,ℎ₃=𝑎²₁𝑐₁+𝑎₁𝑐²₁𝑦_∗+2𝑏𝑐²₁𝑥_∗𝑦_∗−𝑏²𝑐₁𝑥²_∗𝑎₁+𝑏𝑥_∗+𝑐₁𝑦_∗⁴,ℎ₄𝑐=−²¹𝑥_∗𝑎₁+𝑏𝑥_∗𝑎₁+𝑏𝑥_∗+𝑐₁𝑦_∗⁴,𝑔^′₁=−𝛽𝑦²_∗𝑥³_∗,𝑔^′₂=𝛽𝑦_∗𝑥²_∗,𝑔^′₃=𝛽𝑦_∗𝑥²_∗,𝑔^′₄𝛽=−𝑥_∗,ℎ^′₁=𝛽𝑦²_∗𝑥⁴_∗,ℎ^′₂=−𝛽𝑦_∗𝑥³_∗,ℎ^′₃=−𝛽𝑦_∗𝑥³_∗. (3.4) Using Riesz representation theorem, there exists a 2×2 matrix function 𝜂(𝜃,𝜇),𝜃∈[−1,0] whose elements are of bounded variation, such that 𝐿_𝜇𝜙=⁰₋₁𝑑𝜂𝜃,𝜇𝜙(𝜃),𝜙∈𝐶[−1,0],𝑅₂. (3.5) In fact, choosing 𝜂=𝜏𝜃,𝜇^1∗+𝜇𝐴^′+𝐶^′+𝐵^′𝜏,𝜃=0,_1∗+𝜇𝐶^′+𝐵^′−𝜏,𝜃∈^2∗𝜏₁,𝜏,0_1∗+𝜇𝐵^′𝜏,𝜃∈−1,−^2∗𝜏₁,0,𝜃=−1. (3.6) For 𝜙∈𝐶([−1,0]), we define 𝐴𝜇(𝜃)𝜙=𝑑𝜙∫𝑑𝜃,−1≤𝜃<0,₀₋₁𝜙(𝜃)𝑅𝜇𝐹𝑑𝜂𝜃,𝜇,𝜃=0,𝜙=0,−1≤𝜃<0,𝜇,𝜙,𝜃=0. (3.7) Then system (3.1) can be transformed into the following operator equation 𝑢̇(𝑡)=𝐴𝜇𝑢_𝑡+𝑅𝜇𝑢_𝑡, (3.8) where 𝑢_𝑡=𝑢(𝑡+𝜃)=(𝑢₁(𝑡+𝜃),𝑢₂(𝑡+𝜃)).

For 𝜑∈𝐶¹([0,1],(𝑅²)^∗), where (𝑅²)^∗ is the 2-dimensional space of row vectors, we further define the adjoint operator 𝐴^∗ of A(0): 𝐴^∗𝜑=−𝑑𝜑(𝑠)∫𝑑𝑠,0<𝑠≤1,₀₋₁𝜑(−𝜉)𝑑𝜂(𝜉,0),𝑠=0, (3.9) and a bilinear inner product: 𝜑(𝑠),𝜙(𝜃)=𝜑^𝑇(0)𝜙(0)−⁰_𝜃=−1^𝜃_𝜉=0𝜑^𝑇(𝜉−𝜃)𝑑𝜂(𝜃)𝜙(𝜉)𝑑𝜉, (3.10) where 𝜂(𝜃)=𝜂(𝜃,0).

Since ±𝑖𝜔_∗𝜏_1∗ are eigenvalues of 𝐴(0), they are also eigenvalues of 𝐴^∗. Let 𝑞(𝜃)=(1,𝑞₂)^𝑇𝑒^{𝑖𝜔∗𝜏1∗𝜃} be the eigenvectors of 𝐴(0) corresponding to 𝑖𝜔_∗𝜏_1∗ and 𝑞^∗(𝑠)=(1/𝜌)(1,𝑞^∗₂)^𝑇𝑒^{𝑖𝜔∗𝜏1∗𝑠} be the eigenvectors of 𝐴^∗ corresponding to −𝑖𝜔_∗𝜏_1∗. By a simple computation, we can get 𝑞₂=𝑖𝜔_∗−𝑎₁₁−𝑏₁₁𝑒^{𝑖𝜔∗𝜏1∗}𝑎₁₂,𝑞^∗₂=−𝑖𝜔_∗+𝑎₁₁+𝑏₁₁𝑒^{𝑖𝜔∗𝜏1∗}𝑐₂₁𝑒^{𝑖𝜔∗𝜏2∗},𝜌=1+𝑞₂𝑞^∗₂+𝑏₁₁𝜏_1∗𝑒^{−𝑖𝜔∗𝜏1∗}+𝑐₂₁𝜏_2∗𝑞^∗₂𝑒^{−𝑖𝜔∗𝜏2∗}+𝑐₂₂𝜏_2∗𝑞₂𝑞^∗₂𝑒^{−𝑖𝜔∗𝜏2∗}. (3.11) Then ⟨𝑞^∗,𝑞⟩=1, ⟨𝑞^∗,𝑞⟩=0.

In the remainder of this section, Following the algorithms given in [26] and using similar computation process to that in [16], we can get that the coefficients which will be used to determine the important qualities of the bifurcating periodic solutions, 𝑔₂₀=2𝜏_1∗𝜌𝑔₁+𝑔₂𝑞⁽²⁾(0)+𝑔₃𝑞⁽²⁾(0)²+𝑔₄𝑞⁽¹⁾()+−1𝑞^∗₂𝑔^′₁𝑞⁽¹⁾−𝜏_2∗𝜏_1∗²+𝑔^′₂𝑞⁽¹⁾−𝜏_2∗𝜏₁^∗𝑞⁽²⁾(0)+𝑔^′₃𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾−𝜏_2∗𝜏_1∗+𝑔^′₄𝑞⁽²⁾(0)𝑞⁽²⁾−𝜏_2∗𝜏_1∗,𝑔₁₁=𝜏_1∗𝜌2𝑔₁+𝑔₂𝑞⁽²⁾(0)+𝑞⁽²⁾(0)+2𝑔₃𝑞⁽²⁾(0)𝑞⁽²⁾(0)+𝑔₄𝑞⁽¹⁾(−1)+𝑞⁽¹⁾()+−1𝑞^∗₂2𝑔^′₁𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽¹⁾−𝜏_2∗𝜏_1∗+𝑔^′₂𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑔^′₃𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾−𝜏_2∗𝜏_1∗+𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾−𝜏_2∗𝜏_1∗+𝑔^′₄𝑞⁽²⁾(0)𝑞⁽²⁾−𝜏_2∗𝜏_1∗+𝑞⁽²⁾(0)𝑞⁽²⁾−𝜏_2∗𝜏_1∗,𝑔₀₂=2𝜏_1∗𝜌𝑔₁+𝑔₂𝑞⁽²⁾(0)+𝑔₃𝑞⁽²⁾(0)²+𝑔₄𝑞⁽¹⁾()+−1𝑞^∗₂𝑔^′₁𝑞⁽¹⁾−𝜏_2∗𝜏_1∗²+𝑔^′₂𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑔^′₃𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾−𝜏_2∗𝜏_1∗+𝑔^′₄𝑞⁽²⁾(0)𝑞⁽²⁾−𝜏_2∗𝜏_1∗,𝑔₂₁=2𝜏_1∗𝜌𝑔₁𝑊⁽¹⁾₂₀(0)+2𝑊⁽¹⁾₁₁(0)+𝑔₂12𝑊⁽²⁾₂₀(0)+𝑊⁽²⁾₁₁(0)+12𝑊⁽¹⁾₂₀(0)𝑞⁽²⁾(0)+𝑊⁽¹⁾₁₁(0)𝑞⁽²⁾(0)+𝑔₃𝑊⁽²⁾₂₀(0)𝑞⁽²⁾(0)+2𝑊⁽²⁾₁₁(0)𝑞⁽²⁾(0)+𝑔₄12𝑊⁽¹⁾₂₀(−1)+𝑊⁽¹⁾₁₁()+1−12𝑊⁽¹⁾₂₀(0)𝑞⁽¹⁾(−1)+𝑊⁽¹⁾₁₁(0)𝑞⁽¹⁾(−1)+3ℎ₁+ℎ₂2𝑞⁽²⁾(0)+𝑞⁽²⁾(0)+ℎ₃2𝑞⁽²⁾(0)𝑞⁽²⁾(0)+𝑞⁽²⁾(0)²+3ℎ₄𝑞⁽²⁾(0)²𝑞⁽²⁾(0)+𝑞^∗₂𝑔^′₁𝑊⁽¹⁾₂₀−𝜏_2∗𝜏_1∗𝑞⁽¹⁾−𝜏_2∗𝜏_1∗+2𝑊⁽¹⁾₁₁−𝜏_2∗𝜏_1∗𝑞⁽¹⁾−𝜏_2∗𝜏_1∗+𝑔^′₂12𝑊⁽¹⁾₂₀−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑊⁽¹⁾₁₁−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+12𝑊⁽²⁾₂₀(0)𝑞⁽¹⁾−𝜏_2∗𝜏_1∗+𝑊⁽²⁾₁₁(0)𝑞⁽¹⁾−𝜏_2∗𝜏_1∗+𝑔^′₃12𝑊⁽¹⁾₂₀−𝜏_2∗𝜏_1∗𝑞⁽²⁾−𝜏_2∗𝜏_1∗+𝑊⁽¹⁾₁₁−𝜏_2∗𝜏_1∗𝑞⁽²⁾−𝜏_2∗𝜏_1∗+12𝑊⁽²⁾₂₀−𝜏_2∗𝜏_1∗𝑞⁽¹⁾−𝜏_2∗𝜏_1∗+𝑔^′₄12𝑊⁽²⁾₂₀(0)𝑞⁽²⁾−𝜏_2∗𝜏_1∗+𝑊⁽²⁾₁₁(0)𝑞⁽²⁾−𝜏_2∗𝜏_1∗+12𝑊⁽²⁾₂₀−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑊⁽²⁾₁₁−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+3ℎ^′₁𝑞⁽¹⁾−𝜏_2∗𝜏_1∗²𝑞⁽¹⁾−𝜏_2∗𝜏_1∗+ℎ^′₂𝑞⁽¹⁾−𝜏_2∗𝜏_1∗²𝑞⁽²⁾(0)+2𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+ℎ^′₃𝑞⁽¹⁾−𝜏_2∗𝜏_1∗²𝑞⁽²⁾−𝜏_2∗𝜏_1∗+2𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾−𝜏_2∗𝜏_1∗, (3.12) with 𝑊₂₀(𝜃)=𝑖𝑔₂₀𝑞(0)𝜏_1∗𝜔_∗𝑒^{𝑖𝜏1∗𝜔∗𝜃}+𝑖𝑔₀₂𝑞(0)3𝜏_1∗𝜔_∗𝑒^{−𝑖𝜏1∗𝜔∗𝜃}+𝐸₂₀𝑒^{2𝑖𝜏1∗𝜔∗𝜃},𝑊₁₁(𝜃)=−𝑖𝑔₁₁𝑞(0)𝜏_1∗𝜔_∗𝑒^{𝑖𝜏1∗𝜔∗𝜃}+𝑖𝑔₁₁𝑞(0)𝜏_1∗𝜔_∗𝑒^{−𝑖𝜏1∗𝜔∗𝜃}+𝐸₁₁, (3.13) where 𝐸₂₀ and 𝐸₁₁ can be computed as the following equations, respectively, 𝐸₂₀𝐸=2⁽¹⁾²⁰𝐸⁽²⁾₂₀×2𝑖𝜔_∗−𝑎₁₁−𝑏₁₁𝑒^{−2𝑖𝜏1∗𝜔∗}−𝑎₁₂−𝑐₂₁𝑒^{−2𝑖𝜏2∗𝜔∗}2𝑖𝜔_∗−𝑐₂₂𝑒^{−2𝑖𝜏2∗𝜔∗−1},𝐸₁₁𝐸=−⁽¹⁾¹¹𝐸⁽²⁾₁₁×𝑎₁₁+𝑏₁₁𝑎₁₂𝑐₂₁𝑐₂₂⁻¹, (3.14) with 𝐸⁽¹⁾₂₀=𝑔₁+𝑔₂𝑞⁽²⁾(0)+𝑔₃𝑞⁽²⁾(0)²+𝑔₄𝑞⁽¹⁾(),𝐸−1₍₂₎₂₀=𝑔^′₁𝑞⁽¹⁾−𝜏_2∗𝜏_1∗²+𝑔^′₂𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑔^′₃𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾−𝜏_2∗𝜏_1∗+𝑔^′₄𝑞⁽²⁾(0)𝑞⁽²⁾−𝜏_2∗𝜏_1∗,𝐸⁽¹⁾₁₁=2𝑔₁+𝑔₂𝑞⁽²⁾(0)+𝑞⁽²⁾(0)+2𝑔₃𝑞⁽²⁾(0)𝑞⁽²⁾(0)+𝑔₄𝑞⁽¹⁾(−1)+𝑞⁽¹⁾(),𝐸−1₍₂₎₁₁=2𝑔^′₁𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽¹⁾−𝜏_2∗𝜏_1∗+𝑔^′₂𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑔^′₃𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾−𝜏_2∗𝜏_1∗+𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾−𝜏_2∗𝜏_1∗+𝑔^′₄𝑞⁽²⁾(0)𝑞⁽²⁾−𝜏_2∗𝜏_1∗+𝑞⁻⁽²⁾(0)𝑞⁽²⁾−𝜏_2∗𝜏_1∗. (3.15) Therefore, we can calculate the following values: 𝐶₁(0)=𝑖2𝜏_1∗𝜔_∗𝑔₁₁𝑔₂₀−2𝑔₁₁²−𝑔₀₂²3+𝑔₂₁2,𝜇₂𝐶=−Re₁(0)Re𝜆^′𝜏_1∗,𝛽₂=2Re𝐶₁(0,𝑇₂)𝐶=−Im₁(0)+𝜇₂Im𝜆^′𝜏_1∗𝜏_1∗𝜔_∗. (3.16) Based on the discussion above, we can obtain the following results.

Theorem 3.1. The direction of the Hopf bifurcation is determined by the sign of 𝜇₂: if 𝜇₂>0(𝜇₂<0), the Hopf bifurcation is supercritical (subcritical). The stability of bifurcating periodic solutions is determined by the sign of 𝛽₂: if 𝛽₂<0(𝛽₂>0), the bifurcating periodic solutions are stable (unstable). The period of the bifurcating periodic solutions is determined by the sign of 𝑇₂: if 𝑇₂>0(𝑇₂<0), the period of the bifurcating periodic solutions increases (decreases).

4. Numerical Example

In order to support the analytic results obtained above, we give some numerical simulations in this section. We only study the most important steady state, namely, the positive steady state. We consider the following system by taking the same coefficients as in [16]: 𝑑𝑥𝑑𝑡=𝑥1−𝑥𝑡−𝜏₁−𝑥𝑦,0.01+3𝑥+𝑦𝑑𝑦𝑦𝑑𝑡=𝑦3.5−2𝑡−𝜏₂𝑥𝑡−𝜏₂, (4.1) where we get the positive equilibrium 𝐸_∗(0.6328,1.1074). For system (4.1), we can get that 𝐴+𝐵+𝐷=3.9017>0, namely, the condition (𝐻₁) holds.

For 𝜏₁=𝜏₂=𝜏≠0. By a simple computation, we obtain that (2.18) has two positive roots: 𝑣₁=12.0256, 𝑣₂=0.3915. Thus, we know that the condition (𝐻₂) holds. Further, we get 𝜔₀=3.4678, 𝜏₀=0.4292. In addition, we have Λ₁×Λ₃+Λ₂×Λ₄=94.4826>0. Therefore, the condition (𝐻₃) is satisfied. Hence, from Theorem 2.1, we conclude that the positive equilibrium 𝐸_∗(0.6328,1.1074) is asymptotically stable when 𝜏∈[0,𝜏₀). The corresponding waveform and the phase plot are illustrated by Figure 1. When the delay 𝜏 passes through the critical value 𝜏₀ the positive equilibrium 𝐸_∗(0.6328,1.1074) will lose its stability and a Hopf bifurcation occurs, and a family of periodic solution bifurcates from the positive equilibrium 𝐸_∗(0.6328,1.1074). This property is illustrated by the numerical simulation in Figure 2.

Figure 1 

𝐸∗ is stable for 𝜏1=𝜏2=𝜏=0.4<𝜏0=0.4292 with initial value 0.62,2.35.

Figure 2 

𝐸∗ is unstable for 𝜏1=𝜏2=𝜏=0.44>𝜏0=0.4292 with initial value 0.62,2.35.

For 𝜏₁≠𝜏₂, 𝜏₁>0 and 𝜏₂>0. Regard 𝜏₁ as a parameter and let 𝜏₂=0.36∈[0,𝜏₂₀), and from (2.29) we can obtian 𝜏₂₀=0.4037. Then by a simple computation, we can obtain (2.34) has a positive root 𝜔_∗=0.6792. Thus, the condition (𝐻₄) holds. Further, we we can get 𝜏_1∗=2.3211 and Δ₁×Δ₃+Δ₂×Δ₄=3.7131>0. Namely, the condition (𝐻₅) holds. Thus, by Theorem 2.2, the positive equilibrium 𝐸_∗(0.6328,1.1074) is asymptotically stable when 𝜏₁∈[0,𝜏_1∗) and unstable when 𝜏₁>𝜏_1∗, and a family of periodic solution bifurcates from the positive equilibrium 𝐸_∗(0.6328,1.1074). The corresponding waveform and the phase plot are illustrated by Figures 3 and 4, respectively. In addition, from (3.16), we can get 𝐶₁(0)=−26.2632−7.5707𝑖, then 𝜇₂=12.1471, 𝛽₂=−52.5264, 𝑇₂=5.8586. From Theorem 3.1, we know that the Hopf bifurcation is supercritical and the bifurcating periodic solutions are stable.

Figure 3 

𝐸∗ is stable for 𝜏1=0.8<𝜏1∗=2.3211,𝜏2=0.36 with initial value 0.62,2.35.

Figure 4 

𝐸∗ is unstable for 𝜏1=2.4>𝜏1∗=2.3211,𝜏2=0.36 with initial value 0.62,2.35.

5. Conclusion

In this paper, a delayed predator-prey system with Beddington-DeAngelis functional response has been investigated. The bifurcation of a predator-prey system with single delay has been studied by many researchers [16, 27–30]. However, there are few papers considering the bifurcation of a predator-prey system with multiple delays (see [31–33]). Compared with the system considered in [16], the system in this paper accounts for not only the feedback delay of the prey density but also the feedback delay of the predator. The sufficient conditions for the stability of the positive equilibrium and the existence of periodic solutions via Hopf bifurcation at the positive equilibrium are obtained when 𝜏₁=𝜏₂ and 𝜏₁≠𝜏₂ with 𝜏₂∈[0,𝜏₂₀). Special attention is paid to the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions. By computation, we find that the feedback delay of the predator is marked because the critical value of 𝜏₂ is much smaller when we only consider it. The feedback delay of the prey is unremarkable because the critical value of 𝜏₁ is much bigger when we consider it with 𝜏₂ in its stable interval. Furthermore, Zhang [16] has obtained that the two species in system (1.4) with only the feedback delay of the predator could coexist. However, we get that the two species could also coexist with some available feedback delays of the prey and the predator. This is valuable from the view of ecology.

Unfortunately, the existence of the periodic solutions remain valid only in a small neighborhood of the critical value. It is definitely an interesting work to investigate whether these nonconstant periodic solutions which are obtained through local Hopf bifurcation can still exist for large values of the corresponding parameter time delay. The global continuation of the local Hopf bifurcation is left as the future work.

Acknowledgments

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