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

Authors: Shuang Guo; Weihua Jiang
Journal: Journal of Applied Mathematics
 (2012)

---

## 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.


---

## Body

## 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)dx(t)dt=xg(x)-yp(x),dy(t)dt=y[-h+ep(x)]-zq(y),dz(t)dt=z[-s+dq(y)],
where x(t), y(t), and z(t) are the population densities of prey, predator, and top predator at time t, respectively. g(x) is the intrinsic growth rate of prey; p(x) and q(y) are the specific growth rates of predator and top predator; h, s>0 are the death rates of y(t) and z(t); e, d>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 thatg(x), p(x), and q(y) satisfy the following conditions.(H1)
g(0)>0, there exists K>0 such that g(K)=0 (K is called the carrying capacity of prey species), and g′(x)≤0 for 0≤x<K.(H2)
p(0)=0 and p′(x)>0.(H3)
q(0)=0 and q′(y)>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 g(x)=α(1-x/K), p(x)=βx/(1+px), and q(y)=ry. For simplicity, we nondimensionalize the system (1.1) with the following scaling:
(1.2)x→xK,t→αt,y→y,z→z.
Then the system (1.1) takes the form(1.3)dx(t)dt=x(1-x)-axyx+b,dy(t)dt=y(-l+cxx+b-rz),dz(t)dt=z(-s+dy),
which satisfies x(0)=α1>0, y(0)=α2>0, and z(0)=α3>0, where(1.4)a=βpkα,b=1pk,r=hα,c=eβpα,d=mr.Through simple analysis, we know that (1.3) has four equilibriums: E1(0,0,0), E2(1,0,0), E3(x0,y0,0), and E¯(x*,y*,z*), where
(1.5)x0=blc-l,y0=(1-x0)(x0+b)a,x*=(1-b)+(1-b)2+4b-4as/d2,y*=sd,z*=-lr+cx*r(x*+b).
There is no obvious biological significance for E1(0,0,0) and E2(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 g(x(t)), a time delay τ in the predator response term p(x(t)), a time delay τ in the interaction term y(t)p(x(t)). 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 p(x), that is,(1.6)dx(t)dt=x(1-x)-ayx(t-τ)x(t-τ)+b,dy(t)dt=y(-l+cx(t-τ)x(t-τ)+b-rz),dz(t)dt=z(-s+dy).
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 Section2, we construct a Lyapunov function to prove the global stability of the equilibrium E3(x0,y0,0) indicating the extinction of top-predator. In Section 3, we first investigate the stability of coexisting equilibrium E¯(x*,y*,z*) 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 E3. The characteristic equation at the equilibrium E3 is given by(2.1)D(λ,τ)=λ3-(m33+m11+n11)λ2+(m11m33+n11m33-m12n21)λ+m12m33n21=0,
where
(2.2)m11=1-2x0,m12=-ax0x0+b<0,m23=-ry0<0,m33=-s+dy0,n11=-aby0(x0+b)2<0,n21=bcy0(x0+b)2>0.Whenm11<0 and m33≤0, E3 is locally asymptotically stable by the Routh-Hurwitz criterion, so we have the following.Lemma 2.1.
Ifbl/(c-l)>1/2 and  (c-(b+1)l)/(c-l)2≤as/bcd,  E3(x0,y0,0) is locally asymptotically stable.To show thatE3(x0,y0,0) is globally stable, we construct a Lyapunov function:
(2.3)V(x,y,z)=∫x0xcp(ξ)-lp(ξ)dξ+a∫y0yη-y0ηdη+radz.It is obvious thatV(x,y,z)∈C1(ℝ+3,ℝ), V(x0,y0,0)=0, and V(x,y,z)>0 for (x,y,z)∈ℝ+3-{(x0,y0,0)}. For simplicity, let g(x)=1-x, and p(x)=x/(b+x); then x′=0 if and only if y=xg(x)/ap(x); y′=0 if and only if x=x0, where x0 satisfies p(x0)=l/c:(2.4)V̇(x,y,z)|(1.3)=cp(x)-lp(x)[xg(x)-ap(x)y]+a(y-y0)(-l+cp(x)-rz)+radz(-s+dy)=cp(x)-lp(x)[xg(x)-ap(x)y+ap(x)y0-ap(x)y0]+a(y-y0)[-l+cp(x)]+a(y-y0)(-rz)+radz(-s+dy)=cp(x)-lp(x)[xg(x)-ap(x)y0+ap(x)(y0-y)]+a(y-y0)[-l+cp(x)]+a(y-y0)(-rz)+radz(-s+dy)=(cp(x)-l)[xg(x)p(x)-ay0]+razy0-rasdz=ca(p(x)-p(x0))[xg(x)p(x)-y0]+raz(y0-sd)≤0.We havey0-s/d≤0 when m33≤0, and it follows that(2.5)V̇(x,y,z)|(1.3)≤0,for0<x<1,y>0,z>0.
Hence, we have the following theorem by Lasalle’s invariance principle [13].Theorem 2.2.
Ifbl/(c-l)>1/2 and (c-(b+1)l)/(c-l)2≤as/bcd, E3(x0,y0,0) is global asymptotically stable.


## 3. Stability and Hopf Bifurcation of Coexisting Equilibrium

In this part, we mainly study the stability of the coexisting equilibriumE¯(x*,y*,z*). If as/d≤b<1, then E¯(x*,y*,z*) is a unique nontrivial equilibrium of (1.6). We consider the linearized system of (1.6) at E¯. The equation of Jacobian determinant at the equilibrium E¯ is given by(3.1)|λ-m11-n11e-λτ-m120-n21e-λτλ-m230-m32λ|=0,
where(3.2)m11=1-2x*-aby*(x*+b)2,m12=-ax*x*+b<0,m23=-ry*<0,m32=dz*>0,n11=-aby*(x*+b)2<0,n21=bcy*(x*+b)2>0.
The eigenvalue λ satisfies the characteristic equation(3.3)D(λ,τ)=λ3+a2λ2+a1λ+a0+(b2λ2+b1λ+b0)e-λτ=0,
where a2=-m11, a1=-m32m23>0, a0=-m32m23m11, b2=-n11>0, b1=-m12m21>0, and b0=m32m23n11>0.Ifm11<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.
Ifas/d≤b<1 and  m11<0, then the coexisting equilibrium  E¯(x*,y*,z*) of (1.3) is locally asymptotically stable.It is known thatE¯(x*,y*,z*) 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 E¯(x*,y*,z*) 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λ=iω into (3.3):(1) whenω=0, D(0,τ)=a0+b0=m23m32(m11+n11)≠0;(2) whenω≠0, D(iω,τ)=(iω)3+a2(iω)2+a1iω+a0+(-b2ω2+b1iω+b0)e-iωτ=0. Separating the real and imaginary parts gives(3.4)-a2ω2+a0-b2ω2cos⁡ωτ+b1ωsin⁡ωτ+b0cos⁡ωτ=0,-ω3+a1ω+b2ω2sin⁡ωτ+b1ωcos⁡ωτ-b0sin⁡ωτ=0.
We get(3.5)ω6+(a22-2a1-b22)ω4+(a12-2a0a2+2b0b2-b12)ω2+a02-b02=0.
Let ω2=l, P1=a22-2a1-b22, P2=a12-2a0a2+2b0b2-b12, and P3=a02-b02; then (3.5) becomes(3.6)l3+P1l2+P2l+P3=0.
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)
P3≥0,P12-3P2<0,(b)
P3≥0,P2>0,P1>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)
as/d≤b<1 and m11<0,(b)
eitherP3≥0,P12-3P2<0 or P3≥0,P2>0,P1>0.Then the equilibrium E¯(x*,y*,z*) of the delay model (1.6) is absolutely stable; that is, E¯(x*,y*,z*) 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 Lemma3.2 are not satisfied.Lemma 3.4.
Denote(3.7)h(l)=l3+P1l2+P2l+P3,l0=-P1+P12-3P23.
If one of the followings holds:(a)
P3<0,(b)
P3≥0,P2<0,(c)
P3≥0,l0>0,h(l0)≤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 l1,l2, and l3, respectively. Then (3.6) has three positive roots ωi=li(i=1,2,3). Let τk(0) be the unique root of (3.4) such that τk(0)ωk∈[0,2π). Also denote (3.8)τk(j)=τk(0)+2jπωk,
for k=1,2,3.So(±ωk,τk(j)) is the solution of (3.3). Clearly,(3.9)lim⁡τk(j)=∞,k=1,2,3.
We can define(3.10)τ0=τk0=min⁡τk(0),ω0=ωk0;
that is, ±iω0 are the purely imaginary roots of (3.3) for τ=τ0.Denotingλ(τ)=α(τ)+β(τ) be the root of (3.3) satisfying α(τ0)=0,ω(τ0)=ω0, we have the following lemma.Lemma 3.5.
Ifh′(ω02)≠0, one has  (dReλ(τ)/dτ)|τ=τ0>0.Proof.
Assume that(3.11)h(ω02)=ω06+(a22-2a1-b22)ω04+(a12-2a0a2+2b0b2-b12)ω02+a02-b02;
then
(3.12)h′(ω02)=3ω04+2(a22-2a1-b22)ω02+(a12-2a0a2+2b0b2-b12).
Because
(3.13)sign⁡[dReλ(τ)dτ]|τ=τ0=sign⁡[dReλ(τ)dτ]-1|τ=τ0,[dλ(τ)dτ]-1=(3λ2+2a2λ+a1)+(2b2λ+b1)e-λτ-τ(b2λ2+b1λ+b0)e-λτλ(b2λ2+b1λ+b0)e-λτ=(3λ2+2a2λ+a1)eλτλ(b2λ2+b1λ+b0)+2b2λ+b1λ(b2λ2+b1λ+b0)-τλ,
we have
(3.14)sign⁡[dReλ(τ)dτ]-1|τ=τ0=sign⁡[Re[(3λ2+2a2λ+a1)eλτλ(b2λ2+b1λ+b0)]τ=τ0+Re[2b2λ+b1λ(b2λ2+b1λ+b0)]|τ=τ0]=sign⁡[3ω04+2(a22-2a1-b22)ω02+(a12-2a0a2+2b0b2-b12)b12ω02+(b0-b2ω02)2]=sign⁡[h′(ω02)b12ω02+(b0-b2ω02)2]=sign⁡(h′(ω02)).
Ifh′(ω02)≠0, then (dReλ(τ)/dτ)|τ=τ0≠0. There must be (dReλ(τ)/dτ)|τ=τ0>0. This is because (3.3) has the positive real part roots as τ<τ0 if (dReλ(τ)/dτ)|τ=τ0<0. This contradicts to the fact when τ∈[0,τ0) and E¯(x*,y*,z*) is asymptotically stable.By Lemma3.5 we have the following theorem.Theorem 3.6.
Suppose thatas/d≤b<1 and  m11<0. If Lemma 3.4 holds, then the equilibrium  E¯(x*,y*,z*) of the delay model (1.6) is asymptotically stable when  τ<τ0, and unstable when  τ>τ0, where  τ0 is defined by (3.10). In addition, if  h′(ω02)≠0, then Hopf bifurcation occurs when  τ=τ0.


## 4. Direction and Stability of Hopf Bifurcation

Letx1(t)=x(t)-x*, x2(t)=y(t)-y*, x3(t)=z(t)-z*, Xi(t)=xi(τt)(i=1,2,3), τ=τ0+μ, μ∈ℝ, and(4.1)B=(m11-1m1200m230m320),C=(n1100n2100000).
The system (1.6) is transformed into a functional differential equation (FDE) in C=C([-1,0],ℝ3)[28], defining(4.2)Lμ(ϕ)=(τ0+μ)Bϕ(0)+(τ0+μ)Cϕ(-1),
where ϕ=(ϕ1,ϕ2,ϕ3)T∈C([-1,0],ℝ3). And the nonlinear term is(4.3)h(μ,ϕ)=(τ0+μ)(-2ϕ12(0)-2bn11x*+bϕ12(-1)-2ϕ2(0)ϕ3(0)-n21y*ϕ2(0)ϕ1(-1)+2n21x*+bϕ12(-1)2drϕ2(0)ϕ3(0)).
Obviously μ=0 is a Hopf bifurcation point. So the system (1.6) can transform into an abstract functional differential equation:(4.4)v(t)̇=Lμ(vt)+h(μ,vt),
where v(t)=(x1(t),x2(t),x3(t))T∈ℝ3.There exists a3×3 matrix η(θ,μ) (-1≤θ≤0), whose elements are of bounded variation functions such that (4.5)Lμ(ϕ)=∫-10[dη(θ,μ)]ϕ(θ),forϕ∈C([-1,0],R3).
In fact, we can choose(4.6)η((θ),μ)={(τ0+μ)B,θ=0,0,θ∈(-1,0),(τ0+μ)C,θ=-1.
Then (4.5) is satisfied. For ϕ∈C1([-1,0],ℝ3), we define(4.7)A(μ)ϕ(θ)={dϕ(θ)dθ,θ∈[-1,0),∫-10[dη(ξ,μ)]ϕ(ξ),θ=0,R(μ)ϕ(θ)={0,θ∈[-1,0),h(μ,ϕ),θ=0.
So (4.5) is equivalent to the following abstract equation:(4.8)ẋt=A(μ)xt+R(μ)xt,
where x=(x1,x2,x3)T and xt=x(t+θ) for θ∈[-1,0].Forψ∈C1([0,1],ℝ3*), we define(4.9)A*ψ(s)={-dψ(s)ds,s∈(0,1],∫-10ψ(-ξ)dη(ξ,0),s=0,
and a bilinear form:(4.10)〈ψ(s),ϕ(θ)〉=ψ̅(0)ϕ(0)-∫-10∫ξ=0θψ̅(ξ-θ)dη(θ)ϕ(ξ)dξ,
where η(θ)=η(θ,0). Then A(0) and A* are adjoint operators. We know that ±iω0τ0 are eigenvalues of A(0) and therefore they are also eigenvalues of A*(0). The vectors q(θ)=(1,α,β)Teiω0τ0θ(θ∈[-1,0]) and q*(s)=D(α*,β*,1)eiω0τ0s(s∈[0,1]) are the eigenvectors of A(0) and A* corresponding to the eigenvalue iω0τ0 and -iω0τ0, respectively, satisfying 〈q*(s),q(θ)〉=1,〈q*(s),q̅(θ)〉=0 with α=(iω0-m11)/m12, β=(iω0(iω0)-m11)/m23m12, α*=((iω02)-m23m32)/m12m23, β*=-iω0/m23, and D̅=1/(α̅*+αβ̅*+β-(α̅*n11+β̅*n21)τ0eiω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: (4.11)g02=2D̅τ0[-2α̅*-2rβ̅*β+2dα̅β̅+2(β̅*n21-α̅*bn11)x*+be2iω0-n21β̅*α̅y*eiω0],g20=2D̅τ0[-2α̅*-2rβ̅*β+2dαβ+2(β̅*n21-α̅*bn11)x*+be-2iω0-n21β̅*αy*e-iω0],g11=2D̅τ0[-2α̅*-2rβ̅*Reβ+2dReαβ̅+2(β̅*n21-α̅*bn11)x*+b-n21β̅*α̅y*Reαeiω0],g21=2D̅τ0[(2β*n21eiω0x*+b-2α̅*n11eiω0x*+b-2α̅β̅*n21y*)-4α̅*W111(0)+(-2α̅*-2rβ̅*β̅)W201(0)+(2dα-2rβ̅*)W113(0)+2dα̅W203(0)+(2β*n21e-iω0x*+b-2α̅*n11e-iω0x*+b-2αβ̅*n21y*)W111(-1)+(2β*n21eiω0x*+b-2α̅*n11eiω0x*+b-2α̅β̅*n21y*)W201(-1)].
Since there are W20(θ) and W11(θ) in g21, we still need to compute them. From [29], we have(4.12)W20(θ)=ig20ω0τ0q(θ)+ig̅023ω0τ0q̅(θ)+E1e2iω0τ0θ.
According to(4.13)[2iω0τ0I-∫-10dη(θ)e2iω0θ]E1=hz2,
where(4.14)hz2=(-2α̅*-2rβ̅*β+2α̅β̅d2(β̅*n21-α̅*bn11)x*+be-2iω0-n21β̅*α̅y*e-iω00),
we have E1=(E1(1),E1(2),E1(3))T, where E1(i)=2Δ1(i)/Δ1(i=1,2,3) with(4.15)Δ1=(2iω0-m11)(-4ω02-m23m32),Δ1(1)=(8ω02+2m32m23)(α̅*+rβ̅*β-dα̅β̅)+2m12iω0(2(β̅*n21-α̅*bn11)x*+b)e-2iω0-n21β̅*α̅y*e-iω0,Δ1(2)=(-4ω02-2m11iω0)(2(β̅*n21-α̅*bn11)x*+b)e-2iω0-n21β̅*α̅y*e-iω0,Δ1(3)=(-m11m32+2m32iω0)(2(β̅*n21-α̅*bn11)x*+b)e-2iω0-n21β̅*α̅y*e-iω0.
And similarly,(4.16)W11(θ)=-ig11ω0τ0q(θ)+ig̅11ω0τ0q̅(θ)+E2.
According to(4.17)(∫-10dη(θ))E2=-hzz̅,
where(4.18)hzz̅=(-2α*-2rβ̅*Re{β}+2dRe{αβ̅}+2(β̅*n21-α̅*bn11)x*+b-n21β̅*y*Re{αeiω0}0),
we have E2=(E2(1),E2(2),E2(3))T, where E2(i)=2Δ2(i)/Δ2(i=1,2,3) with(4.19)Δ2=m23m32(m11+n11),Δ2(1)=-2m23m32(α*-rβ̅*Re{β}+2dRe{αβ̅}+β̅*n21-α̅*bn11x*+b),Δ2(2)=0,Δ2(3)=-2n21m32(α*+rβ̅*Re{β}-dRe{αβ̅}+2dRe{αβ̅}+β̅*n21-α̅*bn11x*+b)+m32(m11+n11)n21β̅*y*Re{αeiω0}.Consequently,gij can be expressed explicitly by the parameters and delay in the system (4.2). Thus, we can compute the following values:(4.20)c1(0)=i2ω0τ0(g11g20-2|g11|2-|g02|23)+g212,μ2=-Re(c1(0))Re(λ′(τ0)),T2=-Im⁡c1(0)+μ2Im⁡λ′(τ0)ω0τ0,β2=2Re(c1(0)),
which determine the properties of bifurcating periodic solutions at the critical value τ0. That is, μ2 determines the direction of Hopf bifurcation: if μ2>0(μ2<0), then Hopf bifurcation at τ0 is forward (or backward); β2 determines the stability of bifurcating periodic solutions: β2<0(β2>0); the bifurcating periodic solution is orbitally asymptotically stable (unstable); T2 determines the period of the bifurcating periodic solutions: the period increases (decreases) if T2>0(T2<0).


## 5. Numerical Simulations

In order to check our conclusions, we perform some numerical simulations. We choose the parameters as follows:(1)
a=0.3232, b=0.4017, c=0.399, r=0.282, s=0.235, d=0.303, and l=0.015;(2)
a=0.432, b=0.6700, c=0.32, r=0.12, s=0.235, d=0.333, and l=0.031.Thus, all the conditions in Theorem3.3 are satisfied. E¯(0.7013,0.7755,1.2312) with initial value (0.48, 0.32, 0.25) for (1) and E¯(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.E̅ is asymptotically stable for all τ>0.
(a)(b)We choose another set of parameters which satisfies Theorem3.6.(3)τ=8.6208, a=0.6000, s=0.3000, d=0.6000, b=0.3700, r=0.2200, l=0.1500, and c=0.3300.By (3.10), we have τ0=15.7288, ω0=4.6133, and c1(0)=-0.0643+i0.0031. Theorem 3.6 indicates that the equilibrium of the delay model (1.6) is asymptotically stable when τ<τ0 (see Figure 2).E¯(0.9433,0.5000,0.2871) with initial value (0.985, 0.34, 0.4) is asymptotically stable when τ=8.6208<τ0=15.7288.
(a)(b)(c)(d)Hopf bifurcation occurs whenτ=τ0, and the bifurcating periodic solution is orbitally asymptotically for τ>τ0 (see 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.
(a)(b)(c)(d)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.The stable periodic orbits of system (1.6) when τ=420.7288 and τ=890.7288 with parameters given by (3).
(a)(b)


## 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×1010 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.



---

*Source: 260798-2012-03-20.xml*