# Delay-Induced Oscillations in a Competitor-Competitor-Mutualist Lotka-Volterra Model

Authors: Changjin Xu
Journal: Complexity
 (2017)

---

## Abstract

This paper deals with a competitor-competitor-mutualist Lotka-Volterra model. A series of sufficient criteria guaranteeing the stability and the occurrence of Hopf bifurcation for the model are obtained. Several concrete formulae determine the properties of bifurcating periodic solutions by applying the normal form theory and the center manifold principle. Computer simulations are given to support the theoretical predictions. At last, biological meaning and a conclusion are presented.


---

## Body

## 1. Introduction

The study of dynamical behaviors for tremendous predator-prey models has been a hot issue in population dynamics in the past few decades. Many results have been reported [1–11]. In the real world, any biological or environmental parameters are naturally subject to fluctuation in time. The effects of a periodically varying environment are important for evolutionary theory as the selective forces on systems in a fluctuating environment differ from those in a stable environment. Meanwhile, time delay due to gestation is common example because, generally, the consumption of prey by the predator throughout its past history governs the present birth rate of the predator. Based on all the above point, Lv et al. [12] had investigated the periodic solution of the following competitor-competitor-mutualist Lotka-Volterra model by using Krasnoselskii’s fixed point theorem(1)x˙1t=x1tr1t-a11tx1t-τ11t-a12tx2t-τ12t+a13tx3t-τ13t,x˙2t=x2tr2t-a21tx1t-τ21t-a22tx2t-τ22t+a23tx3t-τ23t,x˙3t=x3tr3t+a31tx1t-τ31t+a32tx2t-τ32t-a33tx3t-τ33t,where x1(t) and x2(t) denote the densities of competing species at time t and x3(t) denotes the density of cooperating species at time t. ri,aij∈C(R,[0,∞)) and τij∈C(R,R) are ω-periodic functions (ω>0). The parameters τij(t)≥0(i=1,2,3;j=1,2,3) are the feedback time delay of different species. In detail, one can see [12].It is well known that the research on the Hopf bifurcation, especially on the stability of bifurcating periodic solutions and direction of Hopf bifurcation, is one of the most important themes on the predator-prey dynamics. There are a great deal of papers which deal with this topic [11, 13–22]. The purpose of this paper is to discuss the stability and the properties of Hopf bifurcation of model (1). To simplify the analysis for model (1), we make the following assumptions: all biological and environmental parameters are constants in time and only the feedback time delay of competing species xi(i=1,2) to the growth of the species itself and the feedback time delay of cooperating species x3 to the growth of the species itself exist and are the same. Then system (1) can be described as the form(2)x˙1t=x1tr1-a11x1t-τ-a12x2t+a13x3t,x˙2t=x2tr2-a21x1t-a22x2t-τ+a23x3t,x˙3t=x3tr3+a31x1t+a32x2t-a33x3t-τ.In this paper, we consider the effect of time delayτ on the dynamics of system (2). We not only give the conditions on the stability of the positive equilibrium of (2) and the existence of periodic solutions but also derive the formulae for determining the properties of a Hopf bifurcation.The remainder of the paper is organized as follows. In Section2, we investigate the stability of the positive equilibrium and the occurrence of local Hopf bifurcations. In Section 3, the direction and stability of the local Hopf bifurcation are established. In Section 4, numerical simulations are carried out to illustrate the validity of the main results. Biological explanations and some main conclusions are drawn in Section 5.


## 2. Stability of the Positive Equilibrium and Local Hopf Bifurcations

Consider the realistic implication and actual application of biological system; in this section, we shall only study the stability of the positive equilibrium and the existence of local Hopf bifurcations. It is easy to see that system (2) has a unique positive equilibrium E0(x1∗,x2∗,x3∗) if the condition H1sign⁡Δ=sign⁡Δ1=sign⁡Δ2=sign⁡Δ3holds, where (3)Δ=det⁡a11a12-a13a21a22-a23a31a31-a33,Δ1=det⁡r1a12-a13r2a22-a23-r3a31-a33,Δ2=det⁡a11r1-a13a21r2-a23a31-r3-a33,Δ3=det⁡a11a12r1a21a22r2a31a31-r3.Let x¯1(t)=x1(t)-x1∗,  x¯2(t)=x2(t)-x2∗,  and x¯3(t)=x3(t)-x3∗ and still denote x¯i(t)(i=1,2,3) by xi(t)(i=1,2,3), and then (2) takes the form(4)x˙1t=b1x1t-τ+b2x2t+b3x3t+b4x1tx1t-τ+b5x1tx2t+b6x1tx3t,x˙2t=c1x1t+c2x2t-τ+c3x3+c4x1tx2t+c5x2tx2t-τ+c6x2tx3t,x˙3t=d1x1t+d2x2t+d3x3t-τ+d4x1tx3t+d5x2tx3t+d6x3tx3t-τ,where(5)b1=-a11x1∗,b2=-a12x1∗,b3=a13x1∗,b4=-a11,b5=-a12,b6=a13,c1=-a21x2∗,c2=-a22x2∗,c3=a23x2∗,c4=-a21,c5=-a22,c6=a23,d1=a31x3∗,d2=a32x3∗,d3=-a33x3∗,d4=a31,d5=a32,d6=-a33.The linearization of (6) near (0,0,0) is given by(6)x˙1t=b1x1t-τ+b2x2t+b3x3t,x˙2t=c1x1t+c2x2t-τ+c3x3,x˙3t=d1x1t+d2x2t+d3x3t-τ,whose characteristic equation takes the form(7)det⁡λ-b1e-λτ-b2-b3-c1λ-c2e-λτ-c3-d1-d2λ-d3e-λτ=0.That is,(8)λ3+m1λ+m2+m3λ2+m4e-λτ+m5λe-2λτ+m6e-3λτ=0,where (9)m1=b2c1+c3d2-b3d1,m2=-b1c3d1-b3d2c1,m3=-b1+c2+d3,m4=b3c2d1-b2c1d3-b1d2c3,m5=b2+c2d3+b1c2,m6=-b1c2d3.Multiplying eλτ on both sides of (8), it is easy to obtain(10)λ3+m1λ+m2eλτ+m3λ2+m4+m5λe-λτ+m6e-2λτ=0.We need the following lemma to discuss the stability of the positive equilibrium.Lemma 1 (see [23]).
For the transcendental equation(11)Pλ,e-λτ1,…,e-λτm=λn+p10λn-1+⋯+pn-10λ+pn0+p11λn-1+⋯+pn-11λ+pn1e-λτ1+⋯+p1mλn-1+⋯+pn-1mλ+pnme-λτm=0,as (τ1,τ2,τ3,…,τm) vary, the sum of orders of the zeros of P(λ,e-λτ1,…,e-λτm) in the open right half complex plane can change, only if a zero appears on or crosses the imaginary axis.Forτ=0, (10) becomes(12)λ3+m3λ2+m1+m5λ+m2+m4+m6=0.Obviously, m3>0. By the Routh-Hurwitz criteria, it follows that all eigenvalues of (12) have negative real parts if and only if the condition H2m2+m4+m6>0,m3m1+m5>m2+m4+m6is fulfilled.Forω>0,  iω is a root of (10) if and only if(13)-ω3i+im1ω+m2eiωτ+m3ω2i+m4+im5ωe-iωτ+m6e-2iωτ=0.Separating the real and imaginary parts, we get(14)m2cos⁡ωτ-m1ω-ω3sin⁡ωτ+m4+m5ωsin⁡ωτ=-m6cos⁡2ωτ,m2sin⁡ωτ+m1ω-ω3cos⁡ωτ+m3ω2+m5ωcos⁡ωτ=m6sin⁡2ωτ.It follows from (14) that(15)m2cos⁡ωτ-m1ω-ω3sin⁡ωτ+m4+m5ωsin⁡ωτ2+m2sin⁡ωτ+m1ω-ω3cos⁡ωτ+m3ω2+m5ωcos⁡ωτ2=m62.According to sin⁡ωτ=±1-cos2⁡ωτ, then (15) takes the form(16)m2cos⁡ωτ-m1ω-ω3±1-cos2⁡ωτ+m4+m5ω±1-cos2⁡ωτ2+m2±1-cos2⁡ωτ+m1ω-ω3cos⁡ωτ+m3ω2+m5ωcos⁡ωτ2=m62.It is easy to see that (16) is equivalent to(17)q1cos4⁡ωτ+q2cos3⁡ωτ+q3cos2⁡ωτ+q4cos⁡ωτ+q5=0,where (18)q1=k12+k52,q2=2k2k4+k3k5,q3=k22+k32-k52-2k4m62-k1,q4=2k2k1-m62-2k3k5,q5=m62-k12-k32,where (19)k1=m22+m1ω-ω32+m52ω2+m42+m32ω4,k2=2m2m4+2m2m5ω+2m3m5ω3+2m3ω2m1ω-ω3,k3=2m4m5ω-2m4m1ω-ω3+2m2m3ω2,k4=2m5ωm1ω-ω3,k5=2m2m5ω.Let cos⁡ωτ=r and denote (20)lr=r4+q2q1r3+q3q1r2+q4q1r+q5q1.It is easy to obtain that l′(r)=4r3+(3q2/q1)r2+(2q3/q1)r+q4/q1. Set(21)4r3+3q2q1r2+2q3q1r+q4q1=0.Let y=r+q2/4q1. Then (21) becomes(22)y3+γ1y+γ2=0,where γ1=q3/2q1-3q22/16q12 and  γ2=q23/32q13-q2q3/8q12+q4/4q1.Defineβ1=γ2/22+γ1/33,  β2=(-1+i3)/2. By (22), then we obtain (23)y1=-γ22+β13+-γ22-β13,y2=-γ22+β1β23+-γ22-β1β223,y3=-γ22+β1β223+-γ22-β1β23.By the discussion above, we can obtain the expression of cos⁡ωτ, say(24)cos⁡ωτ=f1ω,where f1(ω) is a function with respect to ω. Substitute (24) into (15); then we can easily get the expression of sin⁡ωτ, say(25)sin⁡ωτ=f2ω,where f2(ω) is a function with respect to ω. Thus we obtain(26)f12ω+f22ω=1.If all the coefficients of system (2) are given, it is easy to use computer to calculate the roots of (26) (say ω). Then from (24), we derive(27)τk=1ωarccos⁡f1ω+2kπk=0,1,2,….Letλ(τ)=α(τ)+iω(τ) be a root of (10) near τ=τk, α(τk)=0, and ω(τk)=ωk. Due to functional differential equation theory, for every τk,k=0,1,2,3,…, there exists ε>0 such that λ(τ) is continuously differentiable in τ for τ-τk<ε. Substituting λ(τ) into the left hand side of (10) and taking derivative with respect to τ, we have (28)dλdτ-1=-3λ2+m1e-λτ+2m3λ+m5e-λτλλ3+m1λ+m2eλτ-m5λ2e-λτ-2m6λe-2λτ-τλ.Then (29)dRe⁡λτdττ=τkj-1=-Re⁡3λ2+m1e-λτ+2m3λ+m5e-λτλλ3+m1λ+m2eλτ-m5λ2e-λτ-2m6λe-2λττ=τk=-Re⁡A3+iA4A1+iA2=A1A3+A2A4A12+A22,where (30)A1=ωk4-m1ωk2+m5ωk2cos⁡ωkτk-m2ωksin⁡ωkτk-2m6ωksin⁡2ωkτk,A2=ωk4-m1ωk2-m5ωk2sin⁡ωkτk+m2ωkcos⁡ωkτk-2m6ωkcos⁡2ωkτk,A3=m1+m5-3ωk2cos⁡ωkτk,A4=3ωk2-m1-m5sin⁡ωkτk+2m3ωk.In order to obtain the main results in this paper, it is necessary to make the following assumption: H3A1A3+A2A4≠0.In view of Lemma 1, it is easy to obtain the following result on stability and bifurcation of system (2).Theorem 2.
Suppose thatH1–H3 hold; then (i)
for system (2), its positive equilibrium E0(x1∗,x2∗,x3∗) is asymptotically stable for τ∈[0,τ0);(ii)
system (2) undergoes a Hopf bifurcation at the positive equilibrium E0(x1∗,x2∗,x3∗) when τ=τk(k=0,1,2,3,…); that is, system (2) has a branch of periodic solutions bifurcating from the positive equilibrium E0(x1∗,x2∗,x3∗) near τ=τk.


## 3. Direction and Stability of the Hopf Bifurcation

In this section, we shall analyze the direction, stability, and period of these periodic solutions bifurcating from the positive equilibriumE0(x1∗,x2∗,x3∗) at these critical value of τ, by using techniques from normal form and center manifold theory [24]. Throughout this section, we always assume that system (2) undergoes Hopf bifurcation at the positive equilibrium E0(x1∗,x2∗,x3∗) for τ=τk,  k=0,1,2,3,…, and then ±iωk are corresponding purely imaginary roots of the characteristic equation at the positive equilibrium E0(x1∗,x2∗,x3∗).For convenience, letx¯i(t)=xi(τt)(i=1,2,3) and τ=τk+μ, where τk is defined by (27) and μ∈R, drop the bar for the simplification of notations, and then system (4) can be written as an FDE in C=C([-1,0]),R3) as(31)u˙t=Lμut+Fμ,ut,where u(t)=(x1(t),x2(t),x3(t))T∈C and ut(θ)=u(t+θ)=(x1(t+θ),  x2(t+θ),  x3(t+θ))T∈C, and Lμ:C→R and F:R×C→R are given by(32)Lμϕ=τk+μ0b2b3c10c3d1d20ϕ10ϕ20ϕ30+τkj+μb1000c2000d3ϕ1-1ϕ2-1ϕ3-1,(33)fμ,ϕ=τk+μf1μ,ϕf2μ,ϕf3μ,ϕ,respectively, where ϕ(θ)=(ϕ1(θ),ϕ2(θ),ϕ3(θ))T∈C and (34)f1μ,ϕ=b4ϕ10ϕ1-1+b5ϕ10ϕ20+b6ϕ10ϕ30,f2μ,ϕ=c4ϕ10ϕ20+c5ϕ20ϕ2-1+c6ϕ20ϕ30,f3μ,ϕ=d4ϕ10ϕ30+d5ϕ20ϕ30+d6ϕ30ϕ3-1.From the discussion in Section 2, we know that if μ=0, then system (31) undergoes a Hopf bifurcation at the positive equilibrium E0(x1∗,x2∗,x3∗) and the associated characteristic equation of system (31) has a pair of simple imaginary roots ±ωkτk.By the representation theorem, there is a matrix function with bounded variation componentsη(θ,μ),  θ∈[-1,0], such that(35)Lμϕ=∫-10dηθ,μϕθ,forϕ∈C.In fact, we can choose(36)ηθ,μ=τk+μ0b2b3c10c3d1d20δθ-τk+μb1000c2000d3δθ+1,where δ is the Dirac delta function.Forϕ∈C([-1,0],R3), define(37)Aμϕ=dϕθdθ,-1≤θ<0,∫-10dηs,μϕs,θ=0,Rϕ=0,-1≤θ<0,fμ,ϕ,θ=0.Then (31) is equivalent to the abstract differential equation(38)u˙t=Aμut+Rμut,where ut(θ)=u(t+θ),  θ∈[-1,0].Forψ∈C([0,1],(R3)∗), define (39)A∗ψs=-dψsds,s∈0,1,∫-10dηTt,0ψ-t,s=0.Forϕ∈C([-1,0],R3) and ψ∈C([0,1],(R3)∗), define the bilinear form(40)ψ,ϕ=ψ¯0ϕ0-∫-10∫ξ=0θψTξ-θdηθϕξdξ,where η(θ)=η(θ,0), A=A(0), and A∗ are adjoint operators. By the discussions in Section 2, we know that ±iωkτk are eigenvalues of A(0), and they are also eigenvalues of A∗.Suppose thatq(θ)=(1,α,β)Teiωkτkθ is the eigenvector of A(0) corresponding to iωkτk; then A(0)q(θ)=iωkτkq(θ). It follows from definition A(0) and (32) and (35) that(41)τkiωk-b1e-iωkτk-b2-b3-c1iωk-c2e-iωkτk-c3-d1-d2iωk-d3e-iωkτk1αβ=000.Then we can obtain (42)α=c3b1e-iωkτk-iωk+b3c1b1c3-b3iωk-c2e-iωkτk,β=b1e-iωkτk-iωkiωk-c2e-iωkτkb1c3-b3iωk-c2e-iωkτk.Similarly, let q∗(s)=D(1,α∗,β∗)eiωkτks be the eigenvector of A∗ corresponding to -iωkτk(j); then A∗q∗(θ)=-iωkτkq∗(θ). It follows from the definition A∗ that (43)Dτkiωk+b1e-iωkτkc1d1b2iωk+c2e-iωkτkd2b3c3iωk+d3e-iωkτk1α∗β∗=000.Hence (44)α∗=b2c1-d2iωk+b1e-iωkτkc1d2-d1iωk+c2e-iωkτk,β∗=-iωk+c2e-iωkτkiωk+b1e-iωkτk-b2c1c1d2-d1iωk+c2e-iωkτk.In order to assure q∗(s),q(θ)=1 and q∗(s),q¯(θ)=0, we need to determine the value of D. From (40), we have(45)q∗s,qθ=D¯1,α∗¯,β∗¯1,α,βT-∫-10∫ξ=0θD¯1,α∗¯,β∗¯e-iωkτkξ-θdηθ1,α,βTeiωkτkξdξ=D¯1+αα∗¯+ββ∗¯-∫-101,α∗¯,β∗¯θeiωkτkθdηθ1,α,βT=D¯1+αα∗¯+ββ∗¯+τke-iωkτkb1+c2α∗¯α+d3β∗¯α.Thus we can choose (46)D=11+α¯α∗+β¯β∗+τkeiωkτkb1+c2α∗α¯+d3β∗β¯.Next, we use the same notations as those in Hassard et al. [24] and we first compute the coordinates to describe the center manifold C0 at μ=0. Let ut be the solution of (31) when μ=0.Define(47)zt=q∗,ut,Wt,θ=utθ-2Re⁡ztqθ,on the center manifold C0, and we have(48)Wt,θ=Wzt,z¯t,θ,where(49)Wzt,z¯t,θ=Wz,z¯=W20z22+W11zz¯+W02z¯22+⋯,and z and z¯ are local coordinates for center manifold C0 in the direction of q∗ and q¯∗. Noting that W is also real if ut is real, we consider only real solutions. For solutions ut∈C0 of (31)(50)z˙t=iωkτkz+q¯∗θf0,Wz,z¯,θ+2Re⁡zqθ=defiωkτkz+q¯∗0f0,which we write in abbreviated form as (51)z˙t=iωkτkz+gz,z¯,where (52)gz,z¯=g20z22+g11zz¯+g02z¯22+g21z2z¯2+⋯.Hence we have (53)gz,z¯=q¯∗0f0z,z¯=q¯∗0f0,ut=τkD¯1,α∗¯,β∗¯f1x1t,x2t,x3tf2x1t,x2t,x3tf3x1t,x2t,x3t=D¯τkH1+α∗¯H2+β∗¯H3z2+2D¯τkP1+α∗¯P2+β∗¯P3zz¯+D¯τkK1+α∗¯K2+β∗¯K3z¯2+D¯τkM1+α∗¯M2+β∗¯M3z2z¯+h.o.t.,where(54)f1x1t,x2t,x3t=b4x1t0x1t-1+b5x1t0x2t0+b6x1t0x3t0,f2x1t,x2t,x3t=c4x1t0x2t0+c5x2t0x2t-1+c6x2t0x3t0,f3x1t,x2t,x3t=d4x1t0x3t0+d5x2t0x3t0+d6x3t0x3t-1,H1=b4e-iωkτk+b5α+b6β,H2=c4α+c5αe-iωkτk+c6αβ,H3=d4β+d5αβ+d6β2e-iωkτk,P1=b4Re⁡eiωkτk+b5Re⁡α+b6Re⁡β,P2=c4Re⁡α+c5α2Re⁡α+c6Re⁡αβ¯,P3=d4Re⁡βeiωkτk+d5Re⁡αβ¯+d6β2Re⁡eiωkτk,K1=b4Reeiωkτk+b5α-+b6β-,K2=c4α¯+c5α¯2eiωkτk+c6Re⁡α¯β¯,K3=d4β¯+d5α¯β¯+d6β¯2eiωkτk,M1=b4W111-1+12W201-1+W1110e-iωkτk+12W2010eiωkτk+b512W2020+12W2010α¯+W1120+W1110+b6W1130+12W2030+12β¯W2010+W1110β,M2=c4W1120+12W2020+12α¯W2010+W1110+c5αW112-1+12α¯W202-1+αe-iωkτkW1120+12α¯W2020eiωkτk+c6αW1130+12α¯W2030+12β¯W2020+βW1120,M3=d4W1130+12W2030+12β¯W2010+βW2010+d5W1130α+12W2030α¯+12W2020β¯+βW1120+d6W113-1β+12W203-1β¯+βe-iωkτkW1130+α¯eiωkτkW1130.Thus we obtain (55)g20=2D¯τkH1+α∗¯H2+β∗¯H3,g11=2D¯τkP1+α∗¯P2+β∗¯P3,g02=2D¯τkK1+α∗¯K2+β∗¯K3,g21=2D¯τkkM1+α∗¯M2+β∗¯M3.For unknown W20(i)(0),  W11(i)(0),  W20(i)(-1),  W11(i)(-1),  (i=1,2,3) in g21, we still need to compute them.In view of (38) and (47), we have(56)W˙=AW-2Re⁡q¯∗0f¯qθ,-1≤θ<0,AW-2Re⁡q¯∗0f¯qθ+f¯,θ=0=defAW+Hz,z¯,θ,where(57)Hz,z¯,θ=H20θz22+H11θzz¯+H02θz¯22+⋯.Comparing the coefficients, we obtain(58)A-2iτkωkW20=-H20θ,(59)AW11θ=-H11θ.We know that for θ∈[-1,0)(60)Hz,z¯,θ=-q¯∗0f0qθ-q∗0f¯0q¯θ=-gz,z¯qθ-g¯z,z¯q¯θ.Comparing the coefficients of (60) with (57) gives that(61)H20θ=-g20qθ-g¯02q¯θ,(62)H11θ=-g11qθ-g¯11q¯θ.From (58), (61), and the definition of A, we get(63)W˙20θ=2iωkτkW20θ+g20qθ+g02¯q¯θ.Noting that q(θ)=q(0)eiωkτkθ, we have(64)W20θ=ig20ωkτkq0eiωkτkθ+ig¯023ωkτkq¯0e-iωkτkθ+E1e2iωkτkθ,where E1=(E1(1),E1(2),E1(3))T∈R3 is a constant vector.Similarly, from (59), (62), and the definition of A, we have(65)W˙11θ=g11qθ+g11¯q¯θ,(66)W11θ=-ig11ωkτkq0eiωkτkθ+ig¯11ωkτkq¯0e-iωkτkθ+E2,where E2=(E2(1),E2(2),E2(3))T∈R3 is a constant vector.In what follows, we shall seek appropriateE1 and  E2 in (64) and (66), respectively. It follows from the definition of A, (61), and (62) that(67)∫-10dηθW20θ=2iωkτkW200-H200,(68)∫-10dηθW11θ=-H110,where η(θ)=η(0,θ). From (58), we have(69)H200=-g20q0-g02¯q¯0+2τkH1,H2,H3T,(70)H110=-g11q0-g11¯0q¯0+2τkP1,P2,P3T.Noting that (71)iωkτkI-∫-10eiωkτkθdηθq0=0,-iωkτkI-∫-10e-iωkτkθdηθq¯0=0,and substituting (64) and (69) into (67), we have (72)2iωkτkI-∫-10e2iωkτkθdηθE1=2τkjH1,H2,H3T.That is, (73)2iωk-b1e-2iωkτk-b2-b3-c12iωk-c2e-2iωkτkj-c3-d1-d22iωk-d3e-2iωkτkE1=2H1H2H3.It follows that(74)E11=Δ11Δ1,E12=Δ12Δ1,E13=Δ13Δ1,where (75)Δ1=det⁡2iωk-b1e-2iωkτk-b2-b3-c12iωk-c2e-2iωkτk-c3-d1-d22iωk-d3e-2iωkτk,Δ11=2det⁡H1-b2-b3H22iωk-c2e-2iωkτk-c3H3-d22iωk-d3e-2iωkτk,Δ12=2det⁡2iωk-b1e-2iωkτkH1-b3-c1H2-c3-d1H32iωk-d3e-2iωkτk,Δ13=2det⁡2iωk-b1e-2iωkτk-b2H1-c12iωk-c2e-2iωkτkH2-d1-d2H3.Similarly, substituting (65) and (70) into (68), we have (76)∫-10dηθE2=2τkP1,P2,P3T.That is, (77)b1b2b3c1c2c3d1d2d3E2=2-P1-P2-P3.It follows that(78)E21=Δ21Δ2,E22=Δ22Δ2,E23=Δ23Δ2,where(79)Δ2=det⁡b1b2b3c1c2c3d1d2d3,Δ21=-det⁡P1b2b3P2c2c3P3d2d3,Δ22=-det⁡b1P1b3c1P2c3d1P3d3,Δ23=-det⁡b1b2P1c1c2P2d1d2P3.From (64), (66), (74), and (78), we can calculate g21 and derive the following values: (80)c10=i2ωkτkg20g11-2g112-g0223+g212,μ2=-Re⁡c10Re⁡λ˙τk,β2=2Re⁡c10,T2=-Im⁡c10+μ2Im⁡λ′τkωkτk.These formulae give a description of the Hopf bifurcation periodic solutions of (31) at τ=τk,  (k=0,1,2,3,…) on the center manifold. From the discussion above, we have the following result.Theorem 3.
The direction of the Hopf bifurcation is forward (backward) ifμ2>0(μ2<0); the bifurcating periodic solutions are orbitally asymptotically stable with asymptotical phase (unstable) if β2<0(β2>0); the periods of the bifurcating periodic solutions increase (decrease) if T2>0(T2<0).


## 4. Numerical Examples

In this section, numerical simulations are carried out to investigate the effect of the time delay on dynamical behaviors of system (2) as well as illustrate our theoretical results. For convenience, we assume that some parametric values of system (2) are kept as r1=7,  r2=6,  r3=5,  a11=3,  a12=2,  a13=0.2,  a21=2,  a22=3,  a23=0.2,  a31=0.3,  a32=0.3,  and a33=2. Then system (2) takes the form(81)x˙1t=x1t7-3x1t-τ-2x2t+0.2x3t,x˙2t=x2t6-2x1t-3x2t-τ+0.2x3t,x˙3t=x3t5+0.3x1t+0.3x2t-2x3t-τ,which has a positive equilibrium E0(1.9170,0.9170,2.9251) and all the conditions given in Theorem 2 are satisfied. When τ=0, the positive equilibrium E0(1.9170,0.9170,2.9251) is asymptotically stable. Fix k=0 and by means of Matlab 7.0, we get ω0≈0.9775,  τ0≈0.265, and  λ˙(τ0)≈0.6505–8.4322i. Thus we can derive the following values: c1(0)≈-2.3628–4.6636i,  μ2≈3.6323,  β2≈-4.7256,  and T2≈136.2424. Furthermore, it follows that μ2>0 and β2<0. Thus the positive equilibrium E0(1.9170,0.9170,2.9251) is stable when τ<τ0. When τ=0.26<τ0≈0.265, we can easily plot the time series of every population and phase portrait of the system and find that the solution of (81) with initial value x1(0)=1.9, x2(0)=0.85, and x3(0)=2.3 would tend to a positive equilibrium E0(1.9170,0.9170,2.9251) (see Figures 1(a)–1(j), in detail); that is, the positive equilibrium E0(1.9170,0.9170,2.9251) is asymptotically stable. When τ passes through the critical value τ0≈0.265, the positive equilibrium E0(1.9170,0.9170,2.9251) loses its stability and a Hopf bifurcation occurs; that is, a family of periodic solutions bifurcate from the positive equilibrium E0(1.9170,0.9170,2.9251). Figures 2(a)–2(j) are plotted by fixing the time delay τ=1 and initial value x1(0)=1.9,  x2(0)=0.85,  and x3(0)=2.3. It is shown that a Hopf bifurcation occurs from the positive equilibrium E0(1.9170,0.9170,2.9251). In view of μ2>0 and β2<0, we find that the direction of the Hopf bifurcation is τ>τ0≈0.265 and these bifurcating periodic solutions from E0(1.9170,0.9170,2.9251) at τ0≈0.265 are stable.Figure 1
Dynamical behavior of system (81) with τ=0.26<τ0≈0.265. The positive equilibrium E0(1.9170,0.9170,2.9251) is asymptotically stable. The initial value (x1(0),x2(0),x3(0))=(1.9,0.85,2.3).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)Figure 2
Dynamical behavior of system (81) with τ=0.28>τ0≈0.265. Hopf bifurcation occurs from the positive equilibrium E0(1.9170,0.9170,2.9251). The initial value (x1(0),x2(0),x3(0))=(1.9,0.85,2.3).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)


## 5. Biological Meanings and Conclusions

In this paper, the dynamics including the local stability of the positive equilibriumE0(x1∗,x2∗,x3∗) and the local Hopf bifurcation of a competitor-competitor-mutualist Lotka-Volterra model are investigated. It is showed that the positive equilibrium E0(x1∗,x2∗,x3∗) of system (2) is asymptotically stable for all τ∈[0,τ0) when the conditions H1–H3 are satisfied. This reveals that the densities of competing species x1,x2 and cooperating species x3 will tend to stabilization; in other words, the densities of competing species x1,x2 and cooperating species will tend to x1∗,x2∗,x3∗, respectively, and this situation does not vary with the delay τ∈[0,τ0). When the delay τ increases, the positive equilibrium E0(x1∗,x2∗,x3∗) loses its stability and a sequence of Hopf bifurcations occur from the positive equilibrium E0(x1∗,x2∗,x3∗); that is, a family of periodic orbits bifurcate from the positive equilibrium E0(x1∗,x2∗,x3∗). From this we can find that competing species x1,x2 and cooperating species x3 can coexist and keep in periodic mode; that is, the densities of competing species x1,x2 and cooperating species x3 will oscillate in the vicinity of x1∗,x2∗, and x3∗, respectively. Finally, the direction of Hopf bifurcation and the stability of the bifurcating periodic orbits are discussed by applying the normal form theory and the center manifold theorem. A numerical example is also included to verify our theoretical findings.



---

*Source: 2578043-2017-04-26.xml*