# Stability and Bifurcation Analysis for a Predator-Prey Model with Discrete and Distributed Delay

Authors: Ruiqing Shi; Junmei Qi; Sanyi Tang
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

We propose a two-dimensional predatory-prey model with discrete and distributed delay. By the use of a new variable, the original two-dimensional system transforms into an equivalent three-dimensional system. Firstly, we study the existence and local stability of equilibria of the new system. And, by choosing the time delayτ as a bifurcation parameter, we show that Hopf bifurcation can occur as the time delay τ passes through some critical values. Secondly, by the use of normal form theory and central manifold argument, we establish the direction and stability of Hopf bifurcation. At last, an example with numerical simulations is provided to verify the theoretical results. In addition, some simple discussion is also presented.


---

## Body

## 1. Introduction

Since the pioneering theoretical works by Lotka [1] and Volterra [2], there were a lot of authors who studied all kinds of predator-prey models modeled by ordinary differential equations (ODEs). To reflect that the dynamical behavior of the models depends on the past history of the system, it is often necessary to incorporate time delays into the models. Therefore, a more realistic predator-prey model should be described by delayed differential equations (DDEs) [3–11]. In general, delay differential equations exhibit more complicated dynamics on stability, periodic structure, bifurcation, and so on [12–26]. In [27, 28], the authors investigated the effect of the discrete delay on the stability of the model. In [29], the effect of the distributed delay on the stability of the model was investigated. In [11], the authors proposed a Logistic model with discrete and distributed delays:
(1)x′(t)=rx(t)[1-a1x(t-τ)-a2∫-∞t‍f(t-s)x(s)ds],
where the parameters r, τ, a1, a2 are positive constants. The function f in (1) is called the delayed kernel, which is the weight given to the population t time units ago. And it was assumed that f(t)≥0 for all t≥0, together with the normalization condition
(2)∫0∞‍f(t)dt=1,
which ensures that the steady states of the model (1) are unaffected by the delay. They studied the stability of the positive equilibrium and existence of Hopf bifurcations, and direction and stability of the Hopf bifurcation were also analyzed. In [7], the authors proposed and investigated the following predator-prey model with time delay:
(3)x′(t)=x(t)[r1-a11x(t-τ)-a12y(t-τ)],y′(t)=y(t)[-r2+a21x(t-τ)-a22y(t-τ)],
where x(t) and y(t) can be interpreted as the population densities of the prey and the predator at time t, respectively. r1>0 denotes the intrinsic growth rate of the prey, and r2>0 denotes the death rate of the predator. For the convenience of computation, they chose the same τ>0 as delays; the delay τ represents the feedback time delay of the prey species to the growth of itself in term a11x(t-τ), represents the feedback time delay of the predator species to the growth of itself in term a22y(t-τ), represents the hunting delay in term a12y(t-τ), and represents the time of the predator maturation in term a21x(t-τ). The parameters  aij(i,j=1,2) are all positive constants. They studied the stability of the positive equilibrium and existence of Hopf bifurcations.Motivated by [7, 11, 27–29] and the references cited therein, in the present paper, we will consider the following predator-prey model with discrete and distributed delay:
(4)N′(t)=N(t)[α1-a11N(t-τ)-a12P(t-τ)],P′(t)=P(t)[-α2+a21∫-∞t‍G(t-s)N(s)ds-a22P(t-τ)],
where  N(t) and  P(t) can be interpreted as the population densities of the prey and the predator at time t, respectively. α1>0 denotes the intrinsic growth rate of the prey and α2>0 denotes the death rate of the predator; τ>0 represents the feedback time delay of the prey species to the growth of itself in term a11N(t-τ), represents the feedback time delay of the predator species to the growth of itself in term a22P(t-τ), and represents the hunting delay in term a12P(t-τ); aij>0(i,j=1,2). The function G(t) is the same as the function f(s) in system (1). Following the ideas of Cushing [30], we define G(t) as the following weak kernel function:
(5)G(t)=αe-αt,α>0.Next, we define a new variable:(6)u(t)=∫-∞t‍αe-α(t-s)N(s)ds,
then using the linear chain trick technique, system (4) can be transformed into the following equivalent system:
(7)N′(t)=N(t)[α1-a11N(t-τ)-a12P(t-τ)],P′(t)=P(t)[-α2+a21u(t)-a22P(t-τ)],u′(t)=αN(t)-αu(t).The organization of this paper is as follows. In Section2, we will consider the existence and stability of equilibria of system (7). The existence of Hopf bifurcation is also discussed. In Section 3, by use of normal form theory and central manifold argument, we illustrate the direction and stability of Hopf bifurcation. In Section 4, we provide an example with some numerical simulations to verify the theoretical results, and we also give some brief discussion.


## 2. Local Stability of Equilibria and the Existence of Hopf Bifurcations

In this section, we will finish two tasks: (a) investigating the existence and stability of equilibriums of system (7) and (b) studying the effect of time delay on the system (7); that is, we will choose τ as bifurcating parameter to analyze Hopf bifurcation.Let the right equations of system (7) equal zero; we get the following algebraic equations:
(8)N=N(α1-a11N-a12P),P=P(-α2+a21u-a22P),u=αN-αu.By simple computation, we know that the trivial equilibrium and boundary equilibrium of system (7) always exist with values E0=(0,0,0) and E1=(α1/a11,0,α1/a11), respectively. In addition, we have the following results.(i)
The eigenvalues of characteristic equations at the trivial equilibriumE0 are λ1=α1>0, λ2=-α2<0, and λ3=-α<0, which means that this equilibrium is always unstable.(ii)
The eigenvalues of characteristic equations at the boundary equilibriumE1 are λ1=0 and λ2=-α<0, and other eigenvalues are determined by λ+α1e-λτ=0. When τ=0, it is easy to see that λ3=-α1<0, which means that this equilibrium is locally stable; whereas when τ>0, the sign of the real part of the eigenvalues can not be determined, which means that this equilibrium may be locally stable or unstable.In fact, there exists a unique positive equilibriumE2=(N*,u*,P*) for system (7) provided that α1a21-α2a11>0, holds. Here
(9)N*=u*=α1a22+α2a12a11a22+a12a21,P*=α1a21-α2a11a11a22+a12a21.
Next, we always assume that(H1)
α1a21-α2a11>0 holds.The characteristic equation for system (7) at the equilibrium E2 takes the form
(10)λ3+d1λ2+(d2λ2+d3λ+d4)e-λτ+(d5λ+d6)e-2λτ=0,
where
(11)d1=α,d2=a11N*+a22P*,d3=αa11N*+αa22P*,d4=αa21a12N*P*,d5=a11a22N*P*,d6=αa11a22N*P*.Multiplyingeλτ on both sides of (10), we obtain equivalent characteristic equation as
(12)(λ3+d1λ2)eλτ+(d2λ2+d3λ+d4)+(d5λ+d6)e-λτ=0.Whenτ=0, characteristic equation (10) or (12) becomes
(13)λ3+(d1+d2)λ2+(d3+d5)λ+d4+d6=0.
It is easy to confirm that d1+d2>0, d4+d6>0 and (d1+d2)(d3+d5)>d4+d6. By the Routh-Hurwitz criterion we know that all the roots of (13) have negative real parts. Thus, the positive equilibrium E2 is locally asymptotically stable for τ=0.Next, we will consider the eigenvalues of (12) for τ>0. Suppose that there is a pure imaginary root λ=iω, ω>0, then we get
(14)(-iω3-d1ω2)(cosωτ+isinωτ)-d2ω2+d3iω+d4+(d5iω+d6)(cosωτ-isinωτ)=0.Separating the real and imaginary parts, we have(15)(-d1ω2+d6)cosωτ+(ω3+d5ω)sinωτ=d2ω2-d4,(d5ω-ω3)cosωτ-(d1ω2+d6)sinωτ=-d3ω.
By simple calculation, we can obtain the following equations:
(16)sinωτ=d2ω5+(d1d3-d4-d2d5)ω3+(d4d5-d3d6)ωω6+d12ω4-d52ω2-d62,cosωτ=(d3-d1d2)ω4+(d1d4-d2d6+d3d5)ω2+d4d6ω6+d12ω4-d52ω2-d62.Let(17)e1=d12,e2=-d52,e3=-d62,e4=d2,e5=d1d3-d4-d2d5,e6=d4d5-d3d6,e7=d3-d1d2,e8=d1d4-d2d6+d3d5,e9=d4d6,
then sinωτ, cosωτ can be written as
(18)sinωτ=e4ω5+e5ω3+e6ωω6+e1ω4+e2ω2+e3,(19)cosωτ=e7ω4+e8ω2+e9ω6+e1ω4+e2ω2+e3.By adding the square of (18) and (19), we obtain
(20)ω12+f1ω10+f2ω8+f3ω6+f4ω4+f5ω2+f6=0,
where
(21)f1=2e1-e42,f2=e12+2e2-2e4e5-e72,f3=2e3+2e1e2-2e4e6-2e7e8-e52,f4=2e1e3+e22-2e5e6-2e7e9-e82,f5=2e2e3-e62-2e8e9,f6=e32-e92.Denotez=ω2, then (20) becomes
(22)z6+f1z5+f2z4+f3z3+f4z2+f5z+f6=0.Let(23)G(z)=z6+f1z5+f2z4+f3z3+f4z2+f5z+f6,
then the following assumption holds true.(H2)
Equation (22) has at least one positive real root.In fact, if all the parameters of system (7) are given, it is easy to calculate the root of (22) by using a computer. Since limz→∞G(z)=+∞, we conclude that if  f6<0, then (22) has at least one positive real root. Without loss of generality, we assume that (22) has six positive roots, defined by z1, z2, z3, z4, z5, and z6, respectively. Then (20) has six positive roots:
(24)ω1=z1,ω2=z2,ω3=z3,ω4=z4,ω5=z5,ω6=z6.
By (19), we get
(25)cosωkτ=e7ωk4+e8ωk2+e9ωk6+e1ωk4+e2ωk2+e3.
If we denote
(26)τk(j)=1ωk{arccos(e7ωk4+e8ωk2+e9ωk6+e1ωk4+e2ωk2+e3)+2jπ},
where k=1,2,…,6 and j=0,1,2,…, then ±iωk is a pair of purely imaginary roots of (12). Define
(27)τ0=τk0(0)=mink∈{1,2,…,6}{τk(0)},ω0=ωk0.
Here, we use the method by [25–27], which is different from [31].In order to obtain the main result, it is necessary to make the following assumption.(H3)
Re(dλ/dτ)|τ=τ0≠0.Taking the derivative of λ with respect to τ in (12), it is easy to obtain
(28)(3λ2+2d1λ)eλτdλdτ+(λ3+d1λ2)eλτ(λ+τdλdτ)+(2d2λ+d3)dλdτ+d5e-λτdλdτ+(d5λ+d6)e-λτ(-λ-τdλdτ)=0,
which is equivalent to
(29)dλdτ=(λ(d5λ+d6)e-λτ-λ(λ3+d1λ2)eλτ)×((3λ2+2d1λ)eλτ+(λ3+d1λ2)τeλτ+d5e-λτ-(d5λ+d6)τe-λτ+2d2λ+d3)-1.
By (12), we have
(30)(dλdτ)-1=(3λ2+2d1λ)eλτ+d5e-λτ+2d2λ+d3λ(d5λ+d6)e-λτ-λ(λ3+d1λ2)eλτ-τλ(31)=(3λ2+2d1λ)eλτ+d5e-λτ+2d2λ+d3d2λ3+d3λ2+d4λ+2λ(d5λ+d6)e-λτ-τλ.
Take λ=iω into the above equation, we get
(32)((-3ω2+2d1iω)(cosωτ+isinωτ)+d5(cosωτ-isinωτ)+2d2iω+d3((-3ω2+2d1iω)(cosωτ+isinωτ))×(-d2iω3-d3ω2+d4iω+2iω(d5iω+d6)×(cosωτ-isinωτ)(-d2iω3-d3ω2+d4iω+2iω(d5iω+d6))-1-τiω=((-3ω2+d5)cosωτ-2d1ωsinωτ+d3+i[(-3ω2-d5)sinωτ+2d1ωcosωτ+2d2ω])×(+2d5ω2sinωτ)-1i(-d2ω3+d4ω+2d6ωcosωτ(+2d5ω2sinωτ)-1(-3d3ω2-2d5ω2cosωτ+2d6ωsinωτ)+i(-d2ω3+d4ω+2d6ωcosωτ+2d5ω2sinωτ)-1i(-d2ω3+d4ω+2d6ωcosωτ)-τiω.Let(33)Q=(-3d3ω2-2d5ω2cosωτ+2d6ωsinωτ)2+(-d2ω3+d4ω+2d6ωcosωτ+2d5ω2sinωτ)2>0.
Then we get
(34)QRe(dλdτ)-1=[(-3d3ω2+d5)cosωτ-2d1ωsinωτ+d3]×[-3d3ω2-2d5ω2cosωτ+2d6ωsinωτ]+[(-3ω2-d5)sinωτ+2d1ωcosωτ+2d2ω]×[-d2ω3+d4ω+2d6ωcosωτ+2d5ω2sinωτ].Note that(35)Sign{Re(dλdτ)|τ=τ0}=Sign{Re(dλdτ)-1|τ=τ0}.Now, we can use the following lemma to get our result.Lemma 1 (see [32]).
Consider the exponential polynomial(36)P(λ,e-λτ1,e-λτ2,…,e-λτm)=λn+p1(0)λn-1+⋯+pn-1(0)λ+pn(0)+[p1(1)λn-1+⋯+pn-1(1)λ+pn(1)]e-λτ1+⋯+[p1(m)λn-1+⋯+pn-1(m)λ+pn(m)]e-λτm,
where τi≥0(i=1,2,…,m) and  pj(i)(i=0,1,…,m;j=1,2,…,n) are constants. As  (τ1,τ2,…,τm) vary, the sum of the order of the zeros of  P(λ,e-λτ1,…,e-λτm) on the open right half plane can change only if a zero appears on or crosses the imaginary axis.Theorem 2.
Suppose that (H1), (H2), and (H3) hold; then the following results hold.(i)
The positive equilibriumE2 of system (7) (or the positive equilibrium  (N*,P*) of system (4)) is asymptotically stable for τ∈[0,τ0).(ii)
The positive equilibriumE2 of system (7) or system (4) undergoes a Hopf bifurcation when τ=τ0. That is, system (7) has a branch of periodic of solutions bifurcating from the positive equilibrium E2 near τ=τ0.


## 3. Direction and Stability of the Hopf Bifurcation

In this section, following the ideas of [33], we derive the explicit formulae for determining the properties of the Hopf bifurcation at critical value of τ0 by using the normal form and the center manifold theory. Throughout this section, we always assume that system (7) undergoes Hopf bifurcation at the positive equilibrium E2 for τ=τ0, and then ±iω0 is the corresponding purely imaginary roots of the characteristic equation at the positive equilibrium E2.Letx1=N-N*,x2=u-u*,x3=P-P*,x¯i(t)=xi(τt), τ=τ0+μ, dropping the bars for simplification of notations; then system (7) is transformed into functional differential equations in C=C([-1,0],R3) as
(37)x′(t)=Lμ(xt)+f(μ,xt),
where x(t)=(x1(t),x2(t),x3(t))T∈R3, Lμ:C→R, f:R×C→R3, and
(38)Lμ(ϕ)=(τ0+μ)[000α-α00a21P*0][ϕ1(0)ϕ2(0)ϕ3(0)]lll+(τ0+μ)[-a11N*0-a12N*00000-a22P*][ϕ1(-1)ϕ2(-1)ϕ3(-1)],(39)f(μ,ϕ)=(τ0+μ)[-a11ϕ1(0)ϕ1(-1)-a12ϕ1(0)ϕ3(-1)0a21ϕ2(0)ϕ3(0)-a22ϕ3(0)ϕ3(-1)],
where ϕ(θ)=(ϕ1(θ),ϕ2(θ),ϕ3(θ))T∈C. By the Riesz representation theorem, there exists a function  η(θ,μ) of bounded variation for  θ∈[-1,0], such that
(40)Lμϕ=∫-10‍dη(θ,μ)ϕ(θ),forϕ∈C.
In fact, we can choose
(41)η(θ,μ)=(τ0+μ)[000α-α00a21P*0]δ(θ)lll-(τ0+μ)[-a11N*0-a12N*00000-a22P*]δ(θ+1),
where δ is the Dirac delta function. For ϕ∈C1([-1,0],R3), define
(42)A(μ)ϕ={dϕ(θ)dθ,θ∈[-1,0),∫-10‍dη(s,μ)ϕ(s),θ=0,R(μ)ϕ={0,θ∈[-1,0),f(μ,ϕ),θ=0.
Then system (37) is equivalent to
(43)x′(t)=A(μ)xt+R(μ)xt,
where xt=x(t+θ),θ∈[-1,0]. For ψ∈C1([0,1],(R3)*), define
(44)A*ψ(s)={-dψ(s)ds,s∈(0,1],∫-10‍dηT(t,0)ψ(-t),s=0,
and a bilinear inner product
(45)〈ψ(s),ϕ(θ)〉=ψ¯(0)ϕ(0)-∫-10‍∫ξ=0θ‍ψ¯(ξ-θ)dη(θ)ϕ(ξ)dξ,
where η(θ)=η(θ,0). Then  A(0) and  A* are adjoint operators. By the discussion in Section 2, we know  ±iω0τ0 are eigenvalues of  A(0). Thus, they are also eigenvalues of  A*. We need to compute the eigenvector of  A(0) and  A* corresponding to  iω0τ0 and  -iω0τ0, respectively.Suppose thatq(θ)=(1,q1,q2)Teiθω0τ0 is the eigenvalues of A(0) corresponding to iω0τ0: then A(0)q(θ)=iω0τ0q(θ). It follows from the definition of A(0) and η(θ,μ) that
(46)τ0[000α-α00a21P*0]q(0)+τ0[-a11N*0-a12N*00000-a22P*]q(-1)=iω0τ0q(0),
because  q(-1)=q(0)e-iω0τ0; then we get
(47)[iω0+a11N*e-iω0τ00a12N*e-iω0τ0-αiω0+α00-a21P*iω0+a22P*e-iω0τ0][1q1q2]=[000],
where
(48)q1=αiω0+α(a12+a11)N*e-iω0τ0iω0a12N*exp-iω0τ0,q2=-iω0+a11N*e-iω0τ0a12N*exp-iω0τ0.Similarly, letq*(θ)=D(1,q1*,q2*)eiθω0τ0 be the eigenvalues of A* corresponding to -iω0τ0; according to the definition of A* we get
(49)τ0[0α00-αa21P*000]q*(0)+τ0[-a11N*00000-a12N*0-a22P*]q*(-1)=-iω0τ0q*(0),
where  q*(-1)=q*(0)e-iω0τ0, from which we obtain
(50)[-iω0+a11N*e-iω0τ0-α00-iω0+α-a21P*a12N*eiω0τ00-iω0+a22P*e-iω0τ0]×[1q1*q2*]=[000],q1*=-iω0+a11N*eiω0τ0α,q2*=--iω0+a22P*eiω0τ0a12N*eiω0τ0.By (45) we get
(51)〈q*(s),q(θ)〉=D¯(1,q1*¯,q2*¯)(1,q1,q2)T-∫-10‍∫ξ=0θ‍D¯(1,q1*¯,q2*¯)e-iω0τ0(ξ-θ)×dη(θ)(1,q1,q2)Teiω0τ0ξdξ=D¯{(1,q1,q2)T-∫-10‍(1,q1*¯,q2*¯)1+q1q¯1*+q2q¯2*-∫-10‍(1,q1*¯,q2*¯)θeiω0τ0θdη(θ)×(1,q1,q2)T-∫-10‍(1,q1*¯,q2*¯)}=D¯{-a22P*q2q2*¯)(-a11N*-a12N*q2-a22P*q2q2*¯)(P*q2q2*¯-a11N*-a12N*q2-a22P*q2q2*¯)e-iω0τ0(P*q2q2*¯-a11N*-a12N*q2}1+q1q¯1*+q2q¯2*+τ(P*q2q2*¯-a11N*-a12N*q2-a22P*q2q2*¯)e-iω0τ0(P*q2q2*¯-a11N*-a12N*q2},
then we choose
(52)D¯=1×({(-a11N*-a12N*q2-a22P*q2q¯2*)e-iω0τ0}(-a11N*-a12N*q2-a22P*q2q¯2*)1+q1q¯1*+q2q¯2*lll+τ(-a11N*-a12N*q2-a22P*q2q¯2*)e-iω0τ0)-1
such that 〈q*(s),q(θ)〉=1, 〈q*(s),q¯(θ)〉=0.Next, we will use the ideas in [33] to compute the coordinates describing center manifold C0 at μ=0. Define
(53)z(t)=〈q*,xt〉,W(t,θ)=xt(θ)-2Re{z(t)q(θ)}.On the center manifoldC0, we have
(54)W(z(t),z¯(t),θ)=W20(θ)z22+W11(θ)zz¯+W02z¯22+⋯,
where z and z¯ are local coordinates for center manifold C0 in the direction of q*(s) and q¯(θ). Note that W is real if xt is real. We consider only real solutions. For the solution xt∈C0 of (43), since μ=0, we have
(55)z′(t)=iω0τ0z+q¯*(0)f(0,W(z,z¯,0))+2Re{zq(θ)}=iω0τ0z+q¯*(0)f0(z,z¯)=iω0τ0z(t)+g(z,z¯,),
where
(56)g(z,z¯,)=q¯*(0)f0(z,z¯)=g20z22+g11zz¯+g02z¯22+g21z2z¯2+⋯.
From (53) and (54), we have
(57)xt=(x1t(θ),x2t(θ),x3t(θ))=W(t,θ)+2Re{z(t)q(θ)}=W(t,θ)+zq(θ)+zq(θ)¯=W20(θ)z22+W11(θ)zz¯+W02z¯22+(1,q1,q2)Teiθω0τ0z+(1,q¯1,q¯2)Teiθω0τ0z¯+⋯,
then we can obtain
(58)x1t(0)=W20(1)(0)z22+W11(1)(0)zz¯+W02(1)(0)z¯22+z+z¯+O(|(z,z¯)|3),x2t(0)=W20(2)(0)z22+W11(2)(0)zz¯+W02(2)(0)z¯22+q1z+q1¯z¯+O(|(z,z¯)|3),x3t(0)=W20(3)(0)z22+W11(3)(0)zz¯+W02(3)(0)z¯22+q2z+q2¯z¯+O(|(z,z¯)|3),x1t(-1)=W20(1)(-1)z22+W11(1)(-1)zz¯+W02(1)(-1)z¯22+ze-iω0τ0+z¯eiω0τ0+O(|(z,z¯)|3),x3t(-1)=W20(3)(-1)z22+W11(3)(-1)zz¯+W02(3)(-1)z¯22+zq2e-iω0τ0+zq2¯eiω0τ0+O(|(z,z¯)|3).
From the definition of  f(μ,xt), we have
(59)g(z,z¯)=q*¯(0)f0(z,z¯)=τ0D¯(1,q1*¯,q2*¯)×[-a11x1t(0)x1t(-1)-a12x1t(0)x3t(-1)0a21x3t(0)x2t(0)-a22x3t(0)x3t(-1)]=τ0D¯{z2[-a11e-iω0τ0-a12q2e-iω0τ0lll+a21q1q2q2*¯-a22q22q2*¯e-iω0τ0]+2zz¯-[Re{eiω0τ0}a11Re{eiω0τ0}-a12Re{q2e-iω0τ0}+a21q2*¯Re{q1q2¯}+a22q22q2*Re{eiω0τ0}]+z2¯[a22q22¯q2*¯eiω0τ0-a11eiω0τ0-a12q2¯eiω0τ0lll+a21q1q2¯q2*¯-a22q22¯q2*¯eiω0τ0]+12z2z¯[-a11(W20(1)(0)eiω0τ0+2W11(1)(0)e-iω0τ0+2W11(1)(-1)+W20(1)(-1))-a12(W20(1)(0)q2¯eiω0τ0+2W11(1)(0)q2e-iω0τ0+2W11(3)(-1)+W20(3)(-1))+a21q2*¯(W20(2)(0)q2¯+2W11(2)(0)q2+2W11(3)(0)q1+W20(3)(0)q1¯)-a22q2*¯(W20(3)(0)q2¯eiω0τ0+2W11(3)(0)q2e-iω0τ0+2W11(3)(-1)q2+W20(3)(-1)q2¯)]+⋯}.
Comparing the coefficients with (56), we obtain
(60)g20=2τ0D¯[-a11e-iω0τ0-a12q2e-iω0τ0+a21q1q2q2*¯-a22q22q2*¯e-iω0τ0],g11=2τ0D¯[-a11Re{eiω0τ0}-a12Re{q2e-iω0τ0}+a21q2*¯Re{q1q2¯}-a22q2q2¯q2*¯Re{eiω0τ0}],g02=2τ0D¯[a22q22¯q2*¯eiω0τ0-a11eiω0τ0-a12q2¯eiω0τ0+a21q1¯q2¯q2*¯-a22q22¯q2*¯eiω0τ0],g21=τ0D¯[-a11(W20(1)(0)eiω0τ0+2W11(1)(0)e-iω0τ0+2W11(1)(-1)+W20(3)(-1))-a12(W20(1)(0)q2¯eiω0τ0+2W11(1)(0)q2e-iω0τ0+2W11(3)(-1)+W20(3)(-1))+a21q2*¯(W20(2)(0)q2¯+2W11(2)(0)q2+2W11(3)(0)q1+W20(3)(0))-a22q2¯(W20(3)(0)q2*¯eiω0τ0+2W11(3)(0)q2e-iω0τ0+2W11(3)(-1)q2+W20(3)(-1)q2)].In order to determineg21 we need to compute W20(θ) and W11(θ). From (43) and (53), we have
(61)W˙=x˙t-z˙q-zq-˙llll={A(0)W-2Re{q¯*(0)f0q(θ)},θ∈[-1,0),A(0)W-2Re{q¯*(0)f0q(θ)}+f0,θ=0,llll=˙A(0)W-H(z,z¯,θ),
where
(62)H(z,z¯,θ)=H20(θ)z22+H11(θ)zz¯+H02z2¯2+⋯.
Note that on the center manifold C0 near to the origin,
(63)W˙=Wzz˙+Wz¯z¯˙,
thus we obtain
(64)(A(0)-2iω0τ0)W20(θ)=-H20(θ),A(0)W11(θ)=-H11(θ).
By (62) we know that for θ∈[-1,0),
(65)H(z,z¯,θ)=-q*¯(0)f0q(θ)-q0*f0¯q¯(θ)=-g(z,z¯)q(θ)-g(z,z¯)¯q(θ)¯.
Comparing the coefficients with (62), we get that
(66)H20(θ)=-g20q(θ)-g¯02q¯(θ),H11(θ)=-g11q(θ)-g¯11q¯(θ).
From (64), (66), and the definition of A, we have
(67)W˙20(θ)=2iω0τ0W20(θ)+g20q(θ)+g¯02q¯(θ).
Noting  q(θ)=q(0)eiω0τ0θ, hence
(68)W20(θ)=ig20ω0τ0q(0)eiω0τ0θ+ig¯023ω0τ0q¯(0)e-iω0τ0θ+E1e2iω0τ0θ,
where E1=(E1(1),E1(2),E1(3))T∈R3 is a constant vector. Similarly, we have
(69)W11(θ)=-ig11ω0τ0q(0)eiω0τ0θ+ig¯11ω0τ0q¯(0)e-iω0τ0θ+E2,
where E2=(E2(1),E2(2),E2(3))T∈R3 is a constant vector. In the following, we will find out E1 and E2. From the definition of A and (64), we can obtain
(70)∫-10‍dη(θ)W20(θ)=2iω0τ0W20(θ)-H20(0),(71)∫-10‍dη(θ)W11(θ)=-H11(θ),
where  η(θ)=η(θ,0). From (61) and (62) we have
(72)H(z,z¯,0)=-2Re{q¯*(0)f0q(0)}+f0=-q¯*(0)f0q(0)-q*(0)f0¯q¯(0)+f0=-g(z,z¯)q(0)-g¯(z,z¯)q¯(0)+f0.
That is,
(73)H20(θ)z22+H11(θ)zz¯+H02z2¯2+⋯=-q(0)(g20z22+g11zz¯+g02z2¯2+⋯)-q¯(0)(g¯20z¯22+g¯11zz¯+g02¯z22+⋯)+f0.
By (39) and (53) we have
(74)f0=τ0[-a11x1t(0)x1t(-1)-a12x1t(0)x3t(-1)0a21x2t(0)x3t(0)-a22x3t(0)x3t(-1)],xt(θ)=W(t,θ)+2Re{z(t)q(θ)}=W(t,θ)+z(t)q(θ)+z¯(t)q¯(θ)=W20(θ)z22+W11(θ)zz¯+z(t)q(θ)+z¯(t)q¯(θ)+⋯.Thus,(75)H20(0)=-g20q(0)-g¯02q¯(0)+2τ0×[-a11e-iω0τ0-a12q2e-iω0τ00a21q1q2-a22q22e-iω0τ0],(76)H11(0)=-g11q(0)-g¯11q¯(0)+2τ0×[-a11Re{eiω0τ0}-a12Re{q2e-iω0τ0}0a21Re{q1q¯2}-a22q22Re{eiω0τ0}].Sinceiω0τ0 is the eigenvalues of A(0) and q(0) is the corresponding eigenvector, we obtain
(77)(iω0τ0I-∫-10‍eiω0τ0θdη(θ))q(0)=0,(-iω0τ0I-∫-10‍e-iω0τ0θdη(θ))q¯(0)=0.
Thus, substituting (68) and (75) into (70), we have
(78)(2iω0τ0I-∫-10‍e2iω0τ0θdη(θ))E1=2τ0[-a11e-iω0τ0-a12q2e-iω0τ00a21q1q2-a22q22e-iω0τ0],
or
(79)[2iω0+a11N*e-2iω0τ00a12N*e-2iω0τ0-α2iω0+α00-a21P*2iω0+a22P*e-2iω0τ0]E1=2[-a11e-iω0τ0-a12q2e-iω0τ00a21q1q2-a22q22e-iω0τ0].From which we can get(80)E1=2[2iω0+a11N*e-2iω0τ00a12N*e-2iω0τ0-α2iω0+α00-a21P*2iω0+a22P*e-2iω0τ0]-1×[-a11e-iω0τ0-a12q2e-iω0τ00a21q1q2-a22q22e-iω0τ0].Similarly, substituting (69) and (76) into (71), we can get
(81)[a11N*0a12N*-αα00-a21P*a22P*]E2=2[-a11Re{eiω0τ0}-a12Re{q2e-iω0τ0}0a21Re{q1q¯2}-a22q22Re{eiω0τ0}],
or
(82)E2=2[a11N*0a12N*-αα00-a21P*a22P*]-1×[-a11Re{eiω0τ0}-a12Re{q2e-iω0τ0}0a21Re{q1q¯2}-a22q22Re{eiω0τ0}].Thus, we can determineW20(θ) and W11(θ) from (68) and (69). Furthermore, g21 can be expressed by the parameters and delay. Thus, we can compute the following values:
(83)c1(0)=i2ω0τ0(g20g11-2|g11|2-|g02|23)+g212,μ2=-Re{c1(0)}Re{λ′(τ0)},β2=2Re{c1(0)},T2=-Im{c1(0)}+μ2Im{λ′(τ0)}ω0τ0,
which determine the qualities of bifurcating periodic solution in the center manifold at critical value τ0; that is, μ2 determine the direction of the Hopf bifurcation: if μ2>0(μ2<0), then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solution exists for τ>τ0(τ<τ0); β2 determines the stability of the bifurcating periodic solution: the bifurcating periodic solution is stable (unstable) if β2<0(β2>0); and T2 determines the period of the bifurcating periodic solution: the period increases (decreases) if T2>0(<0).


## 4. Numerical Investigations and Discussion

In this paper, we propose a two-dimensional predatory-prey model with discrete and distributed delay. Then, by introducing a new variable, the original system is transformed into an equivalent three-dimensional system. In Section2, we analyze the existence and local stability of the equilibria of the three-dimensional system. The condition for the existence of a Hopf bifurcation is also obtained. In Section 3, by the use of normal form theory and central manifold argument, we establish the formulae for the direction and the stability of the Hopf bifurcation.In order to confirm our main results obtained in this work, we consider the following special system:(84)N′(t)=N(t)[1-0.5N(t-τ)-2P(t-τ)],P′(t)=P(t)[-1+0.8u(t)-0.5P(t-τ)],u′(t)=2N(t)-2u(t).
By simple calculation, it is easy to see that model (84) exists a unique positive equilibrium E2 and E2=(50/37,6/37,50/37). Note that the parameter set provided in model (84) satisfies the conditions of Theorem 2. When τ=1.1, the positive equilibrium E2 is asymptotically stable, as shown in Figures 1(a) and 1(c). It follows from the discussion in Section 2 that ω0≈0.7246, τ0≈1.1746, and λ′(τ0)=0.09345-0.07497i. Thus, E2 is stable when 0≤τ<τ0, as indicated in Figures 1(a) and 1(c).The stability of unique positive equilibriumE2. (a), (c) The equilibrium E2 is stable for τ=1.1. (b), (d) The equilibrium E2 is unstable and a stable periodic solution appears for τ=1.4.
(a)
(b)
(c)
(d)Whenτ passes through the critical value τ0, E2 loses its stability and a Hopf bifurcation occurs, that is, a family of periodic solutions bifurcate from E2, as shown in Figures 1(b) and 1(d). Since μ2>0 and β2<0, the Hopf bifurcation is supercritical and the direction of the bifurcation is τ>τ0 and these bifurcating periodic solutions from E2 are stable; please see Figures 1(b) and 1(d) and Figure 2(a) for τ=1.5. Note that the model (84) may have very complex dynamics if we choose the time delay τ as a bifurcation parameter. It follows from Figure 2 that the period of periodic solution is doubled as τ increase, and around τ=1.6 (Figure 2(b)). If the time delay τ is increasing further, a periodic solution with 4-time period appears around τ=1.75 (Figure 2(c)). Finally, the chaotic solution exists once the time delay reaches around τ=1.8 (Figure 2(d)).Bifurcation of model (84) with bifurcation parameter τ.
(a)
(b)
(c)
(d)Both our theoretical and numerical results show that the positive equilibrium is asymptotically stable ifτ<τ0, which indicates that the dynamical behavior is simple for the considered system. However, if τ>τ0, bifurcation and chaos may occur, which means that the considered system can take on very complex dynamics, and this may explain some complex phenomenon in the natural world.



---

*Source: 454097-2013-06-25.xml*