Global Stability and Hopf Bifurcation for Gause-Type Predator-Prey System

Abstract

A class of three-dimensional Gause-type predator-prey model is considered. Firstly, local stability of equilibrium indicating the extinction of top-predator is obtained. Meanwhile, we construct a Lyapunov function, which is an extension of the Lyapunov functions constructed by Hsu for predator-prey system (2005), to give the global stability of the equilibrium. Secondly, we analyze the stability of coexisting equilibrium of predator-prey system with time delay when the predator catches the prey of pregnancy or with growth time. The delay can lead to periodic solutions, which is consistent with the law of growth for birds and some mammals. Further, an explicit formula is given which determines the stability of the bifurcating periodic solutions theoretically and the existence of periodic solutions is displayed by numerical simulations.

1. Introduction

The predator-prey systems have been extensively studied. The most popular one is Lotka-Volterra type model, which exhibits the well-known “paradox of enrichment” observed by [1, 2] and so on. However, there is often the interaction among multiple populations in nature, whose relationships are more complex than those in two populations. So the dynamics of three-dimensional model may become more complicated [3–8].

In general, the relationship among three groups may be competitive, a predator and two preys, two species catching the same prey, or a food chain. See [6, 9–11]. The Gause-type predator-prey food chain model was proposed by Freedman and Waltman in 1977 [5], which can be described as follows:𝑑𝑥(𝑡)(𝑥)(𝑥),(𝑡)𝑑𝑡=𝑥𝑔−𝑦𝑝𝑑𝑦(𝑥𝑦,(𝑡)𝑑𝑡=𝑦−ℎ+𝑒𝑝)−𝑧𝑞𝑑𝑧𝑦,𝑑𝑡=𝑧−𝑠+𝑑𝑞 (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. For (1.1), there have been a lot of results, including the properties of equilibria and bifurcations. Freedman has argued that the unique interior equilibrium exists and it is locally asymptotically stable [5]. In [7, 12], the authors performed the normal forms for Hopf bifurcations and saddle node bifurcations near a degenerate equilibrium and showed that the model within a certain parameter range could take on chaos by numerical simulations. Many biologists believe that if the unique equilibrium of a predator-prey system is locally asymptotically stable, then it is globally asymptotically stable. For proving the global stability of the equilibrium, one can construct a Lyapunov function. If we are able to construct a Lyapunov function for the system, then the global stability is directly established from the modified LaSalle’s invariant principle [13]. Many authors presented some globally qualitative analysis of solutions of system (1.1) (see [8, 14–17]).

We assume that 𝑔(𝑥), 𝑝(𝑥), and 𝑞(𝑦) satisfy the following conditions.(H1)𝑔(0)>0, there exists 𝐾>0 such that 𝑔(𝐾)=0 (𝐾 is called the carrying capacity of prey species), and 𝑔^′(𝑥)≤0 for 0≤𝑥<𝐾.(H2)𝑝(0)=0 and 𝑝^′(𝑥)>0.(H3)𝑞(0)=0 and 𝑞^′(𝑦)>0.

In general, there are three prototypes of monotone response functions which we often refer to as Holling I, II, and III, respectively. In [9, 15], Chiu and Hsu considered the three-dimensional food chain model that the response function is of Holling II; in [18], Hsu et al. obtained a complete classification about the asymptotic behavior of the solutions of Gause-type predator-prey system with Holling II function. In this paper, we discuss the general Gause-type food chain models, in particular, the case of Holling’s type I and II functional responses for predator and prey, respectively. So we let 𝑔(𝑥)=𝛼(1−𝑥/𝐾), 𝑝(𝑥)=𝛽𝑥/(1+𝑝𝑥), and 𝑞(𝑦)=𝑟𝑦. For simplicity, we nondimensionalize the system (1.1) with the following scaling: 𝑥𝑥→𝐾,𝑡→𝛼𝑡,𝑦→𝑦,𝑧→𝑧. (1.2) Then the system (1.1) takes the form𝑑𝑥(𝑡)()−𝑑𝑡=𝑥1−𝑥𝑎𝑥𝑦,(𝑡)𝑥+𝑏𝑑𝑦𝑑𝑡=𝑦−𝑙+𝑐𝑥,(𝑡)𝑥+𝑏−𝑟𝑧𝑑𝑧,𝑑𝑡=𝑧−𝑠+𝑑𝑦 (1.3) which satisfies 𝑥(0)=𝛼₁>0, 𝑦(0)=𝛼₂>0, and 𝑧(0)=𝛼₃>0, where𝛽𝑎=1𝑝𝑘𝛼,𝑏=ℎ𝑝𝑘,𝑟=𝛼,𝑐=𝑒𝛽𝑝𝛼,𝑑=𝑚𝑟. (1.4)

Through simple analysis, we know that (1.3) has four equilibriums: 𝐸₁(0,0,0), 𝐸₂(1,0,0), 𝐸₃(𝑥₀,𝑦₀,0), and 𝐸(𝑥^∗,𝑦^∗,𝑧^∗), where 𝑥₀=𝑏𝑙,𝑦𝑐−𝑙₀=1−𝑥₀𝑥₀+𝑏𝑎,𝑥^∗=(1−𝑏)+(1−𝑏)²+4𝑏−4𝑎𝑠/𝑑2,𝑦^∗=𝑠𝑑,𝑧^∗𝑙=−𝑟+𝑐𝑥^∗𝑟(𝑥^∗).+𝑏 (1.5) There is no obvious biological significance for 𝐸₁(0,0,0) and 𝐸₂(1,0,0). The behavior of the solutions of (1.3) can be very complicated (see [12, 19–21]). We know that the time delay can be incorporated into (1.1) in three different ways: a time delay 𝜏 in the prey-specific growth term 𝑔(𝑥(𝑡)), a time delay 𝜏 in the predator response term 𝑝(𝑥(𝑡)), a time delay 𝜏 in the interaction term 𝑦(𝑡)𝑝(𝑥(𝑡)). The predators have a gestation period or reaction time such as grass-hare-fox food chain. Thus, a time delay 𝜏 can be incorporated into the predator response term 𝑝(𝑥), that is,𝑑𝑥(𝑡)()𝑥()𝑑𝑡=𝑥1−𝑥−𝑎𝑦𝑡−𝜏𝑥(),(𝑡)𝑡−𝜏+𝑏𝑑𝑦()𝑑𝑡=𝑦−𝑙+𝑐𝑥𝑡−𝜏𝑥(),(𝑡)𝑡−𝜏+𝑏−𝑟𝑧𝑑𝑧.𝑑𝑡=𝑧−𝑠+𝑑𝑦 (1.6) In general, delay differential equations (DDEs) exhibit much more complicated dynamics than ordinary differential equations (ODEs) such as the existence of Bogdanov-Takens bifurcation and even chaos (see [21, 22], etc.). The reason is that a time delay could destabilize a stable equilibrium and induce bifurcations (see [23–26]). The cyclic phenomenon is very significant for the stability and the balance of ecosystems. Using the delay as a bifurcation parameter, we investigate the stability of the positive equilibrium and the existence of Hopf bifurcation of the model with delay. It is shown that Hopf bifurcation can occur as the delay crosses some critical values.

The rest of our paper is organized as follows: in Section 2, we construct a Lyapunov function to prove the global stability of the equilibrium 𝐸₃(𝑥₀,𝑦₀,0) indicating the extinction of top-predator. In Section 3, we first investigate the stability of coexisting equilibrium 𝐸(𝑥^∗,𝑦^∗,𝑧^∗) and the existence of Hopf bifurcation of (1.6). To determine the dynamics of the delay model, we study the characteristic equation of the linearized system at the equilibrium. In Section 4, we derive an explicit formula for determining the stability of bifurcating periodic solutions and the direction of Hopf bifurcation by the normal form method and the center manifold theory. In Section 5, we carry out some numerical simulations to illustrate the results obtained.

2. Global Stability of the Extinction of Top Predator

We consider the linearized system of (1.3) at 𝐸₃. The characteristic equation at the equilibrium 𝐸₃ is given by𝐷(𝜆,𝜏)=𝜆³−𝑚₃₃+𝑚₁₁+𝑛₁₁𝜆²+𝑚₁₁𝑚₃₃+𝑛₁₁𝑚₃₃−𝑚₁₂𝑛₂₁𝜆+𝑚₁₂𝑚₃₃𝑛₂₁=0, (2.1) where 𝑚₁₁=1−2𝑥₀,𝑚₁₂=−𝑎𝑥₀𝑥₀+𝑏<0,𝑚₂₃=−𝑟𝑦₀𝑚<0,₃₃=−𝑠+𝑑𝑦₀,𝑛₁₁=−𝑎𝑏𝑦₀𝑥₀+𝑏²<0,𝑛₂₁=𝑏𝑐𝑦₀𝑥₀+𝑏²>0. (2.2)

When 𝑚₁₁<0 and 𝑚₃₃≤0, 𝐸₃ is locally asymptotically stable by the Routh-Hurwitz criterion, so we have the following.

Lemma 2.1. If 𝑏𝑙/(𝑐−𝑙)>1/2 and (𝑐−(𝑏+1)𝑙)/(𝑐−𝑙)²≤𝑎𝑠/𝑏𝑐𝑑, 𝐸₃(𝑥₀,𝑦₀,0) is locally asymptotically stable.

To show that 𝐸₃(𝑥₀,𝑦₀,0) is globally stable, we construct a Lyapunov function: 𝑉𝑥,𝑦,𝑧=^𝑥_𝑥0𝑐𝑝(𝜉)−𝑙𝑝(𝜉)𝑑𝜉+𝑎^𝑦_𝑦0𝜂−𝑦₀𝜂𝑑𝜂+𝑟𝑎𝑑𝑧. (2.3)

It is obvious that 𝑉(𝑥,𝑦,𝑧)∈𝐶¹(ℝ³₊,ℝ), 𝑉(𝑥₀,𝑦₀,0)=0, and 𝑉(𝑥,𝑦,𝑧)>0 for (𝑥,𝑦,𝑧)∈ℝ³₊−{(𝑥₀,𝑦₀,0)}. For simplicity, let 𝑔(𝑥)=1−𝑥, and 𝑝(𝑥)=𝑥/(𝑏+𝑥); then 𝑥′=0 if and only if 𝑦=𝑥𝑔(𝑥)/𝑎𝑝(𝑥); 𝑦′=0 if and only if 𝑥=𝑥₀, where 𝑥₀ satisfies 𝑝(𝑥₀)=𝑙/𝑐:𝑉̇𝑥,𝑦,𝑧_(1.3)=𝑐𝑝(𝑥)−𝑙𝑝(𝑥)𝑥𝑔(𝑥)−𝑎𝑝(𝑥)𝑦+𝑎𝑦−𝑦₀−𝑙+𝑐𝑝(𝑥)−𝑟𝑧+𝑟𝑎𝑑𝑧=(𝑥)−𝑠+𝑑𝑦𝑐𝑝−𝑙𝑝(𝑥)𝑥𝑔(𝑥)−𝑎𝑝(𝑥)𝑦+𝑎𝑝(𝑥)𝑦₀−𝑎𝑝(𝑥)𝑦₀+𝑎𝑦−𝑦₀−𝑙+𝑐𝑝(𝑥)+𝑎𝑦−𝑦₀(−𝑟𝑧)+𝑟𝑎𝑑𝑧=(𝑥)−𝑠+𝑑𝑦𝑐𝑝−𝑙𝑝(𝑥)𝑥𝑔(𝑥)−𝑎𝑝(𝑥)𝑦₀+𝑎𝑝(𝑥)𝑦₀−𝑦+𝑎𝑦−𝑦₀−𝑙+𝑐𝑝(𝑥)+𝑎𝑦−𝑦₀(−𝑟𝑧)+𝑟𝑎𝑑𝑧=(𝑥)(𝑥)−𝑠+𝑑𝑦𝑐𝑝−𝑙𝑥𝑔𝑝(𝑥)−𝑎𝑦₀+𝑟𝑎𝑧𝑦₀−𝑟𝑎𝑠𝑑𝑧=𝑐𝑎𝑝(𝑥)−𝑝𝑥₀𝑥𝑔(𝑥)𝑝(𝑥)−𝑦₀+𝑟𝑎𝑧𝑦₀−𝑠𝑑≤0. (2.4)

We have 𝑦₀−𝑠/𝑑≤0 when 𝑚₃₃≤0, and it follows that𝑉̇𝑥,𝑦,𝑧_(1.3)≤0,for0<𝑥<1,𝑦>0,𝑧>0. (2.5) Hence, we have the following theorem by Lasalle’s invariance principle [13].

Theorem 2.2. If 𝑏𝑙/(𝑐−𝑙)>1/2 and (𝑐−(𝑏+1)𝑙)/(𝑐−𝑙)²≤𝑎𝑠/𝑏𝑐𝑑, 𝐸₃(𝑥₀,𝑦₀,0) is global asymptotically stable.

3. Stability and Hopf Bifurcation of Coexisting Equilibrium

In this part, we mainly study the stability of the coexisting equilibrium 𝐸(𝑥^∗,𝑦^∗,𝑧^∗). If 𝑎𝑠/𝑑≤𝑏<1, then 𝐸(𝑥^∗,𝑦^∗,𝑧^∗) is a unique nontrivial equilibrium of (1.6). We consider the linearized system of (1.6) at 𝐸. The equation of Jacobian determinant at the equilibrium 𝐸 is given by𝜆−𝑚₁₁−𝑛₁₁𝑒^−𝜆𝜏−𝑚₁₂0−𝑛₂₁𝑒^−𝜆𝜏𝜆−𝑚₂₃0−𝑚₃₂𝜆=0, (3.1) where𝑚₁₁=1−2𝑥^∗−𝑎𝑏𝑦^∗(𝑥^∗+𝑏)²,𝑚₁₂=−𝑎𝑥^∗𝑥^∗+𝑏<0,𝑚₂₃=−𝑟𝑦^∗𝑚<0,₃₂=𝑑𝑧^∗>0,𝑛₁₁=−𝑎𝑏𝑦^∗(𝑥^∗+𝑏)²<0,𝑛₂₁=𝑏𝑐𝑦^∗(𝑥^∗+𝑏)²>0. (3.2) The eigenvalue 𝜆 satisfies the characteristic equation𝐷(𝜆,𝜏)=𝜆³+𝑎₂𝜆²+𝑎₁𝜆+𝑎₀+𝑏₂𝜆²+𝑏₁𝜆+𝑏₀𝑒^−𝜆𝜏=0, (3.3) where 𝑎₂=−𝑚₁₁, 𝑎₁=−𝑚₃₂𝑚₂₃>0, 𝑎₀=−𝑚₃₂𝑚₂₃𝑚₁₁, 𝑏₂=−𝑛₁₁>0, 𝑏₁=−𝑚₁₂𝑚₂₁>0, and 𝑏₀=𝑚₃₂𝑚₂₃𝑛₁₁>0.

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

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

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

Now we substitute 𝜆=𝑖𝜔 into (3.3):

(1) when 𝜔=0, 𝐷(0,𝜏)=𝑎₀+𝑏₀=𝑚₂₃𝑚₃₂(𝑚₁₁+𝑛₁₁)≠0;

(2) when 𝜔≠0, 𝐷(𝑖𝜔,𝜏)=(𝑖𝜔)³+𝑎₂(𝑖𝜔)²+𝑎₁𝑖𝜔+𝑎₀+(−𝑏₂𝜔²+𝑏₁𝑖𝜔+𝑏₀)𝑒^{−𝑖𝜔𝜏}=0. Separating the real and imaginary parts gives−𝑎₂𝜔²+𝑎₀−𝑏₂𝜔²cos𝜔𝜏+𝑏₁𝜔sin𝜔𝜏+𝑏₀cos𝜔𝜏=0,−𝜔³+𝑎₁𝜔+𝑏₂𝜔²sin𝜔𝜏+𝑏₁𝜔cos𝜔𝜏−𝑏₀sin𝜔𝜏=0. (3.4) We get𝜔⁶+𝑎²₂−2𝑎₁−𝑏²₂𝜔⁴+𝑎²₁−2𝑎₀𝑎₂+2𝑏₀𝑏₂−𝑏²₁𝜔²+𝑎²₀−𝑏²₀=0. (3.5) Let 𝜔²=𝑙, 𝑃₁=𝑎²₂−2𝑎₁−𝑏²₂, 𝑃₂=𝑎²₁−2𝑎₀𝑎₂+2𝑏₀𝑏₂−𝑏²₁, and 𝑃₃=𝑎²₀−𝑏²₀; then (3.5) becomes𝑙³+𝑃₁𝑙²+𝑃₂𝑙+𝑃₃=0. (3.6) From Ruan and Wei [27], we have the following results on the distribution of roots of (3.6).

Lemma 3.2. If one of the followings holds: (a)𝑃₃≥0,𝑃²₁−3𝑃₂<0,(b)𝑃₃≥0,𝑃₂>0,𝑃₁>0,then (3.6) has no positive root.

Therefore, the real parts of all the eigenvalues of (3.3) are negative for all delay 𝜏≥0, so we have the following.

Theorem 3.3. Suppose that (a)𝑎𝑠/𝑑≤𝑏<1 and 𝑚₁₁<0,(b)either 𝑃₃≥0,𝑃²₁−3𝑃₂<0 or 𝑃₃≥0,𝑃₂>0,𝑃₁>0.Then the equilibrium 𝐸(𝑥^∗,𝑦^∗,𝑧^∗) of the delay model (1.6) is absolutely stable; that is, 𝐸(𝑥^∗,𝑦^∗,𝑧^∗) is asymptotically stable for all 𝜏≥0.

However, the stability of the steady state depends on the delay and the delay could even induce oscillations if the conditions in Lemma 3.2 are not satisfied.

Lemma 3.4. Denote ℎ(𝑙)=𝑙³+𝑃₁𝑙²+𝑃₂𝑙+𝑃₃,𝑙₀=−𝑃₁+𝑃₁²−3𝑃₂3. (3.7) If one of the followings holds:(a)𝑃₃<0,(b)𝑃₃≥0,𝑃₂<0,(c)𝑃₃≥0,𝑙₀>0,ℎ(𝑙₀)≤0,then (3.6) has at least one positive root. This implies that the characteristic equation (3.3) has at least a pair of purely imaginary roots.

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

So (±𝜔_𝑘,𝜏_𝑘^(𝑗)) is the solution of (3.3). Clearly,lim𝜏_𝑘^(𝑗)=∞,𝑘=1,2,3. (3.9) We can define𝜏₀=𝜏_𝑘0=min𝜏_𝑘⁽⁰⁾,𝜔₀=𝜔_𝑘0; (3.10) that is, ±𝑖𝜔₀ are the purely imaginary roots of (3.3) for 𝜏=𝜏₀.

Denoting 𝜆(𝜏)=𝛼(𝜏)+𝛽(𝜏) be the root of (3.3) satisfying 𝛼(𝜏₀)=0,𝜔(𝜏₀)=𝜔₀, we have the following lemma.

Lemma 3.5. If ℎ^′(𝜔²₀)≠0, one has (𝑑Re𝜆(𝜏)/𝑑𝜏)|_𝜏=𝜏0>0.

Proof. Assume that ℎ𝜔²₀=𝜔₀⁶+𝑎²₂−2𝑎₁−𝑏²₂𝜔₀⁴+𝑎²₁−2𝑎₀𝑎₂+2𝑏₀𝑏₂−𝑏²₁𝜔₀²+𝑎²₀−𝑏²₀; (3.11) then ℎ^′𝜔²₀=3𝜔⁴₀+2𝑎²₂−2𝑎₁−𝑏²₂𝜔²₀+𝑎²₁−2𝑎₀𝑎₂+2𝑏₀𝑏₂−𝑏²₁. (3.12) Because (𝜏)sign𝑑Re𝜆𝑑𝜏_𝜏=𝜏0(𝜏)=sign𝑑Re𝜆𝑑𝜏⁻¹_𝜏=𝜏0,𝑑𝜆(𝜏)𝑑𝜏⁻¹=3𝜆²+2𝑎₂𝜆+𝑎₁+2𝑏₂𝜆+𝑏₁𝑒^−𝜆𝜏−𝜏𝑏₂𝜆²+𝑏₁𝜆+𝑏₀𝑒^−𝜆𝜏𝜆𝑏₂𝜆²+𝑏₁𝜆+𝑏₀𝑒^−𝜆𝜏=3𝜆²+2𝑎₂𝜆+𝑎₁𝑒^𝜆𝜏𝜆𝑏₂𝜆²+𝑏₁𝜆+𝑏₀+2𝑏₂𝜆+𝑏₁𝜆𝑏₂𝜆²+𝑏₁𝜆+𝑏₀−𝜏𝜆, (3.13) we have (𝜏)sign𝑑Re𝜆𝑑𝜏⁻¹_𝜏=𝜏0=signRe3𝜆²+2𝑎₂𝜆+𝑎₁𝑒𝜆𝜏𝜆𝑏₂𝜆²+𝑏₁𝜆+𝑏_0𝜏=𝜏0+Re2𝑏₂𝜆+𝑏₁𝜆𝑏₂𝜆²+𝑏₁𝜆+𝑏_0𝜏=𝜏0=sign3𝜔⁴₀+2𝑎²₂−2𝑎₁−𝑏²₂𝜔²₀+𝑎²₁−2𝑎₀𝑎₂+2𝑏₀𝑏₂−𝑏²₁𝑏²₁𝜔²₀+𝑏₀−𝑏₂𝜔²₀²ℎ=sign^′𝜔²₀𝑏²₁𝜔²₀+𝑏₀−𝑏₂𝜔²₀²=signℎ^′𝜔²₀. (3.14)
If ℎ′(𝜔²₀)≠0, then (𝑑Re𝜆(𝜏)/𝑑𝜏)|_𝜏=𝜏0≠0. There must be (𝑑Re𝜆(𝜏)/𝑑𝜏)|_𝜏=𝜏0>0. This is because (3.3) has the positive real part roots as 𝜏<𝜏₀ if (𝑑Re𝜆(𝜏)/𝑑𝜏)|_𝜏=𝜏0<0. This contradicts to the fact when 𝜏∈[0,𝜏₀) and 𝐸(𝑥^∗,𝑦^∗,𝑧^∗) is asymptotically stable.

By Lemma 3.5 we have the following theorem.

Theorem 3.6. Suppose that 𝑎𝑠/𝑑≤𝑏<1 and 𝑚₁₁<0. If Lemma 3.4 holds, then the equilibrium 𝐸(𝑥^∗,𝑦^∗,𝑧^∗) of the delay model (1.6) is asymptotically stable when 𝜏<𝜏₀, and unstable when 𝜏>𝜏₀, where 𝜏₀ is defined by (3.10). In addition, if ℎ^′(𝜔²₀)≠0, then Hopf bifurcation occurs when 𝜏=𝜏₀.

4. Direction and Stability of Hopf Bifurcation

Let 𝑥₁(𝑡)=𝑥(𝑡)−𝑥^∗, 𝑥₂(𝑡)=𝑦(𝑡)−𝑦^∗, 𝑥₃(𝑡)=𝑧(𝑡)−𝑧^∗, 𝑋_𝑖(𝑡)=𝑥_𝑖(𝜏𝑡)(𝑖=1,2,3), 𝜏=𝜏₀+𝜇, 𝜇∈ℝ, and𝑚𝐵=¹¹−1𝑚₁₂00𝑚₂₃0𝑚₃₂0𝑛,𝐶=𝑛¹¹00₂₁00000. (4.1) The system (1.6) is transformed into a functional differential equation (FDE) in 𝐶=𝐶([−1,0],ℝ³)[28], defining𝐿_𝜇𝜙=𝜏₀+𝜇𝐵𝜙(0)+𝜏₀+𝜇𝐶𝜙(−1), (4.2) where 𝜙=(𝜙₁,𝜙₂,𝜙₃)^𝑇∈𝐶([−1,0],ℝ³). And the nonlinear term isℎ𝜇,𝜙=𝜏₀+𝜇−2𝜙²₁(0)−2𝑏𝑛₁₁𝑥^∗𝜙+𝑏²¹(−1)−2𝜙₂(0)𝜙₃(0)−𝑛₂₁𝑦^∗𝜙₂(0)𝜙₁(−1)+2𝑛₂₁𝑥^∗𝜙+𝑏²¹(−1)2𝑑𝑟𝜙₂(0)𝜙₃(0). (4.3) Obviously 𝜇=0 is a Hopf bifurcation point. So the system (1.6) can transform into an abstract functional differential equation:𝑣(𝑡)̇=𝐿_𝜇𝑣_𝑡+ℎ𝜇,𝑣_𝑡, (4.4) where 𝑣(𝑡)=(𝑥₁(𝑡),𝑥₂(𝑡),𝑥₃(𝑡))^𝑇∈ℝ³.

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

For 𝜓∈𝐶¹([0,1],ℝ^3∗), we define𝐴^∗𝜓(𝑠)=−𝑑𝜓(𝑠)(],𝑑𝑠,𝑠∈0,1₀₋₁𝜓(−𝜉)𝑑𝜂(𝜉,0),𝑠=0, (4.9) and a bilinear form:𝜓(𝑠),𝜙(𝜃)=𝜓(0)𝜙(0)−⁰₋₁^𝜃_𝜉=0𝜓(𝜉−𝜃)𝑑𝜂(𝜃)𝜙(𝜉)𝑑𝜉, (4.10) where 𝜂(𝜃)=𝜂(𝜃,0). Then 𝐴(0) and 𝐴^∗ are adjoint operators. We know that ±𝑖𝜔₀𝜏₀ are eigenvalues of 𝐴(0) and therefore they are also eigenvalues of 𝐴^∗(0). The vectors 𝑞(𝜃)=(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 𝛼=(𝑖𝜔₀−𝑚₁₁)/𝑚₁₂, 𝛽=(𝑖𝜔₀(𝑖𝜔₀)−𝑚₁₁)/𝑚₂₃𝑚₁₂, 𝛼^∗=((𝑖𝜔²₀)−𝑚₂₃𝑚₃₂)/𝑚₁₂𝑚₂₃, 𝛽^∗=−𝑖𝜔₀/𝑚₂₃, and 𝐷=1/(𝛼^∗+𝛼𝛽^∗+𝛽−(𝛼^∗𝑛₁₁+𝛽^∗𝑛₂₁)𝜏₀𝑒^{𝑖𝜔0𝜏0}). Following the same algorithms as Hassard et al. [29], we can obtain the coefficients which will be used to determine the important quantities: 𝑔₀₂=2𝐷𝜏₀−2𝛼^∗−2𝑟𝛽^∗𝛽+2𝑑𝛼2𝛽+𝛽^∗𝑛₂₁−𝛼^∗𝑏𝑛₁₁𝑥^∗𝑒+𝑏^2𝑖𝜔0−𝑛₂₁𝛽^∗𝛼𝑦^∗𝑒^𝑖𝜔0,𝑔₂₀=2𝐷𝜏₀−2𝛼^∗−2𝑟𝛽^∗2𝛽+2𝑑𝛼𝛽+𝛽^∗𝑛₂₁−𝛼^∗𝑏𝑛₁₁𝑥^∗𝑒+𝑏^{−2𝑖𝜔0}−𝑛₂₁𝛽^∗𝛼𝑦^∗𝑒^−𝑖𝜔0,𝑔₁₁=2𝐷𝜏₀−2𝛼^∗−2𝑟𝛽^∗Re𝛽+2𝑑Re𝛼2𝛽+𝛽^∗𝑛₂₁−𝛼^∗𝑏𝑛₁₁𝑥^∗−𝑛+𝑏²¹𝛽^∗𝛼𝑦^∗Re𝛼𝑒^𝑖𝜔0,𝑔₂₁=2𝐷𝜏₀−4𝛼^∗𝑊¹₁₁(0)+−2𝛼^∗−2𝑟𝛽^∗𝛽𝑊¹₂₀(0)+2𝑑𝛼−2𝑟𝛽^∗𝑊³₁₁(0)+2𝑑𝛼𝑊³₂₀(0)+2𝛽^∗𝑛₂₁𝑒^−𝑖𝜔0𝑥^∗−2+𝑏𝛼^∗𝑛₁₁𝑒^−𝑖𝜔0𝑥^∗−+𝑏2𝛼𝛽^∗𝑛₂₁𝑦^∗𝑊¹₁₁()+−12𝛽^∗𝑛₂₁𝑒^𝑖𝜔0𝑥^∗−2+𝑏𝛼^∗𝑛₁₁𝑒^𝑖𝜔0𝑥^∗−2+𝑏𝛼𝛽^∗𝑛₂₁𝑦^∗𝑊¹₂₀(−1). (4.11) Since there are 𝑊₂₀(𝜃) and 𝑊₁₁(𝜃) in 𝑔₂₁, we still need to compute them. From [29], we have𝑊₂₀(𝜃)=𝑖𝑔₂₀𝜔₀𝜏₀𝑞(𝜃)+𝑖𝑔₀₂3𝜔₀𝜏₀𝑞(𝜃)+𝐸₁𝑒^{2𝑖𝜔0𝜏0𝜃}. (4.12) According to2𝑖𝜔₀𝜏₀𝐼−⁰₋₁𝑑𝜂(𝜃)𝑒^{2𝑖𝜔0𝜃}𝐸₁=ℎ_𝑧2, (4.13) whereℎ_𝑧2=−2𝛼^∗−2𝑟𝛽^∗𝛽+2𝛼2𝛽𝑑𝛽^∗𝑛₂₁−𝛼^∗𝑏𝑛₁₁𝑥^∗𝑒+𝑏^{−2𝑖𝜔0}−𝑛₂₁𝛽^∗𝛼𝑦^∗𝑒^−𝑖𝜔00, (4.14) we have 𝐸₁=(𝐸₁⁽¹⁾,𝐸₁⁽²⁾,𝐸₁⁽³⁾)^𝑇, where 𝐸₁^(𝑖)=2Δ₁^(𝑖)/Δ₁(𝑖=1,2,3) withΔ₁=2𝑖𝜔₀−𝑚₁₁−4𝜔²₀−𝑚₂₃𝑚₃₂,Δ₁⁽¹⁾=8𝜔²₀+2𝑚₃₂𝑚₂₃𝛼^∗+𝑟𝛽^∗𝛽−𝑑𝛼𝛽+2𝑚₁₂𝑖𝜔₀2𝛽^∗𝑛₂₁−𝛼^∗𝑏𝑛₁₁𝑥^∗𝑒+𝑏^{−2𝑖𝜔0}−𝑛₂₁𝛽^∗𝛼𝑦^∗𝑒^−𝑖𝜔0,Δ₁⁽²⁾=−4𝜔²₀−2𝑚₁₁𝑖𝜔₀2𝛽^∗𝑛₂₁−𝛼^∗𝑏𝑛₁₁𝑥^∗𝑒+𝑏^{−2𝑖𝜔0}−𝑛₂₁𝛽^∗𝛼𝑦^∗𝑒^−𝑖𝜔0,Δ₁⁽³⁾=−𝑚₁₁𝑚₃₂+2𝑚₃₂𝑖𝜔₀2𝛽^∗𝑛₂₁−𝛼^∗𝑏𝑛₁₁𝑥^∗𝑒+𝑏^{−2𝑖𝜔0}−𝑛₂₁𝛽^∗𝛼𝑦^∗𝑒^−𝑖𝜔0. (4.15) And similarly,𝑊₁₁(𝜃)=−𝑖𝑔₁₁𝜔₀𝜏₀𝑞(𝜃)+𝑖𝑔₁₁𝜔₀𝜏₀𝑞(𝜃)+𝐸₂. (4.16) According to𝑑𝜂(⁰₋₁𝜃)𝐸₂=−ℎ_𝑧𝑧, (4.17) whereℎ_𝑧𝑧=−2𝛼^∗−2𝑟𝛽^∗Re𝛽+2𝑑Re𝛼𝛽+2𝛽^∗𝑛₂₁−𝛼^∗𝑏𝑛₁₁𝑥^∗−𝑛+𝑏₂₁𝛽^∗𝑦^∗Re𝛼𝑒^𝑖𝜔00, (4.18) we have 𝐸₂=(𝐸₂⁽¹⁾,𝐸₂⁽²⁾,𝐸₂⁽³⁾)^𝑇, where 𝐸₂^(𝑖)=2Δ₂^(𝑖)/Δ₂(𝑖=1,2,3) withΔ₂=𝑚₂₃𝑚₃₂𝑚₁₁+𝑛₁₁,Δ₂⁽¹⁾=−2𝑚₂₃𝑚₃₂𝛼^∗−𝑟𝛽^∗Re𝛽+2𝑑Re𝛼𝛽+𝛽^∗𝑛₂₁−𝛼^∗𝑏𝑛₁₁𝑥^∗,Δ+𝑏₂₍₂₎Δ=0,₂₍₃₎=−2𝑛₂₁𝑚₃₂𝛼^∗+𝑟𝛽^∗Re𝛽−𝑑Re𝛼𝛽+2𝑑Re𝛼𝛽+𝛽^∗𝑛₂₁−𝛼^∗𝑏𝑛₁₁𝑥^∗+𝑏+𝑚₃₂𝑚₁₁+𝑛₁₁𝑛₂₁𝛽^∗𝑦^∗Re𝛼𝑒^𝑖𝜔0. (4.19)

Consequently, 𝑔_𝑖𝑗 can be expressed explicitly by the parameters and delay in the system (4.2). Thus, we can compute the following values:𝑐₁(0)=𝑖2𝜔₀𝜏₀𝑔₁₁𝑔₂₀−2𝑔₁₁²−𝑔₀₂²3+𝑔₂₁2,𝜇₂𝑐=−Re₁(0)Re𝜆^′𝜏₀,𝑇₂=−Im𝑐₁(0)+𝜇₂Im𝜆′𝜏₀𝜔₀𝜏₀,𝛽₂=2Re𝑐₁(0), (4.20) 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: 𝛽₂<0(𝛽₂>0); the bifurcating periodic solution is orbitally asymptotically stable (unstable); 𝑇₂ determines the period of the bifurcating periodic solutions: the period increases (decreases) if 𝑇₂>0(𝑇₂<0).

5. Numerical Simulations

In order to check our conclusions, we perform some numerical simulations. We choose the parameters as follows:(1)𝑎=0.3232, 𝑏=0.4017, 𝑐=0.399, 𝑟=0.282, 𝑠=0.235, 𝑑=0.303, and 𝑙=0.015;(2)𝑎=0.432, 𝑏=0.6700, 𝑐=0.32, 𝑟=0.12, 𝑠=0.235, 𝑑=0.333, and 𝑙=0.031.

Thus, all the conditions in Theorem 3.3 are satisfied. 𝐸(0.7013,0.7755,1.2312) with initial value (0.48, 0.32, 0.25) for (1) and 𝐸(0.9037,0.7057,0.8022)with initial value (0.08, 0.02, 0.25) 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 shock 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 satisfies Theorem 3.6.

(3) 𝜏=8.6208, 𝑎=0.6000, 𝑠=0.3000, 𝑑=0.6000, 𝑏=0.3700, 𝑟=0.2200, 𝑙=0.1500, and 𝑐=0.3300.

By (3.10), we have 𝜏₀=15.7288, 𝜔₀=4.6133, and 𝑐₁(0)=−0.0643+𝑖0.0031. Theorem 3.6 indicates that the equilibrium of the delay model (1.6) is asymptotically stable when 𝜏<𝜏₀ (see Figure 2).

(a)

(b)

(c)

(d)

Figure 2 

𝐸(0.9433,0.5000,0.2871) with initial value (0.985, 0.34, 0.4) is asymptotically stable when 𝜏=8.6208<𝜏0=15.7288.

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 (1.6) when initial value is (0.985, 0.34, 0.4) and 𝜏=28.7206>𝜏0=15.7288.

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

(a)

(b)

Figure 4 

The stable periodic orbits of system (1.6) when 𝜏=420.7288 and 𝜏=890.7288 with parameters given by (3).

6. Conclusion

In this paper, we analyze the dynamics of the equilibria of the extinction of top-predator and coexistence for a class of three-dimensional Gause-type predator-prey model. We obtain that the equilibrium of the extinction of top-predator is not only locally but also globally asymptotically stable for the certain parameters; introducing delay changes the stability of the coexisting equilibrium and Hopf bifurcation occurs with the increase of 𝜏. The existence of periodic solutions for sufficiently large delay has been shown by numerical simulations. We know that the population outbreak may happen for the species with periodic fluctuation. The outbreaks of pests and mice are famous. For example, the number of locusts is estimated as many as 1.6×10¹⁰ and their weight is 5000 t, which appeared in a plague of locusts in Somalia in 1957 [30]. Thus, it is of great significance to research multiple periodic solutions of biological systems for controlling insect pests, preventing epidemics, and maintaining ecological balance.

Acknowledgments

The work is supported in part by the National Natural Science Foundation of China, by the Heilongjiang Provincial Natural Science Foundation (nos. A200806, A200813), and by the Heilongjiang Provincial Young Backbone Project (11251G0011).