Hopf Bifurcation Analysis on General Gause-Type Predator-Prey Models with Delay

Abstract

A class of three-dimensional Gause-type predator-prey model with delay is considered. Firstly, a group of sufficient conditions for the existence of Hopf bifurcation is obtained via employing the polynomial theorem by analyzing the distribution of the roots of the associated characteristic equation. Secondly, the direction of the Hopf bifurcation and the stability of the bifurcated periodic solutions are determined by applying the normal form method and the center manifold theorem. Finally, some numerical simulations are carried out to illustrate the obtained results.

1. Introduction

Multispecies predator-prey models arise frequently on various ecosystems [1–4]. They usually exhibit more rich and complex dynamics as the number of species increases. Take three-species models, for example. There might be a food chain, or two predators or preys, whose relations might be cooperative or competitive. It has attracted extensive studies on the dynamics of various multispecies predator-prey models. Cheng et al. [5] derived some results to ensure the global stability of a predator-prey system with Holling’s type III functional response. klebanoff and Hastings [6] investigated a three-species food chain model. It showed that the system exhibits rich complexity features such as stable, periodic, and chaotic dynamics. Hsu and Huang [7] deals with the question of global stability of the positive locally asymptotically stable equilibrium in a class of predator-prey systems by using the Dulac’s criterion and constructing Lyapunov functions. See [8–14].

In 1977, Freedman and Waltman [12] investigated the following general Gause-type predator-prey model:𝑑𝑥(𝑡)(𝑥)(𝑥),(𝑡)𝑑𝑡=𝑥g−𝑦𝑝𝑑𝑦(𝑥𝑦,(𝑡)𝑑𝑡=𝑦−ℎ+𝑒𝑝)−𝑧𝑞𝑑𝑧𝑦,𝑑𝑡=𝑧−𝑠+𝑚𝑞 (1.1) where 𝑥(𝑡), 𝑦(𝑡), and 𝑧(𝑡) are the population densities of prey, predator, and top predator at time 𝑡, respectively. 𝑔(𝑥) is the intrinsic growth rate of prey; 𝑝(𝑥) and 𝑞(𝑦) are the specific growth rates of predator and top predator; ℎ, 𝑠>0 are the death rates of 𝑦(𝑡) and 𝑧(𝑡); 𝑒, 𝑚>0 are the conversion rates for prey and predator. They established the stability criteria and argued that the unique interior equilibrium exists and is locally asymptotically stable. Ginoux et al. [13] highlighted that this model has several Hopf bifurcations and a period-doubling cascade generating a snail shell-shaped chaotic attractor. They compared the quantity and property of the equilibria, bifurcation structure, and shape of attractors among this model and the so-called Rosenzweig-MacArthur and Hastings-Powell models and also gave the bifurcation analysis for each model. Hastings and Powell [14] discussed a continuous time model of a food chain incorporating nonlinear functional responses, and the model exhibits chaotic in long-term behavior when appropriate biologically reasonable parameter values are chosen. They found that, for different values of the key parameter, the system exhibits several types of asymptotic motions, namely, stable equilibrium point, limit cycles, change in periodicity of these cycles, and the so called teacup chaos.

We know that the population outbreak may happen for the species with periodic fluctuation. The outbreaks of pests and mice are famous. For example, Finerty discovered in Canada the populations of polar rabbits, lynxs and Ondatrazibethicas have 10-year cycle fluctuations, and the populations of lemmings and some murine experience 4-year cycle fluctuations. When the peak arrived, the vegetation would be severely damaged. Some cyclical fluctuations are consistent with the periodicity of intensity of certain infectious diseases. Thus, it is of great significance to study periodic solutions of biological systems for controlling insect. It is known that the delay differential equations (DDEs) will exhibit much more complicated dynamics than ordinary differential equations (ODEs) such as the existence of Bogdanov-Takens bifurcation and even chaos. See [15–17]. Many authors have researched the different delays Gause-type predator-prey systems [18–22].

For two-dimensional Gause-type predator-prey model,𝑑𝑥(𝑡)(𝑡)𝑔(𝑥(𝑡(𝑡)𝑝(𝑥(𝑡,(𝑡)𝑑𝑡=𝑥))−𝑦))𝑑𝑦(𝑡(𝑥(𝑡,𝑑𝑡=𝑦)−𝑑+𝜇𝑝)) (1.2) the time delay can be incorporated into (1.2) at three different terms [4], including the prey specific growth term 𝑔(𝑥(𝑡)), the predator response term 𝑝(𝑥(𝑡)), or the interaction term 𝑦(𝑡)𝑝(𝑥(𝑡)). The detailed analysis on stability and bifurcation can be found in [23]. Because the prey usually has a gestation or reaction time of predator, introducing time delay is necessary and hence changes the model into a DDE. It is well known that delay could cause a stable equilibrium to be unstable and induce bifurcations as well as periodic oscillations. Under the hypothesis that prey 𝑥(𝑡) has a gestation in (1.1), we modify it to be the following one:𝑑𝑥(𝑡)(𝑥)(),(𝑡)𝑑𝑡=𝑥𝑔−𝑦𝑝𝑥−𝜏𝑑𝑦(𝑥𝑦,(𝑡)𝑑𝑡=𝑦−ℎ+𝑒𝑝)−𝑧𝑞𝑑𝑧𝑦,𝑑𝑡=𝑧−𝑠+𝑚𝑞 (1.3) where 𝜏 is the time of gestation.

The purpose of current work is to analyze the effect of delay on the dynamics for (1.3), so we plan to employ bifurcation analysis approach with delay 𝜏 being the parameter. In particular, we choose 𝑔(𝑥)=𝛼(1−(𝑥/𝐾)), 𝑝(𝑥)=𝛽𝑥/(1+𝑝𝑥), and 𝑞(𝑦)=𝑟𝑦. So the following delay model will be obtained:𝑑𝑥(𝑡)𝑥𝑑𝑡=𝛼𝑥1−𝑘−𝛽𝑦𝑥(𝑡−𝜏)(),(𝑡)1+𝑝𝑥𝑡−𝜏𝑑𝑦𝑑𝑡=𝑦−ℎ+𝑒𝛽𝑥(𝑡)1+𝑝𝑥−𝑟𝑧𝑦,𝑑𝑧,𝑑𝑡=𝑧−𝑠+𝑚𝑟𝑦 (1.4) where 𝛼,𝛽,𝑘,𝑝,ℎ,𝑒,𝑟,𝑠,𝑚 are all positive parameters.

Our results reveal that Hopf bifurcation can occur as the delay crosses some critical values which leads to the existence of periodic solution that may conform to certain phenomena in ecosystem system.

The rest of the paper is organized as follows. In Section 2, we first investigate the stability of coexisting equilibrium and the existence of the Hopf bifurcation of (1.4) by analyzing the characteristic equation of the associated linearized system. In Section 3, we derive an explicit formula for determining the stability and the direction of bifurcating periodic solutions by the normal form method and the center manifold theory. In Section 4, we carry out some numerical simulations to illustrate the results obtained and forecast the change of population quantity.

2. Stability and Hopf Bifurcation of Coexisting Equilibrium

For the sake of convenience, we nondimensionalized (1.4) with the following: scaling𝑥𝑥→𝑘,𝑡→𝛼𝑡,𝑦→𝑦,𝑧→𝑧, (2.1) then (1.4) takes the form𝑑𝑥(𝑡)()𝑥()𝑑𝑡=𝑥1−𝑥−𝑎𝑦𝑡−𝜏𝑥(),(𝑡)𝑡−𝜏+𝑏𝑑𝑦𝑑𝑡=𝑦−𝑙+𝑐𝑥,(𝑡)𝑥+𝑏−𝑟𝑧𝑑𝑧,𝑑𝑡=𝑧−𝑠+𝑑𝑦 (2.2) where𝛽𝑎=1𝑝𝑘𝛼,𝑏=ℎ𝑝𝑘,𝑟=𝛼,𝑐=𝑒𝛽𝑝𝛼,𝑑=𝑚𝑟. (2.3)

Obviously, the delay cannot change the number of equilibria and non-dimensionalizations cannot change the properties of system. Through simple analysis, we know (2.2) has four equilibria: 𝐸₁(0,0,0), 𝐸₂(1,0,0), 𝐸₃(𝑥₀,𝑦₀,0), and 𝐸(x^∗,𝑦^∗,𝑧^∗) with𝑥₀=𝑏𝑙,𝑦𝑐−𝑙₀=1−𝑥₀𝑥₀+𝑏𝑎,𝑥^∗=(1−𝑏)+(1−𝑏)²+4𝑏−(4𝑎𝑠/𝑑)2,𝑦^∗=𝑠𝑑,𝑧^∗𝑙=−𝑟+𝑐𝑥^∗𝑟(𝑥^∗).+𝑏 (2.4)

There is no obvious biological significance for 𝐸₁(0,0,0) and 𝐸₂(1,0,0). In this paper, we mainly study the change of stability of coexisting equilibrium 𝐸(𝑥^∗,𝑦^∗,𝑧^∗) with the variation of time delay. If 𝑎𝑠/𝑑≤𝑏<1, then 𝐸(𝑥^∗,𝑦^∗,𝑧^∗) is the uniqueness equilibrium of (2.2). We consider the linearized system of (2.2) at 𝐸. The characteristic equation at 𝐸 is given by𝜆−𝑚₁₁−𝑛₁₁𝑒^−𝜆𝜏−𝑚₁₂0−𝑚₂₁𝜆−𝑚₂₃0−𝑚₃₂𝜆=0, (2.5) where𝑚₁₁=1−2𝑥^∗,𝑚₁₂=−𝑎𝑥^∗𝑥^∗+𝑏<0,𝑚₂₁=𝑏𝑐𝑦^∗(𝑥^∗+𝑏)²𝑚>0,₂₃=−𝑟𝑦^∗<0,𝑚₃₂=𝑑𝑧^∗>0,𝑛₁₁=−𝑎𝑏𝑦^∗(𝑥^∗+𝑏)²<0. (2.6) The eigenvalue 𝜆 satisfies the characteristic equation𝐷(𝜆,𝜏)=𝜆³+𝑎₂𝜆²+𝑎₁𝜆+𝑎₀+𝑏₂𝜆²+𝑏₀𝑒^−𝜆𝜏=0, (2.7) where 𝑎₂=−𝑚₁₁, 𝑎₁=−𝑚₁₂𝑚₂₁−𝑚₃₂𝑚₂₃>0, 𝑎₀=𝑚₃₂𝑚₂₃𝑚₁₁, 𝑏₂=−𝑛₁₁>0, and 𝑏₀=𝑚₃₂𝑚₂₃𝑛₁₁>0.

If 𝑚₁₁<0, then all the eigenvalues of (2.7) have negative real parts when 𝜏=0 by the Routh-Hurwitz criterion. So we have the Lemma 2.1.

Lemma 2.1. If 𝑎𝑠/𝑑≤𝑏<1 and 𝑚₁₁<0, then the coexisting equilibrium 𝐸(𝑥^∗,𝑦^∗,𝑧^∗) of (2.2) is locally asymptotically stable.

It is known that 𝐸(𝑥^∗,𝑦^∗,𝑧^∗) is asymptotically stable if all roots of the corresponding characteristic equation (2.7) have negative real parts. We shall study the distribution of the roots of the (2.7) when 𝜏≠0. We assume that the equilibrium 𝐸(𝑥^∗,𝑦^∗,𝑧^∗) of the ODE model is stable, then we derive some conditions to ensure that the steady state of the delay model is still stable.

We compute the eigenvalues of the Jacobian matrix at the coexisting equilibrium 𝐸(𝑥^∗,𝑦^∗,𝑧^∗). Substituting 𝜆=𝑖𝜔 into (2.7) yields,(1)when 𝜔=0, 𝐷(0,𝜏)=𝑎₀+𝑏₀=𝑚₂₃𝑚₃₂(𝑛₁₁+𝑚₁₁)≠0,(2)when 𝜔≠0, 𝐷(𝑖𝜔,𝜏)=(𝑖𝜔)³+𝑎₂(𝑖𝜔)²+𝑎₁𝑖𝜔+𝑎₀+(𝑏₀−𝑏₂𝜔²)𝑒^{−𝑖𝜔𝜏}=0. Separating the real and imaginary parts, we get−𝑎₂𝜔²+𝑎₀+𝑏₀−𝑏₂𝜔²cos𝜔𝜏=0,−𝜔³+𝑎₁𝜔−𝑏₀−𝑏₂𝜔²sin𝜔𝜏=0. (2.8) Consequently, we get𝜔⁶+𝑎²₂−2𝑎₁−𝑏²₂𝜔⁴+𝑎²₁−2𝑎₀𝑎₂+2𝑏₀𝑏₂𝜔²+𝑎²₀−𝑏²₀=0. (2.9) Let 𝜔²=𝑙, 𝐴=𝑎²₂−2𝑎₁−𝑏²₂, 𝐵=𝑎²₁−2𝑎₀𝑎₂+2𝑏₀𝑏₂ and 𝐶=𝑎²₀−𝑏²₀, then (2.9) becomes𝑙³+𝐴𝑙²+𝐵𝑙+𝐶=0. (2.10) From Ruan and Wei [24] or Li and Wei [25], we have the following results on the distribution of roots of (2.10).

Lemma 2.2. Denote 𝑙₁=−𝐴+√Δ3,𝑙₂=−𝐴−√Δ3,Δ=𝐴²−3𝐵,ℎ(𝑙)=𝑙³+𝐴𝑙²+𝐵𝑙+𝐶. (2.11)(1)If𝐶≥0, then (2.10) has at least one positive root.(2)If𝐶≥0 and Δ<0, then (2.10) has no positive root.(3)If𝐶≥0, then (2.10) has positive roots ⇔ if 𝑙₁>0 and ℎ′(𝑙₁)≤0.

Suppose (2.10) has root with positive real part. Without loss of generality, we assume that it has three positive roots, denoted by 𝑙₁,𝑙₂,𝑙₃, respectively. Then, (2.10) has three positive roots 𝜔_𝑖=𝑙_𝑖(𝑖=1,2,3). Let 𝜏_𝑘⁽⁰⁾ be the unique root of (2.8) such that 𝜏_𝑘⁽⁰⁾𝜔_𝑘∈[0,2𝜋). Also denote by𝜏_𝑘^(𝑗)=𝜏_𝑘⁽⁰⁾+2𝑗𝜋𝜔_𝑘, (2.12) for 𝑘=1,2,3.

So (±𝜔_𝑘,𝜏_𝑘^(𝑗)) is the solution of (2.7). Clearly,lim𝜏_𝑘^(𝑗)=∞,𝑘=1,2,3. (2.13) We can define𝜏₀=𝜏_𝑘0=min𝜏_𝑘⁽⁰⁾,𝜔₀=𝜔_𝑘0, (2.14) that is, ±𝑖𝜔₀ is the purely imaginary roots of (2.7) for 𝜏=𝜏₀. So we have the following.

Lemma 2.3. 𝜏₀ is defined by (2.14), then
(1) if one of the following holds (𝑖)𝐶<0; (𝑖𝑖)𝐶≥0, Δ≥0; (iii) 𝐶≥0,𝑙₁>0; ℎ′(𝑙₁)>0, then all roots of (2.7) have negative real parts for 𝜏∈[0,𝜏₀);
(2) if none of the conditions (i)–(𝑖𝑖𝑖) are satisfied, then all toots of (2.7) have negative real parts for all 𝜏≥0.

So we have the following theorem.

Theorem 2.4. Supposed that(a)𝑎𝑠/𝑑≤𝑏<1 and 𝑚₁₁<0,(b)either𝐶≥0,𝐴²−3𝐵<0, or𝐶≥0,𝐵>0,then the equilibrium 𝐸(𝑥^∗,𝑦^∗,𝑧^∗) of the delay model (2.2) is absolutely stable; that is, 𝐸(𝑥^∗,𝑦^∗,𝑧^∗) is asymptotically stable for all 𝜏≥0.

Denoting 𝜆(𝜏)=𝛼(𝜏)+𝛽(𝜏) to be the root of (2.7) satisfying 𝛼(𝜏₀)=0,𝜔(𝜏₀)=𝜔₀, then we have the following Lemma 2.5.

Lemma 2.5. Suppose 𝐴²−3𝐵<0,(a)±𝑖𝜔₀ is a pair of simple purely imaginary roots of (2.7) when 𝜏=𝜏₀,(b)[𝑑Re𝜆(𝜏)/𝑑𝜏]|_𝜏=𝜏0>0.

Proof. If 𝑖𝜔₀ is not simple, then 𝜔₀ must satisfy 𝑑𝜆𝑑𝜏³+𝑎₂𝜆²+𝑎₁𝜆+𝑎₀+𝑏₂𝜆²+𝑏₀𝑒^−𝜆𝜏|_{𝜆=𝑖𝜔0}=0, (2.15) that is, 𝜔₀ and 𝜏₀ must satisfy 𝜏𝑏₂𝜔²−𝑏₀cos𝜔𝜏+2𝑏₂𝜔sin𝜔𝜏=3𝜔²−𝑎₁,𝜏𝑏₂𝜔²−𝑏₀sin𝜔𝜏−2𝑏₂𝜔cos𝜔𝜏=2𝑎₂𝜔. (2.16) Eliminating 𝜏, we have 𝑎₁−3𝜔²sin𝜔𝜏+2𝑎₂𝜔cos𝜔𝜏=−2𝑏₂𝜔, (2.17) substituting sin𝜔𝜏, cos𝜔𝜏 from (2.8) into (2.17). We have 3𝜔⁴+2𝑎²₂−2𝑎₁−𝑏²₂𝜔²+𝑎²₁−2𝑎₂𝑎₀+2𝑏₂𝑏₀=0. (2.18) Recall 𝜔²₀=𝑙₀, which implies 3𝑙²₀+2𝐴𝑙₀+𝐵=0. (2.19) However if 𝐴²−3𝐵<0, we have ℎ′𝑙₀=3𝑙²₀+2𝐴𝑙₀+𝐵>0, (2.20) a contradiction. So ±𝑖𝜔₀ is a pair of simple purely imaginary roots of (2.7).
Note that 𝑀=𝑎₁−3𝜔²₀cos𝜔₀𝜏₀−2𝑎₂𝜔₀sin𝜔₀𝜏₀+𝑏₂𝜔²₀−𝑏₀𝜏₀,𝑁=𝑎₁−3𝜔²₀sin𝜔₀𝜏₀+2𝑎₂𝜔₀cos𝜔₀𝜏₀+2𝑏₂𝜔₀,𝑄=𝑏₀𝜔₀−𝑏₂𝜔³₀. (2.21) So 𝑑𝜆(𝜏)=𝑑𝜏𝑖𝑄=𝑀+𝑖𝑁𝑁𝑄+𝑖𝑀𝑄𝑀²+𝑁². (2.22) Then, 𝑑𝑅𝑒𝜆𝜏₀=𝑑𝜏𝑁𝑄𝑀²+𝑁² (2.23) has the same sign as that of 𝑁𝑄.
Now, 𝑁𝑄=𝑎₁−3𝜔²₀sin𝜔₀𝜏₀+2𝑎₂𝜔₀cos𝜔₀𝜏₀+2𝑏₂𝜔₀𝑏₀𝜔−𝑏₂𝜔³=𝜔²₀3𝜔⁴₀+2𝑎²₂−2𝑎₁−𝑏²₂𝜔²₀+𝑎²₁−2𝑎₂𝑎₀+2𝑏₂𝑏₀=𝜔²₀3𝜔⁴₀+2𝐴𝜔²₀+𝐵. (2.24) When 𝐴²−3𝐵<0, 𝑁𝑄 is greater than zero, that is [𝑑Re𝜆(𝜏)/𝑑𝜏]|_𝜏=𝜏0>0.

By Lemma 2.5, we have the Theorem 2.6.

Theorem 2.6. Supposed that Lemma 2.3 is right, then, when 𝐴²−3𝐵<0, the equilibrium 𝐸(𝑥^∗,𝑦^∗,𝑧^∗) of the delay model (2.2) is asymptotically stable when 𝜏<𝜏₀ and unstable when 𝜏>𝜏₀, where 𝜏₀ is defined by (2.14). When 𝜏=𝜏₀, a Hopf bifurcation occurs.

3. Direction and Stability of Hopf Bifurcation

Let 𝑥₁(𝑡)=𝑥(𝑡)−𝑥^∗, 𝑥₂(𝑡)=𝑦(𝑡)−𝑦^∗, 𝑥₃(𝑡)=𝑧(𝑡)−𝑧^∗, 𝑋_𝑖(𝑡)=𝑥_𝑖(𝜏𝑡) (𝑖=1,2,3), 𝜏=𝜏₀+𝜇, 𝜇∈ℝ, and𝑚𝐵=¹¹𝑚₁₂0𝑚₂₁0𝑚₂₃0𝑚₃₂0𝑛,𝐶=¹¹00000000. (3.1) The system (1.4) is transformed into a functional differential equation (FDE) in 𝐶=𝐶([−1,0],ℝ³), defining𝐿_𝜇𝜙=𝜏₀+𝜇𝐵𝜙(0)+𝜏₀+𝜇𝐶𝜙(−1), (3.2) where 𝜙=(𝜙₁,𝜙₂,𝜙₃)^𝑇∈𝐶([−1,0],ℝ³), and the nonlinear term isℎ𝜇,𝜙=𝜏₀+𝜇−2𝜙²₁(0)−2𝑎𝑏(𝑥^∗+𝑏)²𝜙₂(0)𝜙₁(−1)−2𝑏𝑐𝑦^∗(𝑏+𝑥^∗)³𝜙²₁(0)+2𝑏𝑐(𝑏+𝑥^∗)²𝜙₁(0)𝜙₂(0)−2𝑟𝜙₂(0)𝜙₃(0)2𝑑𝑟𝜙₂(0)𝜙₃(0). (3.3) So 𝜇=0 is Hopf bifurcation point.

By the Riesz representation theorem, there exists a 3×3 matrix 𝜂(𝜃,𝜇) (−1≤𝜃≤0), whose elements are of bounded variation functions such that𝐿_𝜇𝜙=⁰₋₁𝑑𝜂𝜃,𝜇𝜙(𝜃),for𝜙∈𝐶[−1,0],ℝ³. (3.4) In fact, we can choose𝜂=𝜏𝜃,𝜇⁰(),𝜏+𝜇𝐵,𝜃=0,0,𝜃∈−1,0₀+𝜇𝐶,𝜃=−1. (3.5) Then, (3.4) is satisfied. For 𝜙∈𝐶¹([−1,0],ℝ³), we define𝐴𝜇𝜙(𝜃)=𝑑𝜙(𝜃)[),𝑑𝜃,𝜃∈−1,0₀₋₁𝜙(𝜉)𝑅𝜇𝜙(𝜃)=[),ℎ𝑑𝜂𝜉,𝜇,𝜃=0,0,𝜃∈−1,0𝜇,𝜙,𝜃=0. (3.6) So (3.4) is equivalent to the following abstract equation:𝑥̇_𝑡=𝐴𝜇𝑥_𝑡+𝑅𝜇𝑥_𝑡, (3.7) where 𝑥=(𝑥₁,𝑥₂,𝑥₃)^𝑇, 𝑥_𝑡=𝑥(𝑡+𝜃), for 𝜃∈[−1,0].

For 𝜓∈𝐶¹([0,1],ℝ^3∗), we define𝐴^∗𝜓(𝑠)=−𝑑𝜓(𝑠)(],𝑑𝑠,𝑠∈0,1₀₋₁𝜓(−𝜉)𝑑𝜂(𝜉,0),𝑠=0, (3.8) and a bilinear form𝜓(𝑠),𝜙(𝜃)=𝜓(0)𝜙(0)−⁰₋₁^𝜃_𝜉=0𝜓(𝜉−𝜃)𝑑𝜂(𝜃)𝜙(𝜉)𝑑𝜉, (3.9) where 𝜂(𝜃)=𝜂(𝜃,0). Then, 𝐴(0) and 𝐴^∗ are adjoint operators. We know that ±𝑖𝜔₀𝜏₀ are eigenvalues of 𝐴(0), and, therefore, they are also eigenvalues of 𝐴^∗(0). The vector 𝑞(𝜃)=(𝑞₁,1,𝑞₃)^𝑇𝑒^{𝑖𝜔0𝜏0𝜃}(𝜃∈[−1,0]) and 𝑞^∗(𝑠)=𝐷(𝑞^∗₁,1,𝑞^∗₃)𝑒^{𝑖𝜔0𝜏0𝑠}(𝑠∈[0,1]) are the eigenvectors of 𝐴(0) and 𝐴^∗ corresponding to the eigenvalue 𝑖𝜔₀𝜏₀ and −𝑖𝜔₀𝜏₀, respectively, satisfying ⟨𝑞^∗(𝑠),𝑞(𝜃)⟩=1,⟨𝑞^∗(𝑠),𝑞(𝜃)⟩=0 with 𝑞₁=𝑚₁₂/(𝑖𝜔₀−𝑛₁₁𝑒^−𝑖𝜔0−𝑚₁₁),𝑞₃=𝑚₃₂/𝑖𝜔₀, 𝑞^∗₁=𝑚₂₁/(−𝑖𝜔₀−𝑛₁₁𝑒^𝑖𝜔0−𝑚₁₁),𝑞^∗₃=𝑚₂₃/𝑖𝜔₀, and 𝐷=1/(𝑞^∗₁𝑞₁+1+𝑞^∗₃𝑞₃+𝑞₁𝑞^∗₁𝑛₁₁𝑒^{−𝑖𝜔0𝜏0}). Following the same algorithms as Hassard et al. [26], we can obtain the coefficients which will be used to determine the important quantities:𝑔₀₂=2𝐷𝜏₀−𝑞₁²𝑞^∗₁−𝑎𝑏(𝑏+𝑥^∗)²𝑞₁𝑞^∗₁𝑒^{𝑖𝜔0𝜏0}−𝑏𝑐𝑦^∗(𝑏+𝑥^∗)²𝑞₁−𝑟𝑞₃+𝑑𝑞₃𝑞^∗₃,𝑔₂₀=2𝐷𝜏₀−𝑞²₁𝑞^∗₁−𝑎𝑏𝑞₁𝑞^∗₁(𝑏+𝑥^∗)²𝑒^{−𝑖𝜔0𝜏0}−𝑏𝑐𝑦^∗(𝑏+𝑥^∗)³𝑞²₁+𝑏𝑐(𝑏+𝑥^∗)²𝑞₁−𝑟𝑞₃+𝑑𝑞₃𝑞^∗₃,𝑔₁₁=2𝐷𝜏₀−2𝑞₁𝑞₁𝑞^∗₁−𝑎𝑏(𝑏+𝑥^∗)²2Re𝑞₁𝑒^{−𝑖𝜔0𝜏0}−2𝑏𝑐𝑦^∗(𝑏+𝑥^∗)³𝑞₁𝑞₁+2𝑏𝑐(𝑏+𝑥^∗)²Re𝑞₁−2𝑟Re𝑞₃+2𝑑Re𝑞₃𝑞^∗₃,𝑔₂₁=2𝐷𝜏₀−4𝑞₁𝑞^∗₁−4𝑏𝑐(𝑏+𝑥^∗)³𝑞₁+2𝑏𝑐(𝑏+𝑥^∗)²𝑊¹₁₁(0)+2𝑏𝑐(𝑏+𝑥^∗)²𝑞₁−2𝑎𝑏(𝑏+𝑥^∗)²𝑞₁𝑞^∗₁𝑒^{−𝑖𝜔0𝜏0}+2𝑑𝑞₃𝑞^∗₃−2𝑟𝑞₃𝑊²₁₁(0)+2𝑑𝑞^∗₃−2𝑟𝑊³₁₁(0)−2𝑎𝑏𝑞^∗₁(𝑏+𝑥^∗)²𝑊¹₁₁()+−1−4𝑞^∗₁𝑞₁−2𝑏𝑐(𝑏+𝑥^∗)³𝑞₁+𝑏𝑐(𝑏+𝑥^∗)²𝑊¹₂₀(0)+−𝑎𝑏𝑞^∗₁(𝑏+𝑥^∗)²𝑞₁𝑒^{𝑖𝜔0𝜏0}+𝑏𝑐(𝑏+𝑥^∗)²𝑞₁+𝑑𝑞^∗₃𝑞₃−𝑟𝑞₃𝑊²₂₀(0)+𝑑𝑞^∗₃−𝑟𝑊³₂₀(0)−𝑎𝑏𝑞^∗₁(𝑏+𝑥^∗)²𝑊¹₂₀(−1). (3.10) Since there are 𝑊₂₀(𝜃) and 𝑊₁₁(𝜃) in 𝑔₂₁, we still need to compute them. From [25], we have𝑊₂₀(𝜃)=𝑖𝑔₂₀𝜔₀𝜏₀𝑞(𝜃)+𝑖𝑔₀₂3𝜔₀𝜏₀𝑞(𝜃)+𝐸₁𝑒^{2𝑖𝜔0𝜏0𝜃}. (3.11) According to2𝑖𝜔₀𝜏₀𝐼−⁰₋₁𝑑𝜂(𝜃)𝑒^{2𝑖𝜔0𝜃}𝐸₁=ℎ_𝑧2, (3.12) whereℎ_𝑧2=−2+2𝑎𝑏𝑦^∗(𝑏+𝑥^∗)³𝑞²₁−2𝑎𝑏(𝑏+𝑥^∗)²𝑞₁−2𝑟𝑞₃+𝑏𝑐(𝑏+𝑥^∗)²𝑞₁𝑒^{−𝑖𝜔0𝜏0}−4𝑏𝑐𝑦^∗(𝑏+𝑥^∗)³𝑞²₁𝑒^{−2𝑖𝜔0𝜏0}2𝑑𝑞₃𝑒^{−𝑖𝜔0𝜏0}, (3.13) we have 𝐸₁=(𝐸₁⁽¹⁾,𝐸₁⁽²⁾,𝐸₁⁽³⁾)^𝑇, where 𝐸₁^(𝑖)=2Δ₁^(𝑖)/Δ₁(𝑖=1,2,3), withΔ₁=−4𝜔²₀𝜏²₀−𝑚₂₃𝑚₃₂−𝑚₁₂𝑚₂₁2𝑖𝜔₀𝜏₀+4𝜔²₀𝜏²₀+𝑚₂₃𝑚₃₂𝑛₁₁𝑒^{−2𝑖𝜔0𝜏0}+𝑚₁₁,Δ₁⁽¹⁾=4𝜔²₀𝜏²₀+𝑚₂₃𝑚₃₂2𝑞²₁+2𝑎𝑏(𝑏+𝑥^∗)²𝑞₁𝑒^{−𝑖𝜔0𝜏0}+2𝑚₁₂𝑖𝜔₀𝜏₀−2𝑏𝑐𝑦^∗(𝑏+𝑥^∗)³𝑞²₁+2𝑏𝑐(𝑏+𝑥^∗)²𝑞₁−2𝑟𝑞₃+2𝑚₁₂𝑚₂₃𝑑𝑞₃,Δ₁⁽²⁾=2𝑚₂₁𝑖𝜔₀𝜏₀−2𝑞²₁−2𝑎𝑏(𝑏+𝑥^∗)²𝑞₁𝑒^{−𝑖𝜔0𝜏0}+2𝑖𝜔₀𝜏₀2𝑖𝜔₀𝜏₀−𝑛₁₁𝑒^{−𝑖𝜔0𝜏0}−𝑚₁₁−2𝑏𝑐𝑦^∗(𝑏+𝑥^∗)³𝑞²₁+2𝑏𝑐(𝑏+𝑥^∗)²𝑞₁−2𝑟𝑞₃+2𝑚₂₃2𝑖𝜔₀𝜏₀−𝑛₁₁𝑒^{−2𝑖𝜔0𝜏0}−𝑚₁₁𝑑𝑞₃,Δ₁⁽³⁾=𝑚₂₁𝑚₃₂−2𝑞²₁−2𝑎𝑏(𝑏+𝑥^∗)²𝑞₁𝑒^{−𝑖𝜔0𝜏0}+𝑚₃₂2𝑖𝜔₀𝜏₀−𝑛₁₁𝑒^{−𝑖𝜔0𝜏0}−𝑚₁₁−2𝑏𝑐𝑦^∗(𝑏+𝑥^∗)³𝑞²₁+2𝑏𝑐(𝑏+𝑥^∗)²𝑞₁−2𝑟𝑞₃+22𝑖𝜔₀𝜏₀2𝑖𝜔₀𝜏₀−𝑛₁₁𝑒^{−𝑖𝜔0𝜏0}−𝑚₁₁−𝑚₁₂𝑚₂₁𝑑𝑞₃. (3.14) And, similarly,𝑊₁₁(𝜃)=−𝑖𝑔₁₁𝜔₀𝜏₀𝑞(𝜃)+𝑖𝑔₁₁𝜔₀𝜏₀𝑞(𝜃)+𝐸₂. (3.15) According to𝑑𝜂(⁰₋₁𝜃)𝐸₂=−ℎ_𝑧𝑧, (3.16) whereℎ_𝑧𝑧=−2+2𝑎𝑏𝑦^∗(𝑏+𝑥^∗)³𝑞₁𝑞₁−2𝑎𝑏(𝑏+𝑥^∗)²𝑞₁+𝑞₁−𝑟𝑞₃+𝑞₃+𝑏𝑐(𝑏+𝑥^∗)²𝑞₁𝑒^{𝑖𝜔0𝜏0}+𝑞₁𝑒^{−𝑖𝜔0𝜏0}−4𝑏𝑐𝑦^∗(𝑏+𝑥^∗)³𝑞₁𝑞₁𝑑𝑞₃𝑞^∗₃𝑒^{𝑖𝜔0𝜏0}+𝑞₃𝑞^∗₃𝑒^{−𝑖𝜔0𝜏0}, (3.17) we have 𝐸₂=(𝐸₂⁽¹⁾,𝐸₂⁽²⁾,𝐸₂⁽³⁾)^𝑇, where 𝐸₂^(𝑖)=2Δ₂^(𝑖)/Δ₂(𝑖=1,2,3), withΔ₂=𝑚₂₃𝑚₃₂𝑚₁₁+𝑛₁₁,Δ₂⁽¹⁾=−𝑚₂₃𝑚₃₂−4𝑞₁𝑞₁−2𝑎𝑏(𝑏+𝑥^∗)²2Re𝑞₁𝑒^{−𝑖𝜔0𝜏0}+2𝑚₁₂𝑚₂₃𝑑2Re𝑞₃,Δ₂⁽²⁾=−2𝑚₂₃𝑚₁₁+𝑛₁₁𝑑2Re𝑞₃,Δ₂⁽³⁾=𝑚₂₁𝑚₃₂−4𝑞₁𝑞₁−2𝑎𝑏(𝑏+𝑥^∗)²2Re𝑞₁𝑒^{−𝑖𝜔0𝜏0}−𝑚₃₂𝑚₁₁+𝑛₁₁−4𝑏𝑐𝑦^∗(𝑏+𝑥^∗)³𝑞₁𝑞₁+2𝑏𝑐(𝑏+𝑥^∗)²2Re𝑞₁−2𝑚₁₂𝑚₂₁𝑑2Re𝑞₃. (3.18)

Consequently, 𝑔_𝑖𝑗 can be expressed explicitly by the parameters and delay. Thus, we can compute the following values:𝑐₁(0)=𝑖2𝜔₀𝜏₀𝑔₁₁𝑔₂₀−2𝑔₁₁²−𝑔₀₂²3+𝑔₂₁2,𝜇₂𝑐=−𝑅𝑒₁(0)𝑅𝑒𝜆^′𝜏₀,𝑇₂=−𝐼𝑚𝑐₁(0)+𝜇₂𝐼𝑚𝜆^′𝜏₀𝜔₀𝜏₀,𝛽₂=2𝑅𝑒𝑐₁(0), (3.19) which determine the properties of bifurcating periodic solutions at the critical value 𝜏₀. That is, 𝜇₂ determines the direction of Hopf bifurcation: if 𝜇₂>0 (𝜇₂<0), then Hopf bifurcation at 𝜏₀ is forward (or backward); 𝛽₂ determines the stability of bifurcating periodic solutions: the bifurcating periodic solutions on the center manifold are stable (unstable) if 𝛽₂<0 (𝛽₂>0); 𝑇₂ determines the period of the bifurcating periodic solutions: the period increases (decreases) if 𝑇₂>0 (𝑇₂<0).

4. Numerical Simulations and Discussions

In this part, we perform some numerical simulations. The results not only support our previous parts but also predict the existence of global Hopf bifurcations periodic solutions. However, to prove this observation theoretically is of great challenge.

We choose the parameters as follows:

(1) 𝑎=0.321,𝑏=0.860,𝑐=0.532,𝑟=0.12,𝑠=0.312,𝑑=0.415,𝑙=0.13,

(2) 𝑎=0.432,𝑏=0.677,𝑐=0.321,𝑟=0.12,𝑠=0.235,𝑑=0.332,𝑙=0.013.

Thus, all the conditions in Theorem 2.4 are satisfied. 𝐸(0.855, 0.7519, 1.1272) with initial value (0.812, 0.702, 0.925) for (1) and 𝐸(0.7915, 0.7078, 1.3339) with initial value (0.607, 0.602, 0.895) for (2) are asymptotically stable for all 𝜏>0 (see Figure 1). From a biological sense, the prey, predator, and top predator will have a short-term fluctuation in the initial stage as the effect of 𝜏. But the population would tend to a steady level after a long period of time.

(a)

(b)

Figure 1 

𝐸 is asymptotically stable for all 𝜏>0.

We choose another set of parameters which satisfy the assumptions in Theorem 2.6:(3)𝜏=8.6017, 𝑎=0.036, 𝑠=0.0160, 𝑑=0.296, 𝑏=0.1347, 𝑟=0.9200, 𝑙=0.1500, 𝑐=0.6330.

Through (2.14), we have 𝜏₀=12.6774, 𝜔₀=14.7193, and 𝑐₁(0)=−0.4617+𝑖0.6024. Utilizing Theorem 2.6, we know that the equilibrium of the delay model (2.2) is asymptotically stable when 𝜏<𝜏₀ (see Figure 2).

(a)

(b)

(c)

(d)

Figure 2 

𝐸(0.9977, 0.6954, 0.4432) with initial value (0.985, 0.135, 0.4) is asymptotically stable when 𝜏=8.6017<𝜏0=12.6774.

Hopf bifurcation occurs when 𝜏=𝜏₀, and the bifurcating periodic solution is orbitally asymptotically for 𝜏>𝜏₀ (see Figure 3).

(a)

(b)

(c)

(d)

Figure 3 

A stable periodic orbit of system (2.2) when initial value is (0.985, 0.135, 0.4) and 𝜏=20.7288>𝜏0=12.6774.

In addition, the periodic solution of system (2.2) still exists when 𝜏 is large and its amplitude is larger compared with the solution in Figure 4. The numerical results of Figure 4 show the global existence of periodic solutions generated by the Hopf bifurcation. How to explain the phenomenon theoretically needs further researches.

(a)

(b)

Figure 4 

The stable periodic orbits of system (2.2) when 𝜏=200.7288 and 𝜏=400.8768 with parameters given by (3).

5. Conclusion

In this paper, we analyze the dynamics of the equilibria coexistence for a class of three-dimensional Gause-type predator-prey model. We obtain the stability of this equilibrium and also claim that the introduced delay changes its stability while a Hopf bifurcation occurs. The existence of the bifurcation periodic solutions for sufficiently large delay has been shown by numerical simulations.

Acknowledgments

The work is supported in part by the National Natural Science Foundation of China and by the Heilongjiang Provincial Natural Science Foundation (no. A200806, A200813).