# Bifurcation Analysis in a Delayed Diffusive Leslie-Gower Model

Authors: Shuling Yan; Xinze Lian; Weiming Wang; Youbin Wang
Journal: Discrete Dynamics in Nature and Society
 (2013)

---

## Abstract

We investigate a modified delayed Leslie-Gower model under homogeneous Neumann boundary conditions. We give the stability analysis of the equilibria of the model and show the existence of Hopf bifurcation at the positive equilibrium under some conditions. Furthermore, we investigate the stability and direction of bifurcating periodic orbits by using normal form theorem and the center manifold theorem.


---

## Body

## 1. Introduction

The dynamic relationship between predators and their preys has long been and will continue to be one of dominant themes in both ecology and mathematical ecology due to its universal existence and importance. A major trend in theoretical work on prey-predator dynamics has been to derive more realistic models, trying to keep to maximum the unavoidable increase in complexity of their mathematics [1]. In this optic, recent years, the important Leslie-Gower predator-prey model [2, 3] has been extensively studied in [4–7]. A modified version of Leslie-Gower predator-prey model with Holling-type II functional response takes the form
(1)dHdt=H(a1-bH)-c1HPk1+H,dPdt=P(a2-c2Pk2+H),
where H and P represent prey and predator population densities at time t, respectively. a1,a2,b,c1,c2,k1, and k2 are positive constants. a1 is the growth rate of prey H. a2 describes the growth rate of predator P. b measures the strength of competition among individuals of species H. c1 is the maximum value of the per capita reduction of H due to P, and c2 is the maximum value of the per capita reduction of P due to H, which is not available in abundance. k1 measures the extent to which environment provides protection to prey H. k2 measures the extent to which environment provides protection to the predator P.On the other hand, time delay plays an important role in many biological dynamical systems, being particularly relevant in ecology [1]. For some predator-prey systems, the rate of the prey population depends on the predation of predator in the earlier times [8–14]. The results indicated that delay differential equations exhibit much more complicated dynamics than ordinary differential equations since a time delay could cause a stable equilibrium to become unstable and induce bifurcations.In this paper, we will focus on the complex dynamics of the delay effect in the extended reaction-diffusion model. The reproduction of the individuals is modeled by diffusion with diffusion coefficientsD1>0 and D2>0 for the prey and predator, respectively. This basic model is described by a system of two partial differential equations:
(2)∂H∂t=H(a1-bH)-c1HP(t-τ)k1+H+D1ΔH,x∈Ω,t>0,∂P∂t=P(a2-c2Pk2+H)+D2ΔP,x∈Ω,t>0,∂H∂n=∂P∂n=0,x∈∂Ω,t>0,
where H=H(t,x), P=P(t,x). Δ=∂2/∂x2, Ω is a bounded open domain in ℝ with boundary ∂Ω, n is the outward unit normal vector on ∂Ω, and homogeneous Neumann boundary conditions reflect the situation where the population cannot move across the boundary of the domain. And we incorporate a single discrete delay τ>0 in the negative feedback of the predator's density.The rest of the paper is organized as follows. In Section2, we give the stability property of the equilibria of model (1). In Section 3, we mainly analyze the distribution of the roots of the characteristic equation and show the occurrence of Hopf bifurcation at the positive equilibrium of model (2) under some conditions. In Section 4, we investigate the stability and direction of bifurcating periodic orbits by using normal form of theorem and the center manifold theorem, corresponding to theorems we also give some numerical simulations.


## 2. Equilibria Stability

In this section, we consider the existence and stability of the equilibria of model (1).It is easy to verify that model (1) always has three boundary equilibria:(i)
E1=(0,0) (extinction of prey and predator), which is a nodal source point;(ii)
E2=(a1/b,0) (extinction of the predator), which is a saddle point;(iii)
E3=(0,a2k2/c2) (extinction of the prey), which is a stable node when a2k2/c2>a1k1/c1.For the positive equilibria, we have(3)a1-bH-c1Pk1+H=0,a2-c2Pk2+H=0,
which yields
(4)bc2H2-(a1c2-a2c1-bc2k1)H+a2c1k2-a1c2k1=0.
For simplicity, we define
(5)A=bc2,B=a1c2-a2c1-bc2k1,C=a2c1k2-a1c2k1,
then (4) can be written as
(6)AH2-BH+C=0,
which has two roots given by
(7)h+=B+B2-4AC2A,h-=B-B2-4AC2A.Case 1 (C<0, i.e., k2<a1c2k1/a2c1).
Model (1) has a unique positive equilibrium E*=(h*,p*)=(h+,(a2/c2)(h++k2)).Case 2 (C=0, i.e., k2=a1c2k1/a2c1).
(i)
IfB>0, that is, k1<(a1c2-a2c1)/bc2, model (1) has a unique positive equilibrium E-=(h-,p-)=(h+,(a2/c2)(h++k2));
(ii)
IfB≤0, that is, k1≥(a1c2-a2c1)/bc2, (4) has no positive root; hence model (1) has no positive equilibrium.Case 3 (C>0, i.e., k2>a1c2k1/a2c1).
(i)
Suppose thatB>0, that is, k1<(a1c2-a2c1)/bc2, then
(i1)
ifB2-4AC>0, that is, k2<(a2c1-a1c2-bc2k1)2/4a2bc1c2+k1, model (1) has two positive equilibria E+=(h+,p+)=(h+,(a2/c2)(h++k2)) and E-=(h-,p-)=(h-,(a2/c2)(h-+k2));
(i2)
ifB2-4AC=0, that is, k2=(a2c1-a1c2-bc2k1)2/4a2bc1c2+k1, (4) has a unique positive root of multiplicity 2 given by he=B/2A=h+=h-, then model (1) has a unique positive equilibrium Ee=(he,pe)=(he,(a2/c2)(he+k2));
(i3)
ifB2-4AC=0, that is, k2>(a2c1-a1c2-bc2k1)2/4a2bc1c2+k1, (4) has no positive root; hence model (1) has no positive equilibrium;
(ii)
ifB<0, that is, k1>(a1c2-a2c1)/bc2, (4) has no positive root; hence model (1) has no positive equilibrium.We show the bifurcation diagram to display the distribute of the positive roots; in Figure1, the whole region has been divided into six parts; the number indicates the number of positive equilibria.Figure 1
The bifurcation diagram displays the distribute of the positive roots; the number indicates the number of positive equilibria.Next, we analyze the stability of these positive equilibria. LetE=(h,p) be arbitrary positive equilibrium, and the Jacobian matrix for E=(h,p) is given by
(8)J(E)=(-bh+c1hp(h+k1)2-c1hh+k1a22c2-a2).
The corresponding characteristic equation is
(9)λ2-tr(J(E))λ+det(J(E))=0,
where
(10)tr(J(E))=c1hp(h+k1)2-a2-bh=h(a1-a2-2bh)-(a2+bh)k1h+k1,det(J(E))=a2h(a2c1-a1c2+bc2k1+2bc2h)c2(h+k1).
And the sign of det(J(E)) is determined by
(11)F1(h)=a2c1-a1c2+bc2k1+2bc2h=2Ah-B.
Thus,
(12)F1(h*)=2Ah*-B=B2-4AC>0,F1(h-)=2Ah--B=B>0,F1(h+)=2Ah+-B=B2-4AC>0;F1(h-)=2Ah--B=-B2-4AC<0,F1(he)=2Ahe-B=0.
Hence, det(J(E*))>0, det(J(E-))>0, det(J(E+))>0, det(J(E-))<0, and det(J(Ee))=0. Obviously, the positive equilibrium E- is a saddle point.In the following, we study the stability of other positive equilibria. The sign oftr(J(E)) is determined by
(13)G1(h)=h(a1-a2-2bh)-(a2+bh)k1.
Then we can get
(14)G1(h*)=h*a2(c1-c2)c2-h*B2-4ACc2-a2k1,G1(h-)=h-a2(c1-c2)c2-h-Bc2-a2k1,G1(h+)=h+a2(c1-c2)c2-h+B2-4ACc2-a2k1,G1(he)=hea2(c1-c2)c2-a2k1.
Hence, if c1≤c2, then tr(J(E*))<0, tr(J(E-))<0, and tr(J(E+))<0 are true. Summarizing the above, we can obtain the following theorem.Theorem 1.
For model (1),(i)
ifk2<a1c2k1/a2c1 holds, the unique positive equilibrium E* is locally asymptotically stable for c1≤c2;(ii)
ifk2=a1c2k1/a2c1 and k1<(a1c2-a2c1)/bc2 hold, the unique positive equilibrium E- is locally asymptotically stable for c1≤c2;(iii)
ifa1c2k1/a2c1<k2<(a2c1-a1c2-bc2k1)2/4a2bc1c2+k1 and k1<(a1c2-a2c1)/bc2 hold, model (1) has two positive equilibria, the positive equilibrium E+ is locally asymptotically stable for c1≤c2, and E- is a saddle point.Figure2 shows the dynamics of model (1). In this case, E1 is a nodal source point; E2 is a saddle point; E3 is a nodal sink point, which is locally asymptotically stable; E+ is locally asymptotically stable; E- is a saddle point. There exists a separatrix curve determined by the stable manifold of E-, which divides the behavior of trajectories; that is, the stable manifold of saddle E- splits the feasible region into two parts such that orbits initiating inside tend to the positive equilibrium E+, while orbits initiating outside tend to E3 except for the stable manifolds of E-. This means that, in this situation, the trajectories of the model can have different behavior strongly depending on the initial conditions.Figure 2
Phase portraits of model (1) with the parameters a1=2, a2=0.5, b=1, c1=c2=2, k1=0.4, and k2=2.1. In this case, E1=(0,0) is a nodal source point; E2=(2,0) is a saddle point; E3=(0,0.525) is a nodal sink point, which is locally asymptotically stable; E+=(0.7791,0.7197) is a spiral sink, which is locally asymptotically stable; and E-=(0.3208,0.6052) is a saddle point. There exists a separatrix curve determined by the stable manifold of E-. The dashed curve is the H-nullcline H(a1-bH)-c1HP/(k1+H), and the dotted curve is the P-nullcline P(a2-c2P/(k2+H)).Theorem 2.
For model (1), if the unique positive equilibrium Ee=(he,(a2/c2)(he+k2)) exists, and k1=(a1c2-a2c1)(c1-c2)/bc2(c1+c2), then it is a cusp of codimension 2.Proof.
The Jacobian matrix atEe=(he,pe) is
(15)JEe=(-a2c1(a2c1-a1c2+bc2k1)c2(bc2k1-a2c1+a1c2)c1(a2c1-a1c2+bc2k1)bc2k1-a2c1+a1c2a22c2-a2),
we have know that det(J(Ee))=0. Moreover, tr(J(Ee))=0, if and only if
(16)k1=(a1c2-a2c1)(c1-c2)bc2(c1+c2).
Then
(17)JEe=(a2-c2a22c2-a2)=a2(1-c2a2a2c2-1),
and the associate Jordan matrix is
(18)J^=(0c1(a2c1-a1c2+bc2k1)bc2k1-a2c1+a1c200).
Hence, following [15, 16], we know that the unique positive equilibrium Ee is a cusp of codimension 2.


## 3. Stability and Hopf Bifurcation Analysis in Delayed Reaction-Diffusion Model (2)

According to the previous section, for model (1), we know that E1, E2, and E3 are unstable and E- is a saddle point, and note that a solution of the model (1) is also a solution of the model  (2), so they are also unstable for model (2). In the following, we will focus on the dynamics of the positive equilibria of model (2). As an example, we only give the proof of the unique positive equilibrium E* of model (2).Introducing small perturbationsH~=H-h*, and P~=P-p* and dropping the hats for simplicity of notation, then we have
(19)∂H∂t=(H+h*)(a1-b(H+h*))-c1(H+h*)(P(t-τ)+p*)k1+H+h*+D1ΔH,x∈Ω,t>0,∂P∂t=(P+p*)(a2-c2(P+p*)k2+H+h*)+D2ΔP,x∈Ω,t>0,∂H∂n=∂P∂n=0,x∈∂Ω,t>0.Denote(20)X={H,P∈W2,2(Ω):∂H∂n=∂P∂n=0,x∈∂Ω}.
In the abstract space C([-τ,0],X), model (19) can be regarded as the following abstract functional differential equation. 						(21)∂U∂t=DΔU+L(Ut)+F(Ut),
where U=(H,P)T, Ut=U(t+θ),θ∈[-τ,0], D=(D1D2), dom(Δ)⊂X, and L:C([-τ,0],X)↦X, F:C([-τ,0],X)↦X are given by
(22)L(ϕ)=((-bh*+c1h*p*(h*+k1)2)ϕ1(0)+-c1h*h*+k1ϕ2(-τ)a22c2ϕ1(0)-a2ϕ2(0)),F(ϕ)=((c1k1p*(k1+h*)3-b)ϕ12(0)-c1k1(k1+h*)2ϕ1(0)ϕ2(-τ)-c2p*2(k2+h*)3ϕ12(0)+2c2p*(k2+h*)2ϕ1(0)ϕ2(0)-c2k2+h*ϕ22(0)),
here ϕ=(ϕ1,ϕ2)T=Ut∈C([-τ,0],X). Then the linearization of model (19) near (h*,p*) is
(23)∂U∂t=DΔU+L(Ut).
Following [17], we obtain that the characteristic equation for liner model (23) is
(24)λy-DΔy-L(eλ·y)=0,y∈dom(Δ)⊂X,y≠0.
It is well known that the eigenvalue problem
(25)-Δψ=μψ,x∈Ω,∂ψ∂n=0,x∈∂Ω,
has eigenvalues μn={μni=-Din2,i=1,2,andn=0,1,2,…}, with the corresponding eigenfunctions ψn(x)=cosnx   (n∈ℕ={0,1,2,…}).Substitutingy=∑n=0∞‍cosnx(y1n,y2n)T into characteristic equation (24), we obtain
(26)(-bh*+c1h*p*(h*+k1)2-D1n2-c1h*h*+k1e-λτa22c2-a2-D2n2)(y1ny2n)=λ(y1ny2n).
Therefore the characteristic equation (24) is equivalent to
(27)λ2+Anλ+Bne-λτ+Cn=0,n∈ℕ,
where
(28)An=(D1+D2)n2+a2+bh*-c1h*p*(h*+k1)2,Bn=a22c1h*c2(h*+k1),Cn=D1D2n4+(D1+D2+a2+bh*-c1h*p*(h*+k1)2)n2+a2bh*-a2c1h*p*(h*+k1)2.The stability of the positive equilibriumE* can be determined by the distribution of the roots of (27); that is, the equilibrium E* is locally asymptotically stable if all the roots of (27) have negative real parts. From the result of [18], the sum of the multiplicities of the roots of (27) in the open right half plane changes only if a root appears on or crosses the imaginary axis. It can be verified that λ=0 is not a root of (27) for n∈ℕ.Theorem 3.
Ifk2<min{a1c2k1/a2c1,bc2(h*+k1)2/a2c1-2h*-k1} holds, then the unique positive equilibrium E* of model (2) is locally asymptotically stable.Proof.
Let±iω  (ω>0) be a pair of roots of (27); substituting iω into (27), then we have
(29)-ω2+iAnω+Bne-iωτ+Cn=0.
Separating the real part from image part, we have
(30)-ω2+Bncosωτ+Cn=0,Anω-Bnsinωτ=0,
then
(31)ω4+(An2-2Cn)ω2+Cn2-Bn2=0,
where
(32)An2-2Cn=(bh*-c1h*p*(h*+k1)2+D1n2)2+(a2+D2n2)2>0,Cn2-Bn2=D1D2n4+(D1+D2+a2+bh*-c1h*p*(h*+k1)2)n2-a2c1h*p*(h*+k1)2-a24c12(h*)2c22(h*+k1)2+a2bh*.
Obviously, if k2<bc2(h*+k1)2/a2c1-2h*-k1, Cn2-Bn2>0 is true. Thus (31) has no positive roots for all n∈ℕ. Hence, all the roots of (27) have negative real part. This completes the proof.If there exists an integern0∈ℕ such that for 0≤n≤n0, Cn2-Bn2<0, then (31) has a unique positive real root
(33)ω0n=-(An2-2Cn)+(An2-2Cn)2-4(Cn2-Bn2)2,
and (27) has a pair of pure imaginary roots ±iω0n, and
(34)τnj=τn0+2πjωn,j=0,1,2,…,0≤n≤n0,n0∈ℕ,
where τn0=arccos((ω0n2-Cn)/Bn)/ω0n.Letλ(τ)=γ(τ)+iω(τ) be the root of (24), where γ(τnj)=0 and ω(τnj)=ω0n when τ is close to τnj. Then we have the following transversality condition.Lemma 4.
For0≤n≤n0(n0∈ℕ), if(H1)
D1D2n4+(D1+D2+a2+bh*-c1h*p*/(h*+k1)2)n2<a2c1h*p*/(h*+k1)2 +a24c12(h*)2/c22(h*+k1)2-a2bh* holds, then (dγ/dτ)|τ=τnj>0 for j=0,1,2,….From this transversality condition, we know that whenτ passes through these critical values τnj, the sum of the multiplicities of the roots of (27) in the open right half plane will increases at least two.Summarizing the above results, we can obtain the following theorem.Theorem 5.
For0≤n≤n0(n0∈ℕ) if (H1) holds, the following statements are true:(i)
ifτ∈[0,τn0), then the equilibrium point E* is locally asymptotically stable;(ii)
ifτ>τn0, then the equilibrium E* is unstable;(iii)
τ=τnj(j=0,1,2,…) are Hopf bifurcation values of model (2).


## 4. Direction and Stability of Spatial Hopf Bifurcation

In the previous section, we have obtained the conditions under which model (2) undergoes a Hopf bifurcation at the the equilibrium point E* when τ crosses though the critical value τnj(0≤n≤n0,n0∈ℕ,j=0,1,2,…). In this section, we will study the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions by employing the center manifold theorem and normal form method [17, 19] for partial differential equations with delay. Then we compute the direction and stability of the Hopf bifurcation when τ~=τnj for fixed j∈{0,1,2,…}.Define(35)μ=-Re(c1(0))Re(λ′(τ~)),β2=2Re(c1(0)),T2=-Im(c1(0))+μ2Im(λ′(τ~))ω0nτ~,
where c1(0) is defined in the appendix. Then we can get the following theorem.Theorem 6.
For model (19), if (H1) holds, we have the following:(i)
μ determines the direction of the Hopf bifurcation: if μ>0 (μ<0), then the bifurcating periodic solutions exist for τ>τ~ (τ<τ~);(ii)
β2 determines the stability of bifurcating periodic solutions: the bifurcating periodic solutions are orbitally asymptotically stable (unstable) if β2<0 (β2>0);(iii)
T2 determines the period of the bifurcating periodic solutions: the period increases (decreases) if T2>0 (T2<0).The proof is deferred to the appendix.


## 5. Conclusions and Remarks

In this paper, we have considered a modified version of Leslie-Gower model with Holling-type II functional and delayed diffusive predator-prey model under homogeneous Neumann boundary conditions. The value of this study lies in two folds. First, it presents local asymptotic stability of the equilibria of model with and without delay and the existence of Hopf bifurcation, which indicates that the dynamics induced by time delay are rich and complex. Second, it give the analysis of direction and stability of spatial Hopf bifurcation, from which one can find that small sufficiently delays cannot change the stability of the positive equilibrium and large delays cannot only destabilize the positive equilibrium but also induce oscillatory behaviors near the positive equilibrium.In the following, we give some numerical examples to illustrate the dynamical behaviors of model (2). In Figure 3, τ=2<τ00=3.435144529, the unique positive equilibrium E*=(2,1.5) remains the stability; the population of the predator and the prey will tend to a steady state. However, in Figure 4, τ=4>τ00, the positive equilibrium E* losses its stability and Hopf bifurcation occurs, which means that a family of stable periodic solutions bifurcate from E* and the system goes into oscillations; it means that the predator coexists with the prey with oscillatory behaviors.The solutions of model (2) tends to aperiodically oscillatory orbit. The parameters are taken as D1=0.01, D2=0.1, a1=2, a2=1, b=0.5, c1=c2=2, k1=k2=1, and the initial values H(x,t)=2+0.1tsinx, P(x,t)=1.5+0.1tcosx. t∈[-2,0], x∈[0,π]. In this case, τ=4>τ00.
(a)
(b)The solutions of model (2) tends to aperiodically oscillatory orbit. The parameters are taken as D1=0.01, D2=0.1, a1=2, a2=1, b=0.5, c1=c2=2, k1=k2=1, and the initial values H(x,t)=2+0.1tsinx, P(x,t)=1.5+0.1tcosx. t∈[-4,0], x∈[0,π]. In this case, τ=4>τ00.
(a)
(b)Our results show that time-delay can make a stable equilibrium to become unstable and induce Hopf bifurcation and the system goes into oscillations; that’s to say, the dynamical behaviors of the delay reaction-diffusion equations are much more complex and rich than reaction-diffusion equations.



---

*Source: 170501-2013-02-24.xml*