Hopf Bifurcation Analysis for a Semiratio-Dependent Predator-Prey System with Two Delays

Abstract

This paper is concerned with a semiratio-dependent predator-prey system with nonmonotonic functional response and two delays. It is shown that the positive equilibrium of the system is locally asymptotically stable when the time delay is small enough. Change of stability of the positive equilibrium will cause bifurcating periodic solutions as the time delay passes through a sequence of critical values. The properties of Hopf bifurcation such as direction and stability are determined by using the normal form method and center manifold theorem. Numerical simulations confirm our theoretical findings.

1. Introduction

It is well known that there are many factors which can affect dynamical properties of predator-prey systems. One of the important factors is the functional response describing the number of prey consumed per predator per unit time for given quantities of prey and predators. Numerous laboratory experiment, and observations have shown that a more suitable general predator-prey system should be based on the “ratio-dependent” theory, especially when predators have to search, share, or compete for food [1–3]. And ratio-dependent predator-prey systems have been investigated by many scholars [4–11]. In [4], Zhang and Lu considered the following semi-ratio-dependent predator-prey system with the nonmonotonic functional response 𝑑𝑥₁(𝑡)𝑑𝑡=𝑥₁(𝑡)×𝑟₁(𝑡)−𝑏₁(𝑡)𝑥₁(𝑡)−𝑎₁(𝑡)𝑥₂(𝑡)𝑚²+𝑛𝑥₁(𝑡)+𝑥²₁(𝑡),𝑑𝑥₂(𝑡)𝑑𝑡=𝑥₂(𝑡)𝑟₂(𝑡)−𝑎₂(𝑡)𝑥₂(𝑡)𝑥₁(𝑡), (1) where 𝑥₁(𝑡) and 𝑥₂(𝑡) denote the densities of the prey and the predator, respectively. 𝑟₁(𝑡) and 𝑟₂(𝑡) denote the intrinsic growth rates of the prey and the predator, respectively. 𝑏₁(𝑡) is the intraspecific competition rate of the prey. 𝑎₁(𝑡) is the capturing rate of the predator. 𝑎₂(𝑡) is a measure of the food quality that the prey provided for conversion into predator birth. The predator grows logistically with the carrying capacity 𝑥₁(𝑡)/𝑎₂(𝑡) proportional to the population size of the prey.

𝑚≠0,𝑛≥0, and 𝑎_𝑖(𝑡),𝑟_𝑖(𝑡),𝑏₁(𝑡),(𝑖=1,2) are assumed to be continuously positive periodic functions with period 𝜔. Zhang and Lu [4] established the existence of positive periodic solutions of system (1), and sufficient conditions for the uniqueness and global stability of the positive solutions of system (1) were also obtained by constructing a Lyapunov function.

Time delays of one type or another have been incorporated into predator-prey systems by many researchers since a time delay could cause a stable equilibrium to become unstable and cause the population to fluctuate [5, 6, 12–15]. In [5], Ding et al. incorporated the time delay 𝜏(𝑡)≥0 due to negative feedback of the prey into system (1) and got the special case (𝑛=0) of system (1) with time delay 𝑑𝑥₁(𝑡)𝑑𝑡=𝑥₁(𝑡)𝑟₁(𝑡)−𝑏₁(𝑡)𝑥₁((𝑡−𝑎𝑡−𝜏))¹(𝑡)𝑥₂(𝑡)𝑚²+𝑥²₁(𝑡),𝑑𝑥₂(𝑡)𝑑𝑡=𝑥₂(𝑡)𝑟₂(𝑡)−𝑎₂(𝑡)𝑥₂(𝑡)𝑥₁(𝑡). (2) Ding et al. [5] established the existence of a positive periodic solution for system (2) by using the continuation theorem of coincidence degree theory. Sufficient conditions for the permanence of system (2) were obtained by Li and Yang in [6]. But studies on predator-prey systems not only involve the persistence and the periodic phenomenon but also involve many patterns of other behavior such as global attractivity [16, 17] and bifurcation phenomenon [18–21]. Starting from this point, we will study the bifurcation phenomenon and the properties of periodic solutions of the following predator-prey system with two delays: 𝑑𝑥₁(𝑡)𝑑𝑡=𝑥₁(𝑡)𝑟₁−𝑏₁𝑥₁𝑡−𝜏₁−𝑎₁𝑥₂(𝑡)𝑚²+𝑛𝑥₁(𝑡)+𝑥²₁(𝑡),𝑑𝑥₂(𝑡)𝑑𝑡=𝑥₂(𝑡)𝑟₂−𝑎₂𝑥₂𝑡−𝜏₂𝑥₁𝑡−𝜏₂. (3) Unlike the assumptions in [4–6], we assume that all the parameters of system (3) are positive constants. 𝜏₁≥0 and 𝜏₂≥0 are the feedback delays of the prey and the predator, respectively.

The rest of this paper is organized as follows. In Section 2, sufficient conditions are obtained for the local stability of the positive equilibrium and the existence of Hopf bifurcation for possible combinations of the two delays in system (3). In Section 3, we give the formula determining the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions. Finally, numerical simulations supporting the theoretical analysis are also included.

2. Stability of the Positive Equilibrium and Local Hopf Bifurcations

It is not difficult to verify that system (2) has at least one positive equilibrium 𝐸^∗(𝑥^∗₁,𝑥^∗₂), where 𝑥^∗₂=𝑟₂𝑥^∗₁/𝑎₂ and 𝑥^∗₁ is the root of (4) 𝑏₁𝑥³₁+𝑛𝑏₁−𝑟₁𝑥²₁+𝑎₁𝑟₂𝑎₂−𝑛𝑟₁+𝑚²𝑏₁𝑥₁−𝑚²𝑟₁=0. (4)

Let 𝑥₁(𝑡)=𝑥₁(𝑡)−𝑥^∗₁, 𝑥₂(𝑡)=𝑥₂(𝑡)−𝑥^∗₂. Dropping the bars for the sake of simplicity, system (2) can be rewritten in the following system: 𝑑𝑥₁(𝑡)𝑑𝑡=𝑎₁₁𝑥₁(𝑡)+𝑎₁₂𝑥₂(𝑡)+𝑏₁₁𝑥₁𝑡−𝜏₁+𝑓₁,𝑑𝑥₂(𝑡)𝑑𝑡=𝑐₂₁𝑥₁𝑡−𝜏₂+𝑐₂₂𝑥₂𝑡−𝜏₂+𝑓₂, (5) where 𝑎₁₁=𝑎₁𝑥^∗₁𝑥^∗₂𝑛+2𝑥^∗₁𝑚²+𝑛𝑥^∗₁+𝑥^∗₁²²,𝑎₁₂𝑎=−¹𝑥^∗₁𝑚²+𝑛𝑥^∗₁+𝑥^∗₁²,𝑏₁₁=−𝑏₁𝑥^∗₁,𝑐₂₁=𝑎₂𝑥^∗₂²𝑥^∗₁²,𝑐₂₂𝑎=−²𝑥^∗₂𝑥^∗₁,𝑓₁=𝑎₁₃𝑥²₁(𝑡)+𝑎₁₄𝑥₁(𝑡)𝑥₂(𝑡)+𝑎₁₅𝑥₁(𝑡)𝑥₁𝑡−𝜏₁+𝑎₁₆𝑥²₁(𝑡)𝑥₂(𝑡)+𝑎₁₇𝑥³₁(𝑡)𝑓+⋯₂=𝑐₂₃𝑥²₁𝑡−𝜏₂+𝑐₂₄𝑥₁𝑡−𝜏₂𝑥₂(𝑡)+𝑐₂₅𝑥₂𝑡−𝜏₂𝑥₂(𝑡)+𝑐₂₆𝑥₁𝑡−𝜏₂𝑥₂𝑡−𝜏₂+𝑐₂₇𝑥³₁𝑡−𝜏₂+𝑐₂₈𝑥²₁𝑡−𝜏₂𝑥₂(𝑡)+𝑐₂₉𝑥²₁𝑡−𝜏₂𝑥₂𝑡−𝜏₂+𝑐₂₁₀𝑥₁𝑡−𝜏₂×𝑥₂𝑡−𝜏₂𝑥₂(𝑡)𝑎+⋯,₁₃=𝑎₁𝑥^∗₁𝑥^∗₂𝑚²−𝑛²−3𝑛𝑥^∗₁−3𝑥^∗₁²𝑚²+𝑛𝑥^∗₁+𝑥^∗₁²³,𝑎₁₄=𝑎₁𝑥^∗₁²−𝑚²𝑚²+𝑛𝑥^∗₁+𝑥^∗₁²²,𝑎₁₅=−𝑏₁,𝑎₁₆=𝑎₁𝑥^∗₁𝑚²−𝑛²−3𝑛𝑥^∗₁−3𝑥^∗₁²𝑚²+𝑛𝑥^∗₁+𝑥^∗₁²³,𝑎₁₇=𝑎₁𝑥^∗₂𝑚²−𝑛²−6𝑛𝑥^∗₁−9𝑥^∗₁²3𝑚²+𝑛𝑥^∗₁+𝑥^∗₁²³−𝑎₁𝑥^∗₁𝑥^∗₂𝑚²−𝑛²−3𝑛𝑥^∗₁−3𝑥^∗₁²𝑛+2𝑥^∗₁𝑚²+𝑛𝑥^∗₁+𝑥^∗₁²⁴,𝑐₂₃𝑎=−²𝑥^∗₂²𝑥^∗₁³,𝑐₂₄=𝑎₂𝑥^∗₂𝑥^∗₁²,𝑐₂₅𝑎=−²𝑥^∗₁,𝑐₂₆=𝑎₂𝑥^∗₂𝑥^∗₁²,𝑐₂₇=𝑎₂𝑥^∗₂²𝑥^∗₁⁴,𝑐₂₈𝑎=−²𝑥^∗₂𝑥^∗₁³,𝑐₂₉𝑎=−²𝑥^∗₂𝑥^∗₁³,𝑐₂₁₀=𝑎₂𝑥^∗₁². (6) The linearized system of (5) is 𝑑𝑥₁(𝑡)𝑑𝑡=𝑎₁₁𝑥₁(𝑡)+𝑎₁₂𝑥₂(𝑡)+𝑏₁₁𝑥₁𝑡−𝜏₁,𝑑𝑥₂(𝑡)𝑑𝑡=𝑐₂₁𝑥₁𝑡−𝜏₂+𝑐₂₂𝑥₂𝑡−𝜏₂. (7) The associated characteristic equation of system (7) is 𝜆²+𝐴𝜆+𝐵𝜆𝑒^−𝜆𝜏1+(𝐶𝜆+𝐷)𝑒^−𝜆𝜏2+𝐸𝑒^{−𝜆(𝜏1+𝜏2)}=0, (8) where 𝐴=−𝑎₁₁,𝐵=−𝑏₁₁,𝐶=−𝑐₂₂,𝐷=𝑎₁₁𝑐₂₂−𝑎₁₂𝑐₂₁,𝐸=𝑏₁₁𝑐₂₂. (9)

Case 1 (𝜏₁=𝜏₂=0). Equation (8) becomes 𝜆²+(𝐴+𝐵+𝐶)𝜆+𝐷+𝐸=0. (10) All the roots of (10) have negative real parts if and only if (𝐻₁)𝐴+𝐵+𝐶>0 and 𝐷+𝐸>0.Thus, the positive equilibrium 𝐸^∗(𝑥^∗₁,𝑥^∗₂) is locally stable when the condition (𝐻₁) holds.

Case 2 (𝜏₁>0,𝜏₂=0). Equation (8) becomes 𝜆²+(𝐴+𝐶)𝜆+𝐷+(𝐵𝜆+𝐸)𝑒^−𝜆𝜏1=0. (11) Let 𝜆=𝑖𝜔₁(𝜔₁>0) be the root of (11), and we have 𝐵𝜔₁sin𝜏₁𝜔₁+𝐸cos𝜏₁𝜔₁=𝜔²₁−𝐷,𝐵𝜔cos𝜏₁𝜔₁−𝐸sin𝜏₁𝜔₁=−(𝐴+𝐶)𝜔₁. (12) It follows that 𝜔⁴₁+(𝐴+𝐶)²−𝐵²−2𝐷𝜔²₁+𝐷²−𝐸²=0. (13) If condition (𝐻₂)𝐷−𝐸<0 holds, then (13) has only one positive root 𝜔₁₀=−(𝐴+𝐶)²−𝐵²+−2𝐷(𝐴+𝐶)²−𝐵²−2𝐷²−4𝐷²−𝐸²×2^−11/2. (14) Now, we can get the critical value of 𝜏₁ by substituting 𝜔₁₀ into (12). After computing, we obtain 𝜏_1𝑘=1𝜔₁₀()𝜔arccos𝐸−𝐴𝐵−𝐵𝐶²₁₀−𝐷𝐸𝐵²𝜔²₁₀+𝐸²+2𝑘𝜋𝜔₁₀,𝑘=0,1,2…. (15) Then, when 𝜏=𝜏_1𝑘, the characteristic equation (11) has a pair of purely imaginary roots ±𝑖𝜔₀.

Next, we will verify the transversality condition. Differentiating (11) with respect to 𝜏₁, we get 𝑑𝜆𝑑𝜏₁⁻¹=−2𝜆+𝐴+𝐶𝜆𝜆²+()+𝐵𝐴+𝐶𝜆+𝐷𝜆()−𝜏𝐵𝜆+𝐸𝜆¹. (16) Thus, Re𝑑𝜆𝑑𝜏_1𝜏⁻¹_1=𝜏10=2𝜔²₁₀+(𝐴+𝐶)²−𝐵²−2𝐷𝐵²𝜔²₁₀+𝐸²=(𝐴+𝐶)²−𝐵²−2𝐷²−4𝐷²−𝐸²𝐵²𝜔²₁₀+𝐸². (17) Obviously, if condition (𝐻₂) holds, then Re[𝑑𝜆/𝑑𝜏₁]_𝜏⁻¹_1=𝜏0>0. Note that sgn{Re[𝑑𝜆/𝑑𝜏₁]_𝜏⁻¹_1=𝜏10} = sgn{[𝑑Re(𝜆)/𝑑𝜏₁]_𝜏⁻¹_1=𝜏10}.

Therefore, if condition (𝐻₂) holds, the transversality condition is satisfied. From the analysis above, we have the following results.

Theorem 1. If conditions (𝐻₁) and (𝐻₂) hold, then the positive equilibrium 𝐸^∗(𝑥^∗₁,𝑥^∗₂) of system (3) is asymptotically stable for 𝜏₁∈[0,𝜏₁₀) and unstable when 𝜏₁>𝜏₁₀. System (3) undergoes a Hopf bifurcation at the positive equilibrium 𝐸^∗(𝑥^∗₁,𝑥^∗₂) when 𝜏₁=𝜏₁₀.

Case 3 (𝜏₁=0,𝜏₂>0). Equation (8) becomes 𝜆²+(𝐴+𝐵)𝜆+(𝐶𝜆+𝐷+𝐸)𝑒^−𝜆𝜏2=0. (18) Let 𝜆=𝑖𝜔₂(𝜔₂>0) be the root of (18), and then we have 𝐶𝜔₂sin𝜏₂𝜔₂+(𝐷+𝐸)cos𝜏₂𝜔₂=𝜔²₂,𝐶𝜔₂cos𝜏₂𝜔₂−(𝐷+𝐸)sin𝜏₂𝜔₂=−(𝐴+𝐵)𝜔₂, (19) which leads to 𝜔⁴₂+(𝐴+𝐵)²−𝐶²𝜔²₂−(𝐷+𝐸)²=0. (20) Obviously, if the condition (𝐻₁) holds, and then (20) has only one positive root 𝜔₂₀=−(𝐴+𝐵)²−𝐶²+(𝐴+𝐵)²−𝐶²²+4(𝐷+𝐸)²2. (21) The corresponding critical value of 𝜏₂ is 𝜏_2𝑘=1𝜔₂₀()𝜔arccos𝐷+𝐸−𝐴𝐶−𝐵𝐶²₁₀𝐶²𝜔²₂₀+(𝐷+𝐸)²+2𝑘𝜋𝜔₂₀,𝑘=0,1,2,…. (22) Similar to Case 2, differentiating (18) with respect to 𝜏₂, we get 𝑑𝜆𝑑𝜏₂⁻¹=−2𝜆+𝐴+𝐵𝜆𝜆²+()𝜆+𝐶𝐴+𝐵𝜆()−𝜏𝐶𝜆+𝐷+𝐸𝜆². (23) Then we can get Re𝑑𝜆𝑑𝜏_2𝜏⁻¹_2=𝜏20=2𝜔²₂₀+(𝐴+𝐵)²−𝐶²𝐶²𝜔²₂₀+(𝐷+𝐸)²=(𝐴+𝐵)²−𝐶²²−4(𝐷+𝐸)²𝐶²𝜔²₂₀+(𝐷+𝐸)². (24)

Obviously, if condition (𝐻₁) holds, then Re[𝑑𝜆/𝑑𝜏₂]_𝜏⁻¹_2=𝜏20>0. Namely, if condition (𝐻₁) holds, then the transversality condition is satisfied. Thus, we have the following results.

Theorem 2. If condition (𝐻₁) holds, then the positive equilibrium 𝐸^∗(𝑥^∗₁,𝑥^∗₂) of system (3) is asymptotically stable for 𝜏₂∈[0,𝜏₂₀) and unstable when 𝜏₂>𝜏₂₀. System (3) undergoes a Hopf bifurcation at the positive equilibrium 𝐸^∗(𝑥^∗₁,𝑥^∗₂) when 𝜏₂=𝜏₂₀.

Case 4 (𝜏₁=𝜏₂=𝜏>0). Equation (8) can be transformed into the following form: 𝜆²+𝐴𝜆+((𝐵+𝐶)𝜆+𝐷)𝑒^−𝜆𝜏+𝐸𝑒^−2𝜆𝜏=0. (25) Multiplying 𝑒^𝜆𝜏 on both sides of (25), we get (𝐵+𝐶)𝜆+𝐷+𝜆²+𝐴𝜆𝑒^𝜆𝜏+𝐸𝑒^−𝜆𝜏=0. (26) Let 𝜆=𝑖𝜔(𝜔>0) be the root of (26), then we obtain 𝐴𝜔sin𝜏𝜔+𝜔²𝜔−𝐸cos𝜏𝜔=𝐷,𝐴𝜔cos𝜏𝜔−²+𝐸sin𝜏𝜔=−(𝐵+𝐶)𝜔, (27) which follows that ()𝜔sin𝜏𝜔=𝐵+𝐶³+(𝐴𝐷−𝐵𝐸−𝐶𝐸)𝜔𝜔⁴+𝐴²𝜔²−𝐸²,()𝜔cos𝜏𝜔=𝐷−𝐴𝐵−𝐴𝐶²+𝐷𝐸𝜔⁴+𝐴²𝜔²−𝐸². (28) Then, we have 𝜔⁸+𝑐₃𝜔⁶+𝑐₂𝜔⁴+𝑐₁𝜔²+𝑐₀=0, (29) where 𝑐₀=𝐸⁴,𝑐₃=2𝐴²−(𝐵+𝐶)²,𝑐₁=2𝐷𝐸(𝐴𝐵+𝐴𝐶−𝐷)−2𝐴²𝐸²−(𝐴𝐷−𝐵𝐸−𝐶𝐸)²,𝑐₂=𝐴⁴−2𝐸²()−()−2𝐵+𝐶)(𝐴𝐷−𝐵𝐸−𝐶𝐸𝐷−𝐴𝐵−𝐴𝐶². (30) Denote 𝜔²=𝑣, and then (29) becomes 𝑣⁴+𝑐₃𝑣³+𝑐₂𝑣²+𝑐₁𝑣+𝑐₀=0. (31) Next, we give the following assumption.(𝐻₃) Equation (31) has at least one positive real root.

Suppose that (𝐻₃) holds. Without loss of generality, we assume that (31) has four real positive roots, which are denoted by 𝑣₁,𝑣₂,𝑣₃, and 𝑣₄. Then (29) has four positive roots 𝜔_𝑘=√𝑣_𝑘, 𝑘=1,2,3,4.

For every fixed 𝜔_𝑘, the corresponding critical value of time delay is 𝜏_𝑘^(𝑗)=1𝜔_𝑘()𝜔arccos𝐷−𝐴𝐵−𝐴𝐶²_𝑘+𝐷𝐸𝜔⁴_𝑘+𝐴²𝜔²_𝑘−𝐸²+2𝑗𝜋𝜔_𝑘,𝑘=1,2,3,4,𝑗=0,1,2,…. (32) Let 𝜏₀=min𝜏_𝑘⁽⁰⁾,𝑘=1,2,3,4.𝜔₀=𝜔_𝑘0. (33) Taking the derivative of 𝜆 with respect to 𝜏 in (26), we obtain 𝑑𝜆𝑑𝜏⁻¹=(2𝜆+𝐴)𝑒^𝜆𝜏+𝐵+𝐶𝜆𝐸𝑒^−𝜆𝜏−𝜆²+𝐴𝜆𝑒^𝜆𝜏−𝜏𝜆. (34) Thus, Re𝑑𝜆𝑑𝜏⁻¹_𝜏=𝜏0=𝑃_𝑅𝑄_𝑅+𝑃_𝐼𝑄_𝐼𝑄²_𝑅+𝑄²_𝐼, (35) with 𝑃_𝑅=𝐴cos𝜏₀𝜔₀−2𝜔₀sin𝜏₀𝜔₀𝑃+𝐵+𝐶,_𝐼=𝐴sin𝜏₀𝜔₀+2𝜔₀cos𝜏₀𝜔₀,𝑄_𝑅=𝐴𝜔²₀cos𝜏₀𝜔₀+𝐸𝜔₀−𝜔³₀sin𝜏₀𝜔₀,𝑄_𝐼=𝐴𝜔²₀sin𝜏₀𝜔₀+𝐸𝜔₀+𝜔³₀cos𝜏₀𝜔₀. (36)

It is easy to see that if condition (𝐻₄)𝑃_𝑅𝑄_𝑅+𝑃_𝐼𝑄_𝐼≠0, then Re[𝑑𝜆/𝑑𝜏]⁻¹_𝜏=𝜏0≠0. Therefore, we have the following results.

Theorem 3. If conditions (𝐻₁) and (𝐻₄) holds, then the positive equilibrium 𝐸^∗(𝑥^∗₁,𝑥^∗₂) of system (3) is asymptotically stable for 𝜏∈[0,𝜏₀) and unstable when 𝜏>𝜏₀. System (3) undergoes a Hopf bifurcation at the positive equilibrium 𝐸^∗(𝑥^∗₁,𝑥^∗₂) when 𝜏=𝜏₀.

Case 5 (𝜏₁≠𝜏₂,𝜏₁>0 and 𝜏₂>0). Considered the following: 𝜆²+𝐴𝜆+𝐵𝜆𝑒^−𝜆𝜏1+(𝐶𝜆+𝐷)𝑒^−𝜆𝜏2+𝐸𝑒^{−𝜆(𝜏1+𝜏2)}=0. (37)
We consider (8) with 𝜏₂ in its stable interval and regard 𝜏₁ as a parameter.
Let 𝜆=𝑖𝜔^∗₁(𝜔^∗₁>0) be the root of (8), and then we can obtain 𝐵𝜔^∗₁−𝐸sin𝜏₂𝜔^∗₁sin𝜏₁𝜔^∗₁+𝐸cos𝜏₂𝜔^∗₁cos𝜏₁𝜔^∗₁=𝜔^∗₁²−𝐷cos𝜏₂𝜔^∗₁−𝐶𝜔^∗₁sin𝜏₂𝜔^∗₁,𝐵𝜔^∗₁−𝐸sin𝜏₂𝜔^∗₁cos𝜏₁𝜔^∗₁−𝐸cos𝜏₂𝜔^∗₁cos𝜏₁𝜔^∗₁=−𝐴𝜔^∗₁+𝐷sin𝜏₂𝜔^∗₁−𝐶𝜔^∗₁cos𝜏₂𝜔^∗₁. (38) Then we get 𝜔^∗₁⁴+𝐴²−𝐵²+𝐶²𝜔^∗₁²+𝐷²−𝐸²+2𝐴𝐶−𝐷²𝜔^∗₁²cos𝜏₂𝜔^∗₁−2𝐶𝜔^∗₁³+𝐴𝐷𝜔^∗₁sin𝜏₂𝜔^∗₁=0. (39) Define 𝐹𝜔^∗₁=𝜔^∗₁⁴+𝐴²−𝐵²+𝐶²𝜔^∗₁²+𝐷²−𝐸²+2𝐴𝐶−𝐷²𝜔^∗₁²cos𝜏₂𝜔^∗₁−2𝐶𝜔^∗₁³+𝐴𝐷𝜔^∗₁sin𝜏₂𝜔^∗₁. (40)
If condition (𝐻₁)𝐷+𝐸>0 and (𝐻₂)𝐷−𝐸<0 hold, then 𝐹(0)=𝐷²−𝐸²<0. In addition, 𝐹(+∞)=+∞. Therefore, (39) has at least one positive root. We suppose that the positive roots of (39) are denoted as 𝜔^∗₁₁,𝜔^∗₁₂,…𝜔^∗_1𝑘. For every fixed 𝜔^∗_1𝑖, the corresponding critical value of time delay 𝜏^∗(𝑗)_1𝑖=1𝜔^∗_1𝑖𝑒×arccos¹𝜔^∗_1𝑖+𝑒₂𝜔^∗_1𝑖cos𝜏₂𝜔^∗_1𝑖+𝑒₃𝜔^∗_1𝑖sin𝜏₂𝜔^∗_1𝑖𝑒₄𝜔^∗_1𝑖+𝑒₅𝜔^∗_1𝑖sin𝜏₂𝜔^∗_1𝑖+2𝑗𝜋𝜔^∗_1𝑖𝑖=1,2,3,…,𝑘,𝑗=0,1,2,… (41) with 𝑒₁𝜔^∗_1𝑖=−𝐴𝐵𝜔^∗_1𝑖²−𝐷𝐸,𝑒₂𝜔^∗_1𝑖=(𝐸−𝐵𝐶)𝜔^∗_1𝑖²,𝑒₃𝜔^∗_1𝑖=(𝐴𝐸+𝐵𝐷)𝜔^∗_1𝑖,𝑒₄𝜔^∗_1𝑖=𝐵²𝜔^∗_1𝑖²+𝐸²,𝑒₅𝜔^∗_1𝑖=−2𝐵𝐸𝜔^∗_1𝑖. (42) Let 𝜏_1∗=min{𝜏^∗(0)_1𝑖}, 𝑖=1,2,3,…𝑘. When 𝜏₁=𝜏_1∗, (8) has a pair of purely imaginary roots ±𝑖𝜔_1∗ for 𝜏₂∈(0,𝜏₂₀).

Next, in order to give the main results, we give the following assumption(𝐻₅)𝑑Re𝜆(𝜏^∗₁₀)/𝑑𝜏₁≠0.

Therefore, we have the following results on the stability and bifurcation in system (3).

Theorem 4. If the conditions (𝐻₁) and (𝐻₅) hold and 𝜏₂∈(0,𝜏₂₀), then the positive equilibrium 𝐸^∗(𝑥^∗₁,𝑥^∗₂) of system (3) is asymptotically stable for 𝜏₁∈[0,𝜏_1∗) and unstable when 𝜏₁>𝜏_1∗. System (3) undergoes a Hopf bifurcation at the positive equilibrium 𝐸^∗(𝑥^∗₁,𝑥^∗₂) when 𝜏₁=𝜏_1∗.

3. Direction and Stability of the Hopf Bifurcation

In this section, we will investigate the direction of Hopf bifurcation and stability of the bifurcating periodic solutions of system (3) with respect to 𝜏₁ for 𝜏₂∈(0,𝜏₂₀). We assume that 𝜏_2∗<𝜏_1∗ where 𝜏_2∗∈(0,𝜏₂₀).

Let 𝜏₁=𝜏_1∗+𝜇,𝜇∈𝑅 so that 𝜇=0 is the Hopf bifurcation value of system (3). Rescale the time delay by 𝑡→(𝑡/𝜏₁). Let 𝑢₁(𝑡)=𝑥₁(𝑡)−𝑥^∗₁, 𝑢₂(𝑡)=𝑥₂(𝑡)−𝑥^∗₂, and then system (3) can be rewritten as an PDE in 𝐶=𝐶([−1,0],𝑅²) as 𝑢̇(𝑡)=𝐿_𝜇𝑢_𝑡+𝐹𝜇,𝑢_𝑡, (43) where 𝑢(𝑡)=𝑢₁(𝑡),𝑢₂(𝑡)^𝑇∈𝑅²,𝐿_𝜇∶𝐶⟶𝑅²,𝐹∶𝑅×𝐶⟶𝑅² (44) are given, respectively, by 𝐿_𝜇𝜙=𝜏_1∗+𝜇𝐴^′𝜙(0)+𝐶^′𝜙−𝜏_2∗𝜏_1∗+𝐵^′𝜙(),𝐹−1𝜇,𝑢_𝑡=𝜏_1∗+𝜇𝐹₁,𝐹₂^𝑇 (45) with 𝜙(𝜃)=𝜙₁(𝜃),𝜙₂(𝜃)^𝑇∈𝐶[−1,0],𝑅²,𝐴^′=𝑎₁₁𝑎₁₂00,𝐵^′=𝑏₁₁0,𝐶00_′=𝑐00₂₁𝑐₂₂,𝐹₁=𝑎₁₃𝜙²₁(0)+𝑎₁₄𝜙₁(0)𝜙₂(0)+𝑎₁₅𝜙₁(0)𝜙₁(−1)+𝑎₁₆𝜙²₁(0)𝜙₂(0)+𝑎₁₇𝜙³₁(0)𝐹+⋯₂=𝑐₂₃𝜙²₁−𝜏_2∗𝜏_1∗+𝑐₂₄𝜙₁−𝜏_2∗𝜏_1∗𝜙₂(0)+𝑐₂₅𝜙₂−𝜏_2∗𝜏_1∗𝜙₂(0)+𝑐₂₆𝜙₁−𝜏_2∗𝜏_1∗𝜙₂−𝜏_2∗𝜏_1∗+𝑐₂₇𝜙³₁−𝜏_2∗𝜏_1∗+𝑐₂₈𝜙²₁−𝜏_2∗𝜏_1∗𝜙₂(0)+𝑐₂₉𝜙²₁−𝜏_2∗𝜏_1∗𝜙₂−𝜏_2∗𝜏_1∗+𝑐₂₁₀𝜙₁−𝜏_2∗𝜏_1∗𝜙₂−𝜏_2∗𝜏_1∗𝜙₂(0)+⋯. (46) By the Riesz representation theorem, there is a 2×2 matrix function with bounded variation components 𝜂(𝜃,𝜇), 𝜃∈[−1,0] such that 𝐿_𝜇𝜙=⁰₋₁𝑑𝜂𝜃,𝜇𝜙(𝜃),𝜙∈𝐶. (47) In fact, we have chosen 𝜂=𝜏𝜃,𝜇^1∗+𝜇𝐴^′+𝐵^′+𝐶^′𝜏,𝜃=0,_1∗+𝜇𝐵^′+𝐶^′−𝜏,𝜃∈^2∗𝜏_1∗,𝜏,0_1∗+𝜇𝐵^′𝜏,𝜃∈−1,−^2∗𝜏_1∗,0,𝜃=−1. (48) For 𝜙∈𝐶([−1,0]),𝑅², we define 𝐴𝜇(𝜃)𝜙=𝑑𝜙∫𝑑𝜃,−1≤𝜃<0,₀₋₁𝜙(𝜃)𝑅𝜇𝐹𝑑𝜂𝜃,𝜇,𝜃=0,𝜙=0,−1≤𝜃<0,𝜇,𝜙,𝜃=0. (49) Then system (43) is equivalent to the following operator equation: 𝑢̇(𝑡)=𝐴𝜇𝑢_𝑡+𝑅𝜇𝑢_𝑡, (50) where 𝑢_𝑡=𝑢(𝑡+𝜃) for 𝜃∈[−1,0].

For 𝜑∈𝐶([−1,0]),(𝑅²)^∗, we define the adjoint operator 𝐴^∗ of 𝐴 as 𝐴^∗𝜇−(𝑠)𝑠=𝑑𝜑∫𝑑𝑠,0<𝑠≤1,₀₋₁𝑑𝜂^𝑇𝑠,𝜇𝜑(−𝑠),𝑠=0, (51) and a bilinear inner product 𝜑,𝜙=𝜑(0)𝜙(0)−⁰_𝜃=−1^𝜃_𝜉=0𝜑(𝜉−𝜃)𝑑𝜂(𝜃)𝜙(𝜉)𝑑𝜉, (52) where 𝜂(𝜃)=𝜂(𝜃,0).

Then, 𝐴(0) and 𝐴^∗(0) are adjoint operators. From the discussion above, we can know that ±𝑖𝜔_1∗𝜏_1∗ are the eigenvalues of 𝐴(0) and they are also eigenvalues of 𝐴^∗(0).

Let 𝑞(𝜃)=(1,𝑞₂)^𝑇𝑒^{𝑖𝜔1∗𝜏1∗𝜃} be the eigenvector of 𝐴(0) corresponding to 𝑖𝜔_1∗𝜏_1∗, and 𝑞^∗(𝜃) = 𝐷(1,𝑞^∗₂)𝑒^{𝑖𝜔1∗𝜏1∗𝑠} be the eigenvector of 𝐴^∗(0) corresponding to −𝑖𝜔_1∗𝜏_1∗. Then we have 𝐴(0)𝑞(𝜃)=𝑖𝜔_1∗𝜏_1∗𝑞(𝜃),𝐴^∗(0)𝑞^∗(𝜃)=−𝑖𝜔_1∗𝜏_1∗𝑞^∗(𝜃). (53) By a simple computation, we can get 𝑞₂=𝑖𝜔_1∗−𝑎₁₁−𝑏₁₁𝑒^{−𝑖𝜔1∗𝜏1∗}𝑎₁₂,𝑞^∗₂=𝑖𝜔_1∗−𝑎₁₁−𝑏₁₁𝑒^{𝑖𝜔1∗𝜏1∗}𝑐₂₁𝑒^{𝑖𝜔1∗𝜏2∗}𝐷=1+𝑞₂𝑞^∗₂+𝑏₁₁𝜏_1∗𝑒^{−𝑖𝜔1∗𝜏1∗}+𝑞^∗₂𝑐₂₁+𝑐₂₂𝑞₂𝜏_2∗𝑒^{−𝑖𝜔1∗𝜏2∗−1}. (54) Then, 𝑞^∗,𝑞=1,𝑞^∗,𝑞=0. (55)

Following the algorithms explained in the work of Hassard et al. [22] and using a computation process similar to that in [18], we can get the following coefficients which can determine the properties of Hopf bifurcation: 𝑔₂₀=2𝜏_1∗𝐷𝑎₁₃+𝑎₁₄𝑞⁽²⁾(0)+𝑎₁₅𝑞⁽¹⁾()+−1𝑞^∗₂𝑐₂₃𝑞⁽¹⁾−𝜏_2∗𝜏_1∗²+𝑐₂₄𝑞⁽¹⁾−𝜏_2∗𝜏_1∗×𝑞⁽²⁾(0)+𝑐₂₅𝑞⁽²⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑐₂₆𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾−𝜏_2∗𝜏_1∗,𝑔₁₁=𝜏_1∗𝐷2𝑎₁₃+𝑎₁₄𝑞⁽²⁾(0)+𝑞⁽²⁾(0)+𝑎₁₅𝑞⁽¹⁾(−1)+𝑞⁽¹⁾(+−1)𝑞^∗₂2𝑐₂₃𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽¹⁾−𝜏_2∗𝜏_1∗+𝑐₂₄𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑐₂₅𝑞⁽²⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑞⁽²⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑐₂₆𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾−𝜏_2∗𝜏_1∗+𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾−𝜏_2∗𝜏_1∗,𝑔₀₂=2𝜏_1∗𝐷𝑎₁₃+𝑎₁₄𝑞⁽²⁾(0)+𝑎₁₅𝑞⁽¹⁾()+−1𝑞^∗₂𝑐₂₃𝑞⁽¹⁾−𝜏_2∗𝜏_1∗²+𝑐₂₄𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑐₂₅𝑞⁽²⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑐₂₆𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾−𝜏_2∗𝜏_1∗,𝑔₂₁=2𝜏_1∗𝐷𝑎₁₃𝑊⁽¹⁾₁₁(0)+𝑊⁽¹⁾₂₀(0)+𝑎₁₄𝑊⁽¹⁾₁₁(0)𝑞⁽²⁾(0)+12𝑊⁽¹⁾₂₀(0)𝑞⁽²⁾(0)+𝑊⁽²⁾₁₁(0)+12𝑊⁽²⁾₂₀(0)+𝑎₁₅𝑊⁽¹⁾₁₁(0)𝑞⁽¹⁾()+1−12𝑊⁽¹⁾₂₀(0)𝑞⁽¹⁾(−1)+𝑊⁽¹⁾₁₁()+1−12𝑊⁽¹⁾₂₀(−1)+𝑎₁₆𝑞⁽²⁾(0)+2𝑞⁽²⁾(0)+3𝑎₁₇𝑞⁽¹⁾(0)+𝑞^∗₂𝑐₂₃2𝑊⁽¹⁾₁₁−𝜏_2∗𝜏_1∗𝑞⁽¹⁾−𝜏_2∗𝜏_1∗+𝑊⁽¹⁾₂₀−𝜏_2∗𝜏_1∗𝑞⁽¹⁾−𝜏_2∗𝜏_1∗+𝑐₂₄𝑊⁽¹⁾₁₁−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+12𝑊⁽¹⁾₂₀−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑊⁽²⁾₁₁(0)𝑞⁽¹⁾−𝜏_2∗𝜏_1∗+12𝑊⁽²⁾₂₀(0)𝑞⁽¹⁾−𝜏_2∗𝜏_1∗+𝑐₂₅𝑊⁽²⁾₁₁−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+12𝑊⁽²⁾₂₀−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑊⁽²⁾₁₁(0)𝑞⁽²⁾−𝜏_2∗𝜏_1∗+12𝑊⁽²⁾₂₀(0)𝑞⁽²⁾−𝜏_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∗²𝑞⁽²⁾(0)+2𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)𝑞⁽¹⁾−𝜏_2∗𝜏_1∗+𝑐₂₉𝑞⁽¹⁾−𝜏_2∗𝜏_1∗²𝑞⁽²⁾−𝜏_2∗𝜏_1∗+2𝑞⁽¹⁾×−𝜏_2∗𝜏_1∗𝑞⁽²⁾−𝜏_2∗𝜏_1∗𝑞⁽¹⁾−𝜏_2∗𝜏_1∗+𝑐₂₁₀𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾−𝜏_2∗𝜏_1∗×𝑞⁽²⁾(0)+𝑞⁽¹⁾−𝜏_2∗𝜏_1∗×𝑞⁽²⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0), (56) 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∗𝜃}+𝐸₂, (57) where 𝐸₁ and 𝐸₂ can be determined by the following equations, respectively: 2𝑖𝜔_1∗−𝑎₁₁−𝑏₁₁𝑒^{−2𝑖𝜔1∗𝜏1∗}−𝑎₁₂−𝑐₂₁𝑒^{−2𝑖𝜔1∗𝜏2∗}2𝑖𝜔_1∗−𝑐₃₃𝑒^{−2𝑖𝜔1∗𝜏2∗}𝐸₁𝐸=2¹⁽¹⁾𝐸₁⁽²⁾,𝑎₁₁+𝑏₁₁𝑎₁₂𝑐₂₁𝑐₂₂𝐸₂𝐸=−²⁽¹⁾𝐸₂⁽²⁾, (58) with 𝐸₁⁽¹⁾=𝑎₁₃+𝑎₁₄𝑞⁽²⁾(0)+𝑎₁₅𝑞⁽¹⁾(),𝐸−1₁₍₁₎=𝑐₂₃𝑞⁽¹⁾−𝜏_2∗𝜏_1∗²+𝑐₂₄𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑐₂₅𝑞⁽²⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑐₂₆𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾−𝜏_2∗𝜏_1∗,𝐸₂⁽¹⁾=2𝑎₁₃+𝑎₁₄𝑞⁽²⁾(0)+𝑞⁽²⁾(0)+𝑎₁₅𝑞⁽¹⁾(−1)+𝑞⁽¹⁾(),𝐸−1₂₍₂₎=2𝑐₂₃𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽¹⁾−𝜏_2∗𝜏_1∗+𝑐₂₄𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑐₂₅𝑞⁽²⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑞⁽²⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾(0)+𝑐₂₆𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾−𝜏_2∗𝜏_1∗+𝑞⁽¹⁾−𝜏_2∗𝜏_1∗𝑞⁽²⁾−𝜏_2∗𝜏_1∗. (59) Then, we can get the following coefficients: 𝐶₁(0)=𝑖2𝜔_1∗𝜏_1∗𝑔₁₁𝑔₂₀−2𝑔₁₁²−𝑔₀₂²3+𝑔₂₁2,𝜇₂𝐶=−Re₁(0)Re𝜆^′𝜏_1∗,𝛽₂=2Re𝐶₁(0,𝑇₂)𝐶=−Im₁(0)+𝜇₂Im𝜆^′𝜏_1∗𝜔_1∗𝜏_1∗. (60)

In conclusion, we have the following results.

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

4. Numerical Simulation and Discussions

To demonstrate the algorithm for determining the existence of Hopf bifurcation in Section 2 and the properties of Hopf bifurcation in Section 3, we give an example of system (3) in the following form: 𝑑𝑥₁(𝑡)𝑑𝑡=𝑥₁(𝑡)4−𝑥₁𝑡−𝜏₁−0.5𝑥₂(𝑡)4+𝑥₁(𝑡)+𝑥²₁(𝑡),𝑑𝑥₂(𝑡)𝑑𝑡=𝑥₂(𝑡)3−1.5𝑥₂𝑡−𝜏₂𝑥₁𝑡−𝜏₂, (61) which has a positive equilibrium 𝐸^∗(3.8298,7.6596). By calculation, we have 𝐴+𝐵+𝐶=6.5788, 𝐷+𝐸=11.2470, and 𝐷−𝐸=−11.7318. Namely, conditions (𝐻₁) and (𝐻₂) hold.

For 𝜏₁>0, 𝜏₂=0. We can get 𝜔₁₀=3.9074, 𝜏₁₀=0.3788. By Theorem 1, we know that when 𝜏₁∈[0,𝜏₁₀), the positive equilibrium 𝐸^∗(3.8298,7.6596) is locally asymptotically stable and is unstable if 𝜏₁>𝜏₁₀. As can be seen from Figure 1, when we let 𝜏₁=0.35∈[0,𝜏₁₀), the positive equilibrium 𝐸^∗(3.8298,7.6596) is locally asymptotically stable. However, when 𝜏₁=0.395>𝜏₁₀, the positive equilibrium 𝐸^∗(3.8298,7.6596) is unstable and a Hopf bifurcation occurs and a family of periodic solutions bifurcate from the positive equilibrium 𝐸^∗(3.8298,7.6596), which can be illustrated by Figure 2. Similarly, we can get 𝜔₂₀=3.0827,𝜏₂₀=0.4912. The corresponding waveform and the phase plots are shown in Figures 3 and 4.

(a)

(b)

(c)

Figure 1 

𝐸∗ is asymptotically stable for 𝜏1=0.350<𝜏10=0.3788.

(a)

(b)

(c)

Figure 2 

𝐸∗ is unstable for 𝜏1=0.395>𝜏10=0.3788.

(a)

(b)

(c)

Figure 3 

𝐸∗ is asymptotically stable for 𝜏2=0.450<𝜏20=0.4912.

(a)

(b)

(c)

Figure 4 

𝐸∗ is unstable for 𝜏2=0.505>𝜏20=0.4912.

For 𝜏₁=𝜏₂=𝜏>0. We can obtain 𝜔₀=3.4601,𝜏₀=0.3361. By Theorem 3, we know that when the time delay 𝜏 increases from zero to the critical value 𝜏₀=0.3361, the positive equilibrium 𝐸^∗(3.8298,7.6596) is locally asymptotically stable. It will lose its stability and a Hopf bifurcation occurs once 𝜏>𝜏₀=0.3361. This property can be illustrated by Figures 5 and 6.

(a)

(b)

(c)

Figure 5 

𝐸∗ is asymptotically stable for 𝜏=0.320<𝜏0=0.3361.

(a)

(b)

(c)

Figure 6 

𝐸∗ is unstable for 𝜏=0.355>𝜏0=0.3361.

Finally, consider 𝜏₁ as a parameter and let 𝜏₂=0.3∈(0,𝜏₂₀). We can get 𝜔_1∗=2.3761,𝜏_1∗=0.3490,𝜆^′(𝜏_1∗)=8.2230+2.6054𝑖. From Theorem 4, the positive equilibrium 𝐸^∗(3.8298,7.6596) is locally asymptotically stable for 𝜏₁∈[0,0.3490) and unstable when 𝜏₁>𝜏_1∗=0.3490, which can be seen from Figures 7 and 8.

(a)

(b)

(c)

Figure 7 

𝐸∗ is asymptotically stable for 𝜏1=0.325<𝜏1∗=0.3490 at 𝜏2=0.3.

(a)

(b)

(c)

Figure 8 

𝐸∗ is unstable for 𝜏1=0.360>𝜏1∗=0.3490 at 𝜏2=0.3.

In addition, from (60), we can get 𝐶₁(0)=−41.5020−17.8366𝑖. Furthermore, we have 𝜇₂=5.0471>0, 𝛽₂=−83.0040<0,𝑇₂=5.6519>0. Therefore, from Theorem 5, we know that the Hopf bifurcation is supercritical. The bifurcating periodic solutions are stable and increase. If the bifurcating periodic solutions are stable, then the prey and the predator may coexist in an oscillatory mode. Therefore, the two species in system (3) may coexist in an oscillatory mode under certain conditions.

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.