# Stability and Hopf Bifurcation in a Modified Holling-Tanner Predator-Prey System with Multiple Delays

Authors: Zizhen Zhang; Huizhong Yang; Juan Liu
Journal: Abstract and Applied Analysis
 (2012)

---

## Abstract

A modified Holling-Tanner predator-prey system with multiple delays is investigated. By analyzing the associated characteristic equation, the local stability and the existence of periodic solutions via Hopf bifurcation with respect to both delays are established. Direction and stability of the periodic solutions are obtained by using normal form and center manifold theory. Finally, numerical simulations are carried out to substantiate the analytical results.


---

## Body

## 1. Introduction

Predator-prey dynamics has long been and will continue to be of interest to both applied mathematicians and ecologists due to its universal existence and importance [1, 2]. Many population models investigating the dynamic relationship between predators and their preys have been proposed and studied. For example, Lotka-Volterra model [3–5], Leslie-Gower model [6–10], and Holling-Tanner model [11–16]. Among these widely used models, Holling-Tanner model plays a special role in view of the interesting dynamics it possesses. Holling-Tanner model for predator-prey interaction is governed by the following nonlinear coupled ordinary differential equations:
(1.1)dXdT=rX(1-XK)-mXYa+X,dYdT=Y[s(1-hYX)],
where X and Y denote the population densities of prey species and predator species at time T, respectively. The first equation in system (1.1) shows that the prey grows logistically with the carrying capacity K and the intrinsic growth rate r in the absence of the predator. And the growth of the prey is hampered by the predator at a rate proportional to the functional response mX/(a+X) in the presence of the predator. The second equation shows that the predator consumes the prey according to the functional response mX/(a+X) and grows logistically with the intrinsic growth rate s and carrying capacity X/h proportional to the number of the prey. The parameter m denotes the maximal predator per capita consumption rate. a is a saturation value; it corresponds to the number of prey necessary to achieve one half the maximum rate m. The parameter h denotes the number of prey required to support one predator at equilibrium when y equals X/h.Recently, there has been considerable interest in predator-prey systems with the Beddington-DeAngelis functional response. And it has been shown that the predator-prey systems with the Beddington-DeAngelis functional response have rich but biologically reasonable dynamics. For more details about this functional response one can refer to [17–21]. Zhang [16], Lu and Liu [22] considered the following modified Holling-Tanner delayed predator-prey system:
(1.2)dXdT=rX(1-XK)-αXYa+bX+cY,dYdT=Y[s(1-hY(t-τ-)X(t-τ-))],
where τ- is incorporated in the negative feedback of the predator density. αXY/(a+bX+cY) is the Beddington-DeAngelis functional response. The parameters α, a, b, and c are assumed to be positive. α is the maximum value at which per capita reduction rate of the prey can attain. a measures the extent to which environment provides protection to the prey. b describes the effect of handling time on the feeding rate. c describes the magnitude of interference among predators. Zhang [16] investigated the local Hopf bifurcation of system (1.2). Lu and Liu [22] proved that system (1.2) is permanent under some conditions and obtained the sufficient conditions of local and global stability of system (1.2). Since both the species are growing logistically, it is reasonable to assume delay in prey species as well. Based on this consideration, we incorporate the negative feedback of the prey density into system (1.2) and obtain the following system:
(1.3)dXdT=rX(1-X(T-T1)K)-αXYa+bX+cY,dYdT=Y[s(1-hY(T-T2)X(T-T2))],
where T1 and T2 are the feedback time delays of the prey density and the predator density respectively. Let X=Kx, Y=(rK/α)y, t=rT, τ1=rT1, τ2=rT2, system (1.3) can be transformed into the following nondimensional form:
(1.4)dxdt=x(1-x(t-τ1))-xya1+bx+c1y,dydt=y[δ-βy(t-τ2)x(t-τ2)],
where a1=a/K,c1=cr/α,δ=s/r,β=sh/α are the non-dimensional parameters and they are positive.The main purpose of this paper is to consider the effect of multiple delays on system (1.4). The local stability of the positive equilibrium and the existence of Hopf bifurcation are investigated. By employing normal form and center manifold theory, the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions are determined. Finally, some numerical simulations are also included to illustrate the theoretical analysis.


## 2. Local Stability and the Existence of Hopf Bifurcation

In this section, we study the local stability of each of feasible equilibria and the existence of Hopf bifurcation at the positive equilibrium. Obviously, system (1.4) has a unique boundary equilibrium E1(1,0) and a unique positive equilibrium E*(x*,y*), where
(2.1)x*=-[(a1-b)β+(1-c1)δ]+[(a1-b)β+(1-c1)δ]2+4a1β(bβ+c1δ)2(bβ+c1δ),y*=δβx*.
The Jacobian matrix of system (1.4) at E1 takes the form
(2.2)J(E1)=(-e-λτ1-1a1+b0δ).
The characteristic equation of system (1.4) at E1 is of the form
(2.3)(λ-δ)(λ+e-λτ1)=0.
Clearly, the boundary equilibrium E1(1,0) is unable.Next, we discuss the existence of Hopf bifurcation at the positive equilibriumE(x*,y*). Let x(t)=z1(t)+x*, y(t)=z2(t)+y*, and still denote z1(t) and z2(t) by x(t) and y(t), respectively, then system (1.4) becomes
(2.4)dxdt=a11x(t)+a12y(t)+b11x(t-τ1)+∑i+j+k≥2‍f1ijkxiyjxk(t-τ1),dydt=c21x(t-τ2)+c22y(t-τ2)+∑i+j+k≥2‍f2ijkyixj(t-τ2)yk(t-τ2),
where
(2.5)a11=bx*y*(a1+bx*+c1y*)2,a12=-(a1+bx*)x*(a1+bx*+c1y*)2,b11=-x*,c21=δ2β,c22=-δ,f1=x(1-x(t-τ1))-xya1+bx+c1y,f2=y(δ-βy(t-τ2)x(t-τ2)),f1ijkxiyjxk(t-τ1)=1i!j!k!∂i+j+kf1∂xi∂yj∂xk(t-τ1)|(x*,y*),f2ijkyixj(t-τ2)yk(t-τ2)=1i!j!k!∂i+j+kf2∂yi∂xj(t-τ2)∂yk(t-τ2)∣(x*,y*).
Then we can obtain the linearized system of system (1.4)
(2.6)dxdt=a11x(t)+a12y(t)+b11x(t-τ1),dydt=c21x(t-τ2)+c22y(t-τ2).
The characteristic equation of system (2.6) is
(2.7)λ2-a11λ-b11λe-λτ1+(a11c22-a12c21-c22λ)e-λτ2+b11c22e-λ(τ1+τ2)=0.Case 1.
τ
1
=
τ
2
=
τ
=
0.
The characteristic equation of system (1.4) becomes
(2.8)λ2+(A+B+D)λ+C+E=0,
where
(2.9)A=-a11,B=-b11,C=a11c22-a12c21,D=-c22,E=b11c22.
It is easy to verify that
(2.10)C+E=δx*+a1δ2x*β(a1+bx*+c1y*)2>0.
Therefore, if (H1):A+B+D>0, the roots of (2.8) must have negative real parts. Then, we know that the positive equilibrium E*(x*,y*) of system (1.4) is locally stable in the absence of delay, if (H1) holds.Case 2.
τ
1
=
τ
2
=
τ
>
0.
The associated characteristic equation of the system is(2.11)λ2+A1λ+(B1+C1λ)e-λτ+D1e-2λτ=0,
where
(2.12)A1=-a11,B1=a11c22-a12c21,C1=-(b11+c22),D1=b11c22.
Multiplying eλτ on both sides of (2.11), we can obtain
(2.13)(λ2+A1λ)eλτ+(B1+C1λ)+D1e-λτ=0.
Now, for τ>0, if λ=iω(ω>0) be a root of (2.13). Then, we have
(2.14)(D1-ω2)cosτω-A1ωsinτω=-B1,(D1+ω2)sinτω-A1ωcosτω=C1ω.
It follows from (2.14) that
(2.15)sinτω=C1ω2+(A1B1-C1D1)ωω4+A12ω2-D12,cosτω=(B1-A1C1)ω2+B1D1ω4+A12ω2-D12.
Then we have
(2.16)ω8+e3ω6+e2ω4+e1ω2+e0=0,
where
(2.17)e3=2A12-C12,e2=A14+2C12D1-A12C12-B12-2D12,e1=4A1B1C1D1-A12B12-C12D12-2A12D12-2B12D1,e0=D14-B12D12.
Let v=ω2, then (2.16) becomes
(2.18)v4+e3v3+e2v2+e1v+e0=0.
Next, we give the following assumption. (H2): (2.18) has at least one positive real root.Suppose that(H2) holds. Without loss of generality, we assume that (2.18) has four real positive roots, which are defined by v1, v2, v3, and v4, respectively. Then (2.16) has four positive roots ωk=vk,k=1,2,3,4. Therefore,
(2.19)τk(j)=1ωkarccos((B1-A1C1)ωk2+B1D1ωk4+A12ωk2-D12+2jπ),k=1,2,3,4;j=0,1,2…
Then we can know that ± iωk are a pair of purely imaginary roots of (2.11) with τ=τk(j). Define
(2.20)τ0=τk(0)=min{τk(0)},ω0=ωk0,k=1,2,3,4.
Let λ(τ)=α(τ)+iω(τ) be the root of (2.11) near τ=τ0 which satisfies α(τ0)=0,ω(τ0)=ω0. Taking the derivative of λ with respect to τ in (2.13), we obtain
(2.21)C1dλdτ+(2λ+A1)eλτdλdτ+(λ2+A1λ)eλτ(λ+τdλdτ)-D1e-λτ(λ+τdλdτ)=0.
it follows that
(2.22)dλdτ=λ(D1e-λτ-(λ2+A1λ)eλτ)(2λ+A1)eλτ+C1-τ(D1e-λτ-(λ2+A1λ)eλτ).
Thus
(2.23)[dλdτ]-1=(2λ+A1)eλτ+C1D1λe-λτ-(λ3+A1λ2)eλτ-τλ.
Let
(2.24)Λ1=(D1ω0-ω03)sinτ0ω0+A1ω02cosτ0ω0,Λ2=(D1ω0+ω03)cosτ0ω0+A1ω02sinτ0ω0,Λ3=A1cosτ0ω0-2ω0sinτ0ω0+C1,Λ4=A1sinτ0ω0-2ω0cosτ0ω0.
Substitute λ=iω0(ω0>0) into (2.23), we can get
(2.25)Re[dλ(τ0)dτ]-1=Re[(2λ+A1)eλτ+C1D1λe-λτ-(λ3+A1λ2)eλτ]λ=iω0=Λ1×Λ3+Λ2×Λ4Λ12+Λ22.
Noting that
(2.26)signRe[dλ(τ0)dτ]=sign[dReλ(τ0)dτ]-1.
Therefore, we make the following assumption in order to give the main results: (H3):Λ1×Λ3+Λ2×Λ4≠0. Then, by Corollary 2.4 in [23, 24], we have the following theorem.Theorem 2.1.
For system (1.4), if the conditions (H1)–(H3) hold, then the equilibrium E*(x*,y*) of system (1.4) is asymptotically stable for τ∈[0,τ0) and unstable when τ>τ0. And system (1.4) has a branch of periodic solution bifurcation from the zero solution near τ=τ0.Case 3.
τ
1
≠
τ
2, τ1>0 and τ2>0.
The associated characteristic equation of the system is(2.27)λ2+A2λ+B2λe-λτ1+(C2+D2λ)e-λτ2+E2e-λ(τ1+τ2),
where
(2.28)A2=-a11,B2=-b11,C2=a11c22-a12c21,D2=-c22,E2=b11c22.
We consider (2.27) with τ2 in its stable interval, regarding τ1 as a parameter. Without loss of generality, we consider system (1.4) under the case considered in [16], and τ2∈[0,τ20). τ20 is defined as in [16] and can be obtained by
(2.29)τ20=1ω+arccos[(a11c22+b11c22-a12c21)ω+2-(a11+b11)c22ω+2(a11c22+b11c22-a12c21)2+(c22ω+)2],
with
(2.30)ω+=-((a11+b11)2-c222)+((a11+b11)2-c222)2+4(a11c22+b11c22-a12c21)22.
Let λ=iω(ω>0) be a root of (2.27). Then we obtain
(2.31)(B2ω-E2sinτ2ω)sinτ1ω+E2cosτ2ωcosτ1ω=ω2-C2cosτ2ω-D2ωsinτ2ω,(B2ω-E2sinτ2ω)cosτ1ω-E2cosτ2ωsinτ1ω=C2sinτ2ω-D2ωcosτ2ω-A2ω.
It follows from (2.31) that
(2.32)sinτ1ω=M1N1-M2N2M12+M22,cosτ1ω=M1N2+M2N1M12+M22,
With
(2.33)M1=B2ω-E2sinτ2ω,M2=E2cosτ2ω,N1=ω2-C2cosτ2ω-D2ωsinτ2ω,N2=-A2ω+C2sinτ2ω-D2ωcosτ2ω.
Then we have
(2.34)P1(ω)+P2(ω)sinτ2ω+P3(ω)cosτ2ω=0,
where
(2.35)P1(ω)=ω4+(A22+D22-B22)ω2+C22-E22,P2(ω)=-2D2ω3-2A2C2ω+2B2E2ω,P3(ω)=2(A2D2-C2)ω2.
Suppose that (H4): (2.34) has at least finite positive roots. If (H4) holds, we define the roots of (2.34) as ω1,ω2,…,ωk. Then, for every fixed ωi(i=1,2,…,k), there exists a sequence {τ1i(j)∣j=1,2,…} which satisfies (2.34). Let
(2.36)τ1*=min{τ1i(j)∣i=1,2,…,k,j=0,1,…}.
When τ1=τ1*, (2.27) has a pair of purely imaginary roots ±iω* for τ2∈[0,τ20).To verify the transversality condition of Hopf bifurcation, we take the derivative ofλ with respect to τ1 in (2.27), we can obtain
(2.37)dλdτ1=λe-λτ1(B2λ+E2e-λτ2)(2λ+A2+B2e-λτ1+D2e-λτ2)-τ1e-λτ1(B2λ+E2e-λτ2)-τ2e-λτ2(C2+D2λ+E2e-λτ1).
Thus
(2.38)[dλdτ1]-1=2λ+A2+B2e-λτ1+D2e-λτ2-τ2e-λτ2(C2+D2λ+E2e-λτ1)λe-λτ1(B2λ+E2e-λτ2)-τ1λ.
Substitute λ=iω*(ω*>0) into (2.38), we can get
(2.39)Re[dλ(τ1*)dτ1]-1=Δ1×Δ3+Δ2×Δ4Δ12+Δ22,
where
(2.40)Δ1=E2ω*cosτ2ω*sinτ1*ω*-ω*cosτ1*ω*(B2ω*-E2ω*sinτ2ω*),Δ2=E2ω*cosτ2ω*cosτ1*ω*+ω*sinτ1*ω*(B2ω*-E2ω*sinτ2ω*),Δ3=A2+(D2-τ2C2)cosτ2ω*-τ2D2ω*sinτ2ω*+τ2E2sinτ2ω*sinτ1*ω*+(B2-τ2E2cosτ2ω*)cosτ1*ω*,Δ4=2ω*+(τ2C2+τ2D2ω*-D2)sinτ2ω*+τ2E2sinτ2ω*cosτ1*ω*-(B2-τ2E2cosτ2ω*)sinτ1*ω*.
Next, we make the following assumption: (H5): Δ1×Δ3+Δ2×Δ4≠0.Thus, by the discussion above and by the general Hopf bifurcation theorem for FDEs in Hale [25], we have the following results.Theorem 2.2.
Forτ2∈(0,τ20), τ20 is defined by (2.29). If the conditions (H4)-(H5) hold, then the equilibrium E*(x*,y*) of system (1.4) is asymptotically stable for τ1∈[0,τ1*) and unstable when τ>τ1*. System (1.4) has a branch of periodic solution bifurcation from the zero solution near τ=τ1*.


## 3. Direction and Stability of Bifurcated Periodic Solutions

In this section, we shall investigate the direction of the Hopf bifurcation and the stability of bifurcating periodic solution of system (1.4) w.r. to τ1 for τ2∈(0,τ20), and τ20 is defined by (2.29). The idea employed here is the normal form and center manifold theory described in Hassard et al. [26]. Throughout this section, it is considered that system (1.4) undergoes the Hopf bifurcation at τ1=τ1*, τ2∈(0,τ20) at E*(x*,y*). Let τ1=τ1*+μ,μ∈R so that the Hopf bifurcation occurs at μ=0. Without loss of generality, we assume that τ2*<τ1*, where τ2*∈(0,τ20).Letu1(t)=x(t)-x*, u2(t)=y(t)-y*, and rescaling the time delay t→(t/τ1), Then system (1.4) can be transformed into an FDE in C=C([-1,0],R2) as:
(3.1)u˙(t)=Lμut+F(μ,ut),
where u(t)=(u1(t),u2(t))T∈R2 and Lμ:C→R2, F:R×C→R2 are given, respectively, by
(3.2)Lμϕ=(τ1*+μ)(A′ϕ(0)+C′ϕ(-τ2*τ1*)+B′ϕ(-1)),F(μ,ϕ)=(τ1*+μ)(f1,f2)T,
with
(3.3)A′=(a11a1200),B′=(b11000),C′=(00c21c22),f1=g1ϕ12(0)+g2ϕ1(0)ϕ2(0)+g3ϕ22(0)+g4ϕ1(0)ϕ1(-1)+h1ϕ13(0)+h2ϕ12(0)ϕ2(0)+h3ϕ1(0)ϕ22(0)+h4ϕ23(0)+⋯,f2=g1′ϕ12(-τ2*τ1*)+g2′ϕ1(-τ2*τ1*)ϕ2(0)+g3′ϕ1(-τ2*τ1*)ϕ2(-τ2*τ1*),+g4′ϕ2(0)ϕ2(-τ2*τ1*)+h1′ϕ13(-τ2*τ1*)+h2′ϕ12(-τ2*τ1*)ϕ2(0),+h3′ϕ12(-τ2*τ1*)ϕ2(-τ2*τ1*)+h4ϕ23(0)+⋯,
where
(3.4)g1=by*(a1+c1y*)(a1+bx*+c1y*)3,g2=-a12+a1bx*+a1c1y*+2bc1x*y*(a1+bx*+c1y*)3,g3=c1x*(a1+bx*)(a1+bx*+c1y*)3,g4=-1,h1=-b2y*(a1+c1y*)(a1+bx*+c1y*)4,h2=a12b+a1b2x*+2b2c1x*y*-bc12y*2(a1+bx*+c1y*)4,h3=a12c1+a1c12y*+2bc12x*y*-b2c1x*2(a1+bx*+c1y*)4,h4=-c12x*(a1+bx*)(a1+bx*+c1y*)4,g1′=-βy*2x*3,g2′=βy*x*2,g3′=βy*x*2,g4′=-βx*,h1′=βy*2x*4,h2′=-βy*x*3,h3′=-βy*x*3.
Using Riesz representation theorem, there exists a 2×2 matrix function η(θ,μ),θ∈[-1,0] whose elements are of bounded variation, such that
(3.5)Lμϕ=∫-10‍dη(θ,μ)ϕ(θ),ϕ∈C([-1,0],R2).
In fact, choosing
(3.6)η(θ,μ)={(τ1*+μ)(A′+C′+B′),θ=0,(τ1*+μ)(C′+B′),θ∈[-τ2*τ1,0),(τ1*+μ)B′,θ∈(-1,-τ2*τ1),0,θ=-1.
For ϕ∈C([-1,0]), we define
(3.7)A(μ)ϕ={dϕ(θ)dθ,-1≤θ<0,∫-10‍dη(θ,μ)ϕ(θ),θ=0,R(μ)ϕ={0,-1≤θ<0,F(μ,ϕ),θ=0.
Then system (3.1) can be transformed into the following operator equation
(3.8)u˙(t)=A(μ)ut+R(μ)ut,
where ut=u(t+θ)=(u1(t+θ),u2(t+θ)).Forφ∈C1([0,1],(R2)*), where (R2)* is the 2-dimensional space of row vectors, we further define the adjoint operator A* of A(0):
(3.9)A*(φ)={-dφ(s)ds,0<s≤1,∫-10‍φ(-ξ)dη(ξ,0),s=0,
and a bilinear inner product:
(3.10)〈φ(s),ϕ(θ)〉=φ-T(0)ϕ(0)-∫θ=-10‍∫ξ=0θ‍φ-T(ξ-θ)dη(θ)ϕ(ξ)dξ,
where η(θ)=η(θ,0).Since ±iω*τ1* are eigenvalues of A(0), they are also eigenvalues of A*. Let q(θ)=(1,q2)Teiω*τ1*θ be the eigenvectors of A(0) corresponding to iω*τ1* and q*(s)=(1/ρ)(1,q2*)Teiω*τ1*s be the eigenvectors of A* corresponding to -iω*τ1*. By a simple computation, we can get
(3.11)q2=iω*-a11-b11eiω*τ1*a12,q2*=-iω*+a11+b11eiω*τ1*c21eiω*τ2*,ρ-=1+q2q-2*+b11τ1*e-iω*τ1*+c21τ2*q-2*e-iω*τ2*+c22τ2*q2q-2*e-iω*τ2*.
Then 〈q*,q〉=1, 〈q*,q-〉=0.In the remainder of this section, Following the algorithms given in [26] and using similar computation process to that in [16], we can get that the coefficients which will be used to determine the important qualities of the bifurcating periodic solutions,
(3.12)g20=2τ1*ρ-{(g1′(q(1)(-τ2*τ1*))2+g2′q(1)(-τ2*τ1*)q(2)(0)+g3′q(1)(-τ2*τ1*)q(2)(-τ2*τ1*)+g4′q(2)(0)q(2)(-τ2*τ1*))g1+g2q(2)(0)+g3(q(2)(0))2+g4q(1)(-1)+q-2*(g1′(q(1)(-τ2*τ1*))2+g2′q(1)(-τ2*τ1*)q(2)(0)+g3′q(1)(-τ2*τ1*)q(2)(-τ2*τ1*)+g4′q(2)(0)q(2)(-τ2*τ1*))(g1′(q(1)(-τ2*τ1*))2+g2′q(1)(-τ2*τ1*)q(2)(0)+g3′q(1)(-τ2*τ1*)q(2)(-τ2*τ1*)},g11=τ1*ρ-{(q(1)(-τ2*τ1*)q-(2)(-τ2*τ1*)+q-(1)(-τ2*τ1*)q(2)(-τ2*τ1*))2g1+g2(q(2)(0)+q-(2)(0))+2g3q(2)(0)q-(2)(0)+g4(q(1)(-1)+q-(1)(-1))+q-2*(2g1′q(1)(-τ2*τ1*)q-(1)(-τ2*τ1*)+g2′(q(1)(-τ2*τ1*)q-(2)(0)+q-(1)(-τ2*τ1*)q(2)(0))+g3′(q(1)(-τ2*τ1*)q-(2)(-τ2*τ1*)+q-(1)(-τ2*τ1*)q(2)(-τ2*τ1*))+g4′(q(2)(0)q-(2)(-τ2*τ1*)+q-(2)(0)q(2)(-τ2*τ1*)))},g02=2τ1*ρ-{(g1′(q-(1)(-τ2*τ1*))2+g2′q-(1)(-τ2*τ1*)q-(2)(0)+g3′q-(1)(-τ2*τ1*)q-(2)(-τ2*τ1*)+g4′q-(2)(0)q-(2)(-τ2*τ1*))g1+g2q-(2)(0)+g3(q-(2)(0))2+g4q-(1)(-1)+q-2*(g1′(q-(1)(-τ2*τ1*))2+g2′q-(1)(-τ2*τ1*)q-(2)(0)+g3′q-(1)(-τ2*τ1*)q-(2)(-τ2*τ1*)+g4′q-(2)(0)q-(2)(-τ2*τ1*))},g21=2τ1*ρ-{((q(1)(-τ2*τ1*))2q-(2)(0)+2q(1)(-τ2*τ1*)q-(1)(-τ2*τ1*)q(2)(0))g1(W20(1)(0)+2W11(1)(0))+g2(12W20(2)(0)+W11(2)(0)+12W20(1)(0)q-(2)(0)+W11(1)(0)q(2)(0))+g3(W20(2)(0)q-(2)(0)+2W11(2)(0)q(2)(0))+g4(12W20(1)(-1)+W11(1)(-1)+12W20(1)(0)q-(1)(-1)+W11(1)(0)q(1)(-1))+3h1+h2(2q(2)(0)+q-(2)(0))+h3(2q(2)(0)q-(2)(0)+(q(2)(0))2)+3h4(q(2)(0))2q-(2)(0)+q-2*(g1′(W20(1)(-τ2*τ1*)q-(1)(-τ2*τ1*)+2W11(1)(-τ2*τ1*)q(1)(-τ2*τ1*))+g2′(12W20(1)(-τ2*τ1*)q-(2)(0)+W11(1)(-τ2*τ1*)q(2)(0)+12W20(2)(0)q-(1)(-τ2*τ1*)+W11(2)(0)q(1)(-τ2*τ1*))+g3′(12W20(1)(-τ2*τ1*)q-(2)(-τ2*τ1*)+W11(1)(-τ2*τ1*)q(2)(-τ2*τ1*)+12W20(2)(-τ2*τ1*)q-(1)(-τ2*τ1*))+g4′(12W20(2)(0)q-(2)(-τ2*τ1*)+W11(2)(0)q(2)(-τ2*τ1*)+12W20(2)(-τ2*τ1*)q-(2)(0)+W11(2)(-τ2*τ1*)q(2)(0))+3h1′(q(1)(-τ2*τ1*))2q-(1)(-τ2*τ1*)+h2′((q(1)(-τ2*τ1*))2q-(2)(0)+2q(1)(-τ2*τ1*)q-(1)(-τ2*τ1*)q(2)(0))+h3′((q(1)(-τ2*τ1*))2q-(2)(-τ2*τ1*)+2q(1)(-τ2*τ1*)q-(1)(-τ2*τ1*)q(2)(-τ2*τ1*)))},
with
(3.13)W20(θ)=ig20q(0)τ1*ω*eiτ1*ω*θ+ig-02q-(0)3τ1*ω*e-iτ1*ω*θ+E20e2iτ1*ω*θ,W11(θ)=-ig11q(0)τ1*ω*eiτ1*ω*θ+ig-11q-(0)τ1*ω*e-iτ1*ω*θ+E11,
where E20 and E11 can be computed as the following equations, respectively,
(3.14)E20=2(E20(1)E20(2))×(2iω*-a11-b11e-2iτ1*ω*-a12-c21e-2iτ2*ω*2iω*-c22e-2iτ2*ω*)-1,E11=-(E11(1)E11(2))×(a11+b11a12c21c22)-1,
with
(3.15)E20(1)=g1+g2q(2)(0)+g3(q(2)(0))2+g4q(1)(-1),E20(2)=g1′(q(1)(-τ2*τ1*))2+g2′q(1)(-τ2*τ1*)q(2)(0)+g3′q(1)(-τ2*τ1*)q(2)(-τ2*τ1*)+g4′q(2)(0)q(2)(-τ2*τ1*),E11(1)=2g1+g2(q(2)(0)+q-(2)(0))+2g3q(2)(0)q-(2)(0)+g4(q(1)(-1)+q-(1)(-1)),E11(2)=2g1′q(1)(-τ2*τ1*)q-(1)(-τ2*τ1*)+g2′(q(1)(-τ2*τ1*)q-(2)(0)+q-(1)(-τ2*τ1*)q(2)(0))+g3′(q(1)(-τ2*τ1*)q-(2)(-τ2*τ1*)+q-(1)(-τ2*τ1*)q(2)(-τ2*τ1*))+g4′(q(2)(0)q-(2)(-τ2*τ1*)+q-(2)(0)q(2)(-τ2*τ1*)).
Therefore, we can calculate the following values:
(3.16)C1(0)=i2τ1*ω*(g11g20-2|g11|2-|g02|23)+g212,μ2=-Re{C1(0)}Re{λ′(τ1*)},β2=2Re{C1(0)},T2=-Im{C1(0)}+μ2Im{λ′(τ1*)}τ1*ω*.
Based on the discussion above, we can obtain the following results.Theorem 3.1.
The direction of the Hopf bifurcation is determined by the sign ofμ2: if μ2>0(μ2<0), the Hopf bifurcation is supercritical (subcritical). The stability of bifurcating periodic solutions is determined by the sign of β2: if β2<0(β2>0), the bifurcating periodic solutions are stable (unstable). The period of the bifurcating periodic solutions is determined by the sign of T2: if T2>0(T2<0), the period of the bifurcating periodic solutions increases (decreases).


## 4. Numerical Example

In order to support the analytic results obtained above, we give some numerical simulations in this section. We only study the most important steady state, namely, the positive steady state. We consider the following system by taking the same coefficients as in [16]:
(4.1)dxdt=x(1-x(t-τ1))-xy0.01+3x+y,dydt=y[3.5-2y(t-τ2)x(t-τ2)],
where we get the positive equilibrium E*(0.6328,1.1074). For system (4.1), we can get that A+B+D=3.9017>0, namely, the condition (H1) holds.Forτ1=τ2=τ≠0. By a simple computation, we obtain that (2.18) has two positive roots: v1=12.0256, v2=0.3915. Thus, we know that the condition (H2) holds. Further, we get ω0=3.4678, τ0=0.4292. In addition, we have Λ1×Λ3+Λ2×Λ4=94.4826>0. Therefore, the condition (H3) is satisfied. Hence, from Theorem 2.1, we conclude that the positive equilibrium E*(0.6328,1.1074) is asymptotically stable when τ∈[0,τ0). The corresponding waveform and the phase plot are illustrated by Figure 1. When the delay τ passes through the critical value τ0 the positive equilibrium E*(0.6328,1.1074) will lose its stability and a Hopf bifurcation occurs, and a family of periodic solution bifurcates from the positive equilibrium E*(0.6328,1.1074). This property is illustrated by the numerical simulation in Figure 2.Figure 1
E
* is stable for τ1=τ2=τ=0.4<τ0=0.4292 with initial value 0.62,2.35.Figure 2
E
* is unstable for τ1=τ2=τ=0.44>τ0=0.4292 with initial value 0.62,2.35.Forτ1≠τ2, τ1>0 and τ2>0. Regard τ1 as a parameter and let τ2=0.36∈[0,τ20), and from (2.29) we can obtian τ20=0.4037. Then by a simple computation, we can obtain (2.34) has a positive root ω*=0.6792. Thus, the condition (H4) holds. Further, we we can get τ1*=2.3211 and Δ1×Δ3+Δ2×Δ4=3.7131>0. Namely, the condition (H5) holds. Thus, by Theorem 2.2, the positive equilibrium E*(0.6328,1.1074) is asymptotically stable when τ1∈[0,τ1*) and unstable when τ1>τ1*, and a family of periodic solution bifurcates from the positive equilibrium E*(0.6328,1.1074). The corresponding waveform and the phase plot are illustrated by Figures 3 and 4, respectively. In addition, from (3.16), we can get C1(0)=-26.2632-7.5707i, then μ2=12.1471, β2=-52.5264, T2=5.8586. From Theorem 3.1, we know that the Hopf bifurcation is supercritical and the bifurcating periodic solutions are stable.Figure 3
E
* is stable for τ1=0.8<τ1*=2.3211,τ2=0.36 with initial value 0.62,2.35.Figure 4
E
* is unstable for τ1=2.4>τ1*=2.3211,τ2=0.36 with initial value 0.62,2.35.


## 5. Conclusion

In this paper, a delayed predator-prey system with Beddington-DeAngelis functional response has been investigated. The bifurcation of a predator-prey system with single delay has been studied by many researchers [16, 27–30]. However, there are few papers considering the bifurcation of a predator-prey system with multiple delays (see [31–33]). Compared with the system considered in [16], the system in this paper accounts for not only the feedback delay of the prey density but also the feedback delay of the predator. The sufficient conditions for the stability of the positive equilibrium and the existence of periodic solutions via Hopf bifurcation at the positive equilibrium are obtained when τ1=τ2 and τ1≠τ2 with τ2∈[0,τ20). Special attention is paid to the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions. By computation, we find that the feedback delay of the predator is marked because the critical value of τ2 is much smaller when we only consider it. The feedback delay of the prey is unremarkable because the critical value of τ1 is much bigger when we consider it with τ2 in its stable interval. Furthermore, Zhang [16] has obtained that the two species in system (1.4) with only the feedback delay of the predator could coexist. However, we get that the two species could also coexist with some available feedback delays of the prey and the predator. This is valuable from the view of ecology.Unfortunately, the existence of the periodic solutions remain valid only in a small neighborhood of the critical value. It is definitely an interesting work to investigate whether these nonconstant periodic solutions which are obtained through local Hopf bifurcation can still exist for large values of the corresponding parameter time delay. The global continuation of the local Hopf bifurcation is left as the future work.



---

*Source: 236484-2012-11-08.xml*