# Further Study on Dynamics for a Fractional-Order Competitor-Competitor-Mutualist Lotka–Volterra System

Authors: Bingnan Tang
Journal: Complexity
 (2021)

---

## Abstract

On the basis of the previous publications, a new fractional-order prey-predator model is set up. Firstly, we discuss the existence, uniqueness, and nonnegativity for the involved fractional-order prey-predator model. Secondly, by analyzing the characteristic equation of the considered fractional-order Lotka–Volterra model and regarding the delay as bifurcation variable, we set up a new sufficient criterion to guarantee the stability behavior and the appearance of Hopf bifurcation for the addressed fractional-order Lotka–Volterra system. Thirdly, we perform the computer simulations with Matlab software to substantiate the rationalisation of the analysis conclusions. The obtained results play an important role in maintaining the balance of population in natural world.


---

## Body

## 1. Introduction

For a long time, the dynamic characteristics of interaction between predator population and prey population has been a central issue in ecology and biomathematics due to its general appearance and potential importance [1]. General speaking, the interaction between predator population and prey population includes four cases: competition, predation, mutualism, and parasitism [2]. During the past few decades, a great deal of valuable research fruits on dynamical behavior of the above four-type predator-prey models has been covered. For example, Alidousti and Ghafari [3] investigated the Hopf bifurcation and limit cycle of the fractional-order predator-prey model; Sasmal and Takeuchi [4] studied the stability behavior of all equilibria, bifurcation nature, global features, and multistability of a predator-prey model; Ryu and Ko [5] discussed the asymptotic peculiarity for positive solutions for a prey-predator model; Guo et al. [6] proved the appearance of traveling waves in a prey-predator system. Zhang et al. [7] revealed the effect of the fear factor on the periodic solution of a prey-predator model. As for concrete works, we refer the readers to [8–25].In real life, any biological or environmental coefficients will change with time. So, the parameters in predator-prey models are usually not fixed constants. They are often functions with respect to time. In particular, the influence of a periodically varying environment on the dynamics of predator-prey models plays a vital role in maintaining population balance. Furthermore, the capture of the prey from the predator throughout its past time has an important effect on the present birth rate of the predator [26, 27]. Thus, it is of importance to establish the various type predator-prey models with periodic coefficients. Based on this idea, in 2010, Lv et al. [27] established the predator-prey model involving periodic coefficients as follows:(1)u˙1t=u1tγ1t−α11tu1t−σ11t−α12tu2t−σ12t+α13tu3t−σ13t,u˙2t=u2tγ2t−α21tu1t−σ21t−α22tu2t−σ22t+α23tu3t−σ23t,u˙3t=u3tγ3t+α31tu1t−σ31t+α32tu2t−σ32t−α33tu3t−σ33t,where u1t and u2t stand for the densities of competing species at time t, u3t stands for the density of cooperating species at time t, And γl,αlk∈CR,0,∞ and σlk∈CR,R denote ψ-periodic function ψ>0. The parameter σlkt≥0l,k=1,2,3 is the delay. Applying the fixed point theory and constructing Lyapunov functions, Lv et al. [27] set up the condition to ensure the global stability of periodic solutions for model (1). For details, one can see [27].To reveal the Hopf bifurcation nature of model (1), Xu [26] assumes that all biological and environmental coefficients remain constants and only the feedback time delays of all species uii=1,2,3 to the growth of the species themselves exist and are same [26]. Then, model (1) is rewritten as the following form:(2)u˙1t=u1tγ1−α11u1t−σ−α12u2t+α13u3t,u˙2t=u2tγ2−α21u1t−α22u2t−σ+α23u3t,u˙3t=u3tγ3+α31u1t+α32u2t−α33u3t−σ.With the Hopf bifurcation theory, normal form theory, and center manifold principle, Xu [26] obtained a sufficient condition to guarantee the stability behavior and the appearance for bifurcation phenomenon of model (2). Also, they have derived the concrete expression to find the nature of bifurcation periodic solution.However, the mentioned works above are only restricted to the integer-order differential systems. Recently, the dynamics of fractional-order dynamical models has attracted great attention of many authors due to their extensive application in numerous areas, such as electromagnetic waves, medicine, mechanics, network science, biology, and finance [28–30]. Many scholars argue that fractional calculus is a powerful tool to depict real phenomena of the object world due to its owned hereditary and memory properties of various practical dynamical models [31]. In recent several decades, fractional calculus has attracted more and more attention from a large number of researchers in various fields. In particular, a great deal of interesting fruits on various dynamical natures of prey-predator systems has sprung up. For instance, Mondal et al. [32] derived the condition to ensure the stability behavior for a class of fractional-order prey-predator system; El-Saka et al. [33] analyzed the local stability and bifurcation of fractional-order predator-prey models; Li et al. [34] dealt with the dynamical property of the solutions and global asymptotic stability for a class of fractional-order prey-predator system. For more relational publications, one can see [35–37]. Here, we particularly emphasize that the study on Hopf bifurcation for fractional-order dynamical models starts relatively late. Up to now, only a few literatures have been published. For example, Alidousti [38] investigated the bifurcation phenomenon of a fractional-order predator-prey model; Huang et al. [30] proposed a novel control technique of bifurcation for a fractional-order prey-predator model with delays; Xu et al. [39] revealed the effect of two delays on bifurcation for fractional-order neural networks. Xiao et al. [40] put up PD control way for Hopf bifurcations of fractional-order networks. As for more related literatures, one can see [41–44].In terms of the above analysis, we think that it is meaningful for us to study the dynamics (especially Hopf bifurcation) of fractional-order prey-predator models. Based on the previous predator-prey model (2) and assuming that the feedback time delays of all species uii=1,2,3 to the growth of the species themselves and other species exists are same, then we propose the following fractional-order competitor-competitor-mutualist Lotka–Volterra system:(3)dϱu1tdtϱ=u1tγ1−α11u1t−σ−α12u2t−σ+α13u3t−σ,dϱu2tdtϱ=u2tγ2−α21u1t−σ−α22u2t−σ+α23u3t−σ,dϱu3tdtϱ=u3tγ3+α31u1t−σ+α32u2t−σ−α33u3t−σ,where u1t and u2t stand for the densities of competing population and u3t stands for the density of cooperating population, γl,αlk∈CR,0,∞l,k=1,2,3, the parameter σ≥0 is the feedback time delay of different species, and ϱ∈0,1 is a constant. For more concrete meaning of coefficients for system (3), see [26, 27].The initial condition of system (3) takes the form(4)u1ϑ=u1ϑ,−σ≤ϑ≤t0,u2ϑ=u2ϑ,−σ≤ϑ≤t0,u3ϑ=u3ϑ,−σ≤ϑ≤t0,where t0>0 is a constant. The key object of this work focuses on existence, uniqueness, nonnegativity, stability, and bifurcation phenomenon of model (3). Different from the methodology in [26, 27], in this paper, we will mainly discuss the various dynamics by applying fractional-order differential equation theory. Due to the introduction of fractional order, various dynamical behaviors of the predator-prey model (3) are different from the integer-order case. We think that it is necessary to reveal the effect of time delay, parameters of system, and fractional order on dynamics such as the stability and Hopt bifurcation. Thus, this study has significance in theory and practice.We plan the structure of this manuscript as follows. Section2 gives some related knowledge about fractional-order dynamical systems. Section 3 discusses the existence, uniqueness, nonnegativity, local stability, and Hopf bifurcation of model (3). Section 4 gives an example to support the effectiveness of the obtained key conclusions. Section 6 ends our work.


## 2. Basic Knowledge

In this part, we present several related definitions and lemmas on fractional-order dynamical systems that will be applied in the later proof. LetR+ represent the set of all nonnegative real numbers.Definition 1.
(see [45]). Define Caputo fractional-order derivative as follows:(5)Dϱhw=1Γn−ϱ∫w0whlsw−sw−l+1ds,where hw∈w0,∞,R,Γs=∫0∞ws−1e−wdw, w≥w0, and l∈Z+,l−1≤ϱ<l.
The Laplace transform of Caputo fractional-order derivative is defined as follows:(6)ℒDϱht;s=sϱℋs−∑j=0m−1sϱ−j−1hj0,m−1≤ϱ<m∈Z+,where ℋs=ℒht. Especially, if hj0=0,j=1,2,…,m, then ℒℋϱht;s=sϱℋs.Definition 2.
(see [46]). v1∗,v2∗,v3∗ is called an equilibrium point of model (3) provided that the following system(7)v1∗γ1−α11v1∗−α12v2∗+α13v3∗=0,v2∗γ2−α21v1∗−α22v2∗+α23v3∗=0,v3∗γ3+α31v1∗+α32v2∗−α33v3∗=0,holds.Lemma 1.
(see [47]). Given the following system(8)dϱvtdtϱ=ht,vt,vt0=vt0,t0>0,where ϱ∈0,1,h:t0,∞×Λ⟶R+n,Λ⊂R+n. If ht,v satisfies the local Lipschitz condition with respect to v, then model (3) possesses a unique solution defined on t0,∞.Lemma 2.
(see [48]). Suppose that ϱ∈0,1,ht∈Cγ,δ and Dϱht∈Cγ,δ, where γ,δ∈R. If Dϱht≥0,t∈γ,δ, then ht is a nondecreasing function ∀t∈γ,δ. If Dϱht≤0,t∈γ,δ, then ht is a nonincreasing function ∀t∈γ,δ.Lemma 3.
(see [49, 50]). For a given fractional-order model,(9)dϱvtdtϱ=gt,vt,v0=v0,where ϱ∈0,1 and gt,vt:R+×Rm⟶Rm, the equilibrium point of model (9) is locally asymptotically stable if all eigenvalues ς of ∂gt,v/∂v evaluated near the equilibrium point satisfy argς>ϱπ/2.Lemma 4.
(see [51]). For given $m$-dimensional fractional-order system,(10)dϱ1V1tdtϱ1=l11V1t−σ11+l12V2t−σ12+⋯+l1mVmt−σ1m,dϱ2V2tdtϱ2=l21V1t−σ21+l22V2t−σ22+⋯+l2mVmt−σ2m,⋮dϱmVmtdtϱm=lm1V1t−σm1+lm2V2t−σm2+⋯+lmmVmt−σmm,where 0<ϱi<1i=1,2,…,m, the initial values Vit=φit∈C−maxi,lσil,0, and t∈−maxi,lσil,0, i,l=1,2,…,m. Denote(11)Δχ=χϱ1−l11e−χσ11−l12e−χσ12⋯−l1me−χσ1m−l21e−χσ12χϱ2−l22e−χσ22⋯−l2me−χσ2m⋮⋮⋱⋮−lm1e−χσm1−lm2e−χσm2⋯χϱm−lmme−χσmm,and then, the zero solution of equation (10) is Lyapunov asymptotically stable if all roots of detΔχ=0 possess negative real parts.


## 3. Main Results

### 3.1. Existence and Uniqueness

Theorem 1.
LetΘ=u1,u2,u3∈R3:maxu1,u2,u3<Q, where Q is a positive constant. Then, for every Ut0=u1t0,u2t0,u3t0∈Θ and for every t≥t0, there exists a unique solution Ut∈Θ of system (3) with initial value Ut0.Proof.
We discuss this issue for model (3) in Θ×t0,t∗, where t∗<+∞. Let U=u1,u2,u3 and U˜=u˜1,u˜2,u˜3, and define the following mapping: LU=L1U,L2U,L3U, where(12)L1U=u1tγ1−α11u1t−σ−α12u2t−σ+α13u3t−σ,L2U=u2tγ2−α21u1t−σ−α22u2t−σ+α23u3t−σ,L3U=u3tγ3+α31u1t−σ+α32u2t−σ−α33u3t−σ.
For arbitraryU,U˜∈Θ,t1,t2∈R,t1,t2≥t0, one has(13)LU−LU˜=L1U−L1U˜+L2U−L2U˜+L3U−L3U˜=|u1t1γ1−α11u1t1−σ−α12u2t1−σ+α13u3t1−σ−u1t2γ1−α11u1t2−σ−α12u2t2−σ+α13u3t2−σ|+|u2t1γ2−α21u1t1−σ−α22u2t1−σ+α23u3t1−σ−u2t2γ2−α21u1t2−σ−α22u2t2−σ+α23u3t2−σ|+|u3t1γ3+α31u1t1−σ+α32u2t1−σ−α33u3t1−σ−u3t2γ3+α31u1t2−σ+α32u2t2−σ−α33u3t2−σ|≤u1t1−u1t2γ1−α11u1t1−σ−α12u2t1−σ+α13u3t1−σ+maxu1t1,u1t2|γ1−α11u1t1−σ−α12u2t1−σ+α13u3t1−σ−γ1−α11u1t2−σ−α12u2t2−σ+α13u3t2−σ|+u2t1−u2t2γ2−α21u1t1−σ−α22u2t1−σ+α23u3t1−σ+maxu2t1,u2t2|γ2−α21u1t1−σ−α22u2t1−σ+α23u3t1−σ−γ2−α21u1t2−σ−α22u2t2−σ+α23u3t2−σ|+u3t1−u3t2γ3+α31u1t1−σ+α32u2t1−σ−α33u3t1−σ+maxu3t1,u3t2|γ3+α31u1t1−σ+α32u2t1−σ−α33u3t1−σ−γ3+α31u1t2−σ+α32u2t2−σ−α33u3t2−σ|≤u1t1−u1t2γ1+α11Q+α12Q+α13Q+Qα11u1t1−σ−u1t2−σ+Qα12u2t1−σ−u2t2−σ+Qα13u3t1−σ−u3t2−σ+u2t1−u2t2γ2+α21Q+α22Q+α23Q+Qα21u1t1−σ−u1t2−σ+Qα22u2t1−σ−u2t2−σ+Qα23u3t1−σ−u3t2−σ+u3t1−u3t2γ3+α31Q+α32Q+α33Q+Qα31u1t1−σ−u1t2−σ+Qα32u2t1−σ−u2t2−σ+Qα33u3t1−σ−u3t2−σ≤L1u1t1−u1t2+L2u2t1−u2t2+L3u3t1−u3t2,where(14)L1=γ1+Q2α11+α21+α31+α12+α13,L2=γ2+Q2α21+α22+α32+α22+α23,L3=γ3+Q2α31+α23+α33+α32+α33.
It follows from (13) that(15)LU−LU˜≤LU−U˜|,where(16)L=maxL1,L2,L3.
In view of Lemma1, one can conclude that Theorem 1 is true. We finish the proof.


### 3.2. Nonnegativity

DenoteΘ+=u1,u2,u3∈Θ:u1,u2,u3∈R+.Theorem 2.
Each solution of model (3) that begins with Θ+ remains nonnegative.Proof.
Our object is to prove thatuit≥0i=1,2,3∀t≥t0. Let Ut0=u1t0,u1t0,u1t0∈Θ+ be the initial value of system (3). Assume that ∃ is a constant and t0≤t<ρ such that(17)uit=0,t0≤t<ρ,uiρ=0,uiρ+<0,where i=1,2,3. In view of system (3), one has(18)dϱu1tdtϱ|ui=0=0,dϱu2tdtϱ|ui=0=0,dϱu3tdtϱ|ui=0=0,where i=1,2,3. By Lemma 2, one knows that uiρ+=0i=1,2,3, which contradicts uiρ+<0i=1,2,3 in (17). Thus, one can conclude that uit≥0∀t≥t0. The proof is complete.


### 3.3. Stability and Hopf Bifurcation

In the section, we shall focus on the local stability and the appearance of Hopf bifurcation of system (3). Consider the biological implication of system (3); we only seek the sufficient condition to guarantee the local stability of the positive equilibrium and the emergence of Hopf bifurcations of system (3).Clearly, system (3) has a unique positive equilibrium U0u1∗,u2∗,u3∗ provided that the following condition(19)ℳ1signA=signA1=signA2=signA3,is satisfied, where(20)A=detα11α12−α13α21α22−α23α31α31−α33≠0,A1=detγ1α12−α13γ2α22−α23−γ3α31−α33,A2=detα11γ1−α13α21γ2−α23α31−γ3−α33,A3=detα11α12γ1α21α22γ2α31α31−γ3,u1∗=A1A,u2∗=A2A,u3∗=A3A.The linear system of equation (3) near U0u1∗,u2∗,u3∗ takes(21)dϱu1tdtϱ=a11u1t−σ+a12u2t−σ+a13u3t−σ,dϱu2tdtϱ=a21u1t−σ+a22u2t−σ+a23u3t−σ,dϱu3tdtϱ=a31u1t−σ+a32u2t−σ+a33u3t−σ,where(22)a11=−u1∗α11,a12=−u1∗α12,a13=u1∗α13,a21=−u2∗α21,a22=−u2∗α22,a23=u2∗α23,a31=u3∗α31,a32=u3∗α32,a33=−u3∗α33.The characteristic equation of equation (21) is computed as follows:(23)detsϱ−a11e−sσ−a12e−sσ−a13e−sσ−a21e−sσsϱ−a22e−sσ−a23e−sσ−a31e−sσ−a32e−sσsϱ−a33e−sσ=0,which leads to(24)s3ϱ+b1s2ϱe−sσ+b2sϱe−2sσ+b3e−3sσ=0,where(25)b1=−a11+a22+a33,b2=a11a33+a22a33+a11a22−a13a31−a12a21−a23a32,b3=a13a22a31+a12a21a33+a11a23a32−a11a22a33−a13a21a32.In terms of (24), we have(26)s3ϱe2sσ+b1s2ϱesσ+b2sϱ+b3e−sσ=0.Lets=iφ=φcosπ/2+isinπ/2 be the root of (26); then, one has(27)G1ϱcos2φσ+G2ϱsin2φσ+G3ϱcosφσ+G4ϱsinφσ=0,S1ϱcos2φσ+S2ϱsin2φσ+S3ϱcosφσ+S4ϱsinφσ=0,where(28)G1ϱ=φ3ϱcos3ϱπ2+b2φϱcosϱπ2,G2ϱ=b2φϱsinϱπ2−φ3ϱsin3ϱπ2,G3ϱ=b1φ2ϱcosϱπ+b3,G4ϱ=−b1φ2ϱsinϱπ,S1ϱ=φ3ϱsin3ϱπ2+b2φϱsinϱπ2,S2ϱ=φ3ϱcos3ϱπ2−b2φϱcosϱπ2,S3ϱ=b1φ2ϱsinϱπ,S4ϱ=b1φ2ϱcosϱπ−b3.In (28), let(29)c1=cos3ϱπ2,c2=b2cosϱπ2,c3=b2sinϱπ2,c4=−sin3ϱπ2,c5=b1cosϱπ,c6=b3,c7=−b1sinϱπ,c8=sin3ϱπ2,c9=b2sinϱπ2,c10=cos3ϱπ2,c11=−b2cosϱπ2,c12=b1sinϱπ,c13=b1cosϱπ,c14=−b3,and then, (28) takes the form:(30)G1ϱ=c1φ3ϱ+c2φϱ,G2ϱ=c3φϱ+c4φ3ϱ,G3ϱ=c5φ2ϱ+c6,G4ϱ=c7φ2ϱ,S1ϱ=c8φ3ϱ+c9φϱ,S2ϱ=c10φ3ϱ+c11φϱ,S3ϱ=c12φ2ϱ,S4ϱ=c13φ2ϱ+c14.For (27), we consider two cases.Case 1.
Ifsinφσ=1−cos2φσ, then the first equation of (27) becomes(31)G1ϱ2cos2φσ−1+2G2ϱcosφσ1−cos2φσ+G3ϱcosφσ+G4ϱ1−cos2φσ=0.
Hence,(32)2G2ϱcosφσ1−cos2φσ+G4ϱ1−cos2φσ2=G1ϱ2cos2φσ−1+G3ϱcosφσ2,which leads to(33)ϑ1cos4φσ+ϑ2cos3φσ+ϑ3cos2φσ+ϑ4cosφσ+ϑ5=0,where(34)ϑ1=4G12ϱ+4G22ϱ,ϑ2=4G1ϱG3ϱ+4G2ϱG4ϱ,ϑ3=G32ϱ+G42ϱ−4G12ϱ−4G22ϱ,ϑ4=−2G1ϱG3ϱ+4G2ϱG4ϱ,ϑ5=G12ϱ−G42ϱ.
Setcosφσ=r, and let(35)zr=r4+ϑ2ϑ1r3+ϑ3ϑ1r2+ϑ4ϑ1r+ϑ5ϑ1.
Hence,(36)dzrdr=4r3+3ϑ2ϑ1r2+2ϑ3ϑ1r+ϑ4ϑ1.
Denote(37)r3+3ϑ2ϑ1r2+2ϑ3ϑ1r+ϑ4ϑ1=0.
Setp=r+ϑ2/4ϑ1, and then, (37) takes the form(38)p3+f1p+f2=0,where(39)f1=ϑ32ϑ1−3ϑ2216ϑ12,f2=ϑ2332ϑ13−ϑ2ϑ38ϑ12+ϑ44ϑ1.
Let(40)q1=f222+f133,q2=−1+i32.
It follows from (38) that(41)p1=−f22+q13+−f22−q13,p2=−f22+q1q23+−f22−q1q223,p3=−f22+q1q223+−f22−q1q23.
According to the discussion above, one easily obtains the expression ofcosφσ. Then, one easily gets the expression of sinφσ. Here, we suppose that(42)cosφσ=g1φ,sinφσ=g2φ.
Then,(43)g12φ+g22φ=1.
With the aid of Matlab 7.0, one can get the root (sayφ) of equation (43). Then, one obtains(44)σ1k=1φarccosg1φ+2kπ,k=0,1,2,…Case 2.
Ifsinφσ=−1−cos2φσ, applying the same way to this case, we can obtain(45)cosφσ=h1φ,sinφσ=h2φ.
Then,(46)h12φ+h22φ=1.
By Matlab 7.0, one can obtain the root (sayφ) of equation (46). So,(47)σ2k=1φarccosh1φ+2kπ,k=0,1,2,…
Set(48)σ0=minσ1k,σ2k,k=0,1,2,…,
In the sequel, we check the transversality condition of the appearance of Hopf bifurcation. We give the hypothesis as follows:ℳ2Z11Z21+Z12Z22>0, where(49)Z11=3ϱφ03ϱ−1cos3ϱ−1π2cos2φ0σ0−sin3ϱ−1π2sin2φ0σ0+2φ03ϱσ0cos3ϱπ2cos2φ0σ0−sin3ϱπ2sin2φ0σ0+2b1ϱφ02ϱ−1cos2ϱ−1π2cosφ0σ0−sin2ϱ−1π2sinφ0σ0+b1σ0φ02ϱcosϱπcosφ0σ0−sinϱπsinφ0σ0+b2ϱφ0ϱ−1cosϱ−1π2−b3σ0cosφ0σ0,Z12=3ϱφ03ϱ−1cos3ϱ−1π2sin2φ0σ0+sin3ϱ−1π2cos2φ0σ0+2φ03ϱσ0cos3ϱπ2sin2φ0σ0+sin3ϱπ2cos2φ0σ0+2b1ϱφ02ϱ−1cos2ϱ−1π2sinφ0σ0+sin2ϱ−1π2cosφ0σ0+b1σ0φ02ϱcosϱπsinφ0σ0+sinϱπcosφ0σ0+b2ϱφ0ϱ−1sinϱ−1π2+b3σ0sinφ0σ0,Z21=b3φ0sinφ0σ0−2φ03ϱ+1cos3ϱ+1π2cos2φ0σ0−sin3ϱ+1π2sin2φ0σ0−2b2ϱ02σsinσπ+b1ϱ0σsinσπ2ϱ0cos2ϱ0η0,−b1φ02ϱ+1cos2ϱ+1π2cos2φ0σ0−sin2ϱ+1π2sin2φ0σ0,Z22=b3φ0cosφ0σ0−2φ03ϱ+1cos3ϱ+1π2sin2φ0σ0+sin3ϱ+1π2cos2φ0σ0−b1φ02ϱ+1cos2ϱ+1π2sin2φ0σ0+sin2ϱ+1π2cos2φ0σ0.Lemma 5.
Suppose thatsσ=ξσ+iησ is a root of equation (26) near σ=σ0 such that ξσ0=0,ησ0=φ0; then, one has Reds/dσ|σ=σ0,φ=φ0>0.Proof.
In terms of equation (26), one knows that(50)3ϱs3ϱ−1e2sσdsdσ+2s3ϱe2sσdsdσσ+s+2b1ϱs2ϱ−1esσdsdσ+b1esσs2ϱdsdσσ+s+b2ϱsϱ−1dsdσ−b3e−sσdsdσσ+s=0,which leads to(51)3ϱs3ϱ−1e2sσ+2s3ϱe2sσσ+2b1ϱs2ϱ−1esσ+b1esσs2ϱσ+b2ϱsϱ−1−b3e−sσσdsdσ=b3se−sσ−2s3ϱ+1e2sσ−b1s2ϱ+1esσ.
Then,(52)dsdσ−1=Z1sZ2s−σs,where(53)Z1s=3ϱs3ϱ−1e2sσ+2s3ϱe2sσσ+2b1ϱs2ϱ−1esσ+b1esσs2ϱσ+b2ϱsϱ−1−b3e−sσσ,Z2s=b3se−sσ−2s3ϱ+1e2sσ−b1s2ϱ+1esσ.
So,(54)Redsdσ|σ=σ0,φ=φ0=ReZ1sZ2s|σ=σ0,φ=φ0=Z11Z21+Z12Z22Z212+Z222.
Byℳ2, one obtains(55)Redsdσ−1|σ=σ0,φ=φ0>0.
This finishes the proof of Lemma5.
Suppose that(56)ℳ3b1>0,b1b2>b3,b3>0.Lemma 6.
Provided thatσ=0 and ℳ3 is fulfilled, then system (3) is locally asymptotically stable.Proof.
Ifσ=0, then (24) takes the form(57)λ3+b1λ2+b2λ+b3=0.
Byℳ3, one knows that all roots λi of (57) satisfy argλi>ϱπ/2i=1,2,3,4. Thus, one knows that Lemma 6 is right. The proof finishes.
Based on the investigation above, we have the following conclusion.Theorem 3.
Provided thatℳ1–ℳ3 hold true, then U0u1∗,u1∗,u3∗ of system (3) is locally asymptotically stable if σ∈0,σ0, and a Hopf bifurcation exists near U0u1∗,u1∗,u3∗ if σ=σ0.


## 3.1. Existence and Uniqueness

Theorem 1.
LetΘ=u1,u2,u3∈R3:maxu1,u2,u3<Q, where Q is a positive constant. Then, for every Ut0=u1t0,u2t0,u3t0∈Θ and for every t≥t0, there exists a unique solution Ut∈Θ of system (3) with initial value Ut0.Proof.
We discuss this issue for model (3) in Θ×t0,t∗, where t∗<+∞. Let U=u1,u2,u3 and U˜=u˜1,u˜2,u˜3, and define the following mapping: LU=L1U,L2U,L3U, where(12)L1U=u1tγ1−α11u1t−σ−α12u2t−σ+α13u3t−σ,L2U=u2tγ2−α21u1t−σ−α22u2t−σ+α23u3t−σ,L3U=u3tγ3+α31u1t−σ+α32u2t−σ−α33u3t−σ.
For arbitraryU,U˜∈Θ,t1,t2∈R,t1,t2≥t0, one has(13)LU−LU˜=L1U−L1U˜+L2U−L2U˜+L3U−L3U˜=|u1t1γ1−α11u1t1−σ−α12u2t1−σ+α13u3t1−σ−u1t2γ1−α11u1t2−σ−α12u2t2−σ+α13u3t2−σ|+|u2t1γ2−α21u1t1−σ−α22u2t1−σ+α23u3t1−σ−u2t2γ2−α21u1t2−σ−α22u2t2−σ+α23u3t2−σ|+|u3t1γ3+α31u1t1−σ+α32u2t1−σ−α33u3t1−σ−u3t2γ3+α31u1t2−σ+α32u2t2−σ−α33u3t2−σ|≤u1t1−u1t2γ1−α11u1t1−σ−α12u2t1−σ+α13u3t1−σ+maxu1t1,u1t2|γ1−α11u1t1−σ−α12u2t1−σ+α13u3t1−σ−γ1−α11u1t2−σ−α12u2t2−σ+α13u3t2−σ|+u2t1−u2t2γ2−α21u1t1−σ−α22u2t1−σ+α23u3t1−σ+maxu2t1,u2t2|γ2−α21u1t1−σ−α22u2t1−σ+α23u3t1−σ−γ2−α21u1t2−σ−α22u2t2−σ+α23u3t2−σ|+u3t1−u3t2γ3+α31u1t1−σ+α32u2t1−σ−α33u3t1−σ+maxu3t1,u3t2|γ3+α31u1t1−σ+α32u2t1−σ−α33u3t1−σ−γ3+α31u1t2−σ+α32u2t2−σ−α33u3t2−σ|≤u1t1−u1t2γ1+α11Q+α12Q+α13Q+Qα11u1t1−σ−u1t2−σ+Qα12u2t1−σ−u2t2−σ+Qα13u3t1−σ−u3t2−σ+u2t1−u2t2γ2+α21Q+α22Q+α23Q+Qα21u1t1−σ−u1t2−σ+Qα22u2t1−σ−u2t2−σ+Qα23u3t1−σ−u3t2−σ+u3t1−u3t2γ3+α31Q+α32Q+α33Q+Qα31u1t1−σ−u1t2−σ+Qα32u2t1−σ−u2t2−σ+Qα33u3t1−σ−u3t2−σ≤L1u1t1−u1t2+L2u2t1−u2t2+L3u3t1−u3t2,where(14)L1=γ1+Q2α11+α21+α31+α12+α13,L2=γ2+Q2α21+α22+α32+α22+α23,L3=γ3+Q2α31+α23+α33+α32+α33.
It follows from (13) that(15)LU−LU˜≤LU−U˜|,where(16)L=maxL1,L2,L3.
In view of Lemma1, one can conclude that Theorem 1 is true. We finish the proof.


## 3.2. Nonnegativity

DenoteΘ+=u1,u2,u3∈Θ:u1,u2,u3∈R+.Theorem 2.
Each solution of model (3) that begins with Θ+ remains nonnegative.Proof.
Our object is to prove thatuit≥0i=1,2,3∀t≥t0. Let Ut0=u1t0,u1t0,u1t0∈Θ+ be the initial value of system (3). Assume that ∃ is a constant and t0≤t<ρ such that(17)uit=0,t0≤t<ρ,uiρ=0,uiρ+<0,where i=1,2,3. In view of system (3), one has(18)dϱu1tdtϱ|ui=0=0,dϱu2tdtϱ|ui=0=0,dϱu3tdtϱ|ui=0=0,where i=1,2,3. By Lemma 2, one knows that uiρ+=0i=1,2,3, which contradicts uiρ+<0i=1,2,3 in (17). Thus, one can conclude that uit≥0∀t≥t0. The proof is complete.


## 3.3. Stability and Hopf Bifurcation

In the section, we shall focus on the local stability and the appearance of Hopf bifurcation of system (3). Consider the biological implication of system (3); we only seek the sufficient condition to guarantee the local stability of the positive equilibrium and the emergence of Hopf bifurcations of system (3).Clearly, system (3) has a unique positive equilibrium U0u1∗,u2∗,u3∗ provided that the following condition(19)ℳ1signA=signA1=signA2=signA3,is satisfied, where(20)A=detα11α12−α13α21α22−α23α31α31−α33≠0,A1=detγ1α12−α13γ2α22−α23−γ3α31−α33,A2=detα11γ1−α13α21γ2−α23α31−γ3−α33,A3=detα11α12γ1α21α22γ2α31α31−γ3,u1∗=A1A,u2∗=A2A,u3∗=A3A.The linear system of equation (3) near U0u1∗,u2∗,u3∗ takes(21)dϱu1tdtϱ=a11u1t−σ+a12u2t−σ+a13u3t−σ,dϱu2tdtϱ=a21u1t−σ+a22u2t−σ+a23u3t−σ,dϱu3tdtϱ=a31u1t−σ+a32u2t−σ+a33u3t−σ,where(22)a11=−u1∗α11,a12=−u1∗α12,a13=u1∗α13,a21=−u2∗α21,a22=−u2∗α22,a23=u2∗α23,a31=u3∗α31,a32=u3∗α32,a33=−u3∗α33.The characteristic equation of equation (21) is computed as follows:(23)detsϱ−a11e−sσ−a12e−sσ−a13e−sσ−a21e−sσsϱ−a22e−sσ−a23e−sσ−a31e−sσ−a32e−sσsϱ−a33e−sσ=0,which leads to(24)s3ϱ+b1s2ϱe−sσ+b2sϱe−2sσ+b3e−3sσ=0,where(25)b1=−a11+a22+a33,b2=a11a33+a22a33+a11a22−a13a31−a12a21−a23a32,b3=a13a22a31+a12a21a33+a11a23a32−a11a22a33−a13a21a32.In terms of (24), we have(26)s3ϱe2sσ+b1s2ϱesσ+b2sϱ+b3e−sσ=0.Lets=iφ=φcosπ/2+isinπ/2 be the root of (26); then, one has(27)G1ϱcos2φσ+G2ϱsin2φσ+G3ϱcosφσ+G4ϱsinφσ=0,S1ϱcos2φσ+S2ϱsin2φσ+S3ϱcosφσ+S4ϱsinφσ=0,where(28)G1ϱ=φ3ϱcos3ϱπ2+b2φϱcosϱπ2,G2ϱ=b2φϱsinϱπ2−φ3ϱsin3ϱπ2,G3ϱ=b1φ2ϱcosϱπ+b3,G4ϱ=−b1φ2ϱsinϱπ,S1ϱ=φ3ϱsin3ϱπ2+b2φϱsinϱπ2,S2ϱ=φ3ϱcos3ϱπ2−b2φϱcosϱπ2,S3ϱ=b1φ2ϱsinϱπ,S4ϱ=b1φ2ϱcosϱπ−b3.In (28), let(29)c1=cos3ϱπ2,c2=b2cosϱπ2,c3=b2sinϱπ2,c4=−sin3ϱπ2,c5=b1cosϱπ,c6=b3,c7=−b1sinϱπ,c8=sin3ϱπ2,c9=b2sinϱπ2,c10=cos3ϱπ2,c11=−b2cosϱπ2,c12=b1sinϱπ,c13=b1cosϱπ,c14=−b3,and then, (28) takes the form:(30)G1ϱ=c1φ3ϱ+c2φϱ,G2ϱ=c3φϱ+c4φ3ϱ,G3ϱ=c5φ2ϱ+c6,G4ϱ=c7φ2ϱ,S1ϱ=c8φ3ϱ+c9φϱ,S2ϱ=c10φ3ϱ+c11φϱ,S3ϱ=c12φ2ϱ,S4ϱ=c13φ2ϱ+c14.For (27), we consider two cases.Case 1.
Ifsinφσ=1−cos2φσ, then the first equation of (27) becomes(31)G1ϱ2cos2φσ−1+2G2ϱcosφσ1−cos2φσ+G3ϱcosφσ+G4ϱ1−cos2φσ=0.
Hence,(32)2G2ϱcosφσ1−cos2φσ+G4ϱ1−cos2φσ2=G1ϱ2cos2φσ−1+G3ϱcosφσ2,which leads to(33)ϑ1cos4φσ+ϑ2cos3φσ+ϑ3cos2φσ+ϑ4cosφσ+ϑ5=0,where(34)ϑ1=4G12ϱ+4G22ϱ,ϑ2=4G1ϱG3ϱ+4G2ϱG4ϱ,ϑ3=G32ϱ+G42ϱ−4G12ϱ−4G22ϱ,ϑ4=−2G1ϱG3ϱ+4G2ϱG4ϱ,ϑ5=G12ϱ−G42ϱ.
Setcosφσ=r, and let(35)zr=r4+ϑ2ϑ1r3+ϑ3ϑ1r2+ϑ4ϑ1r+ϑ5ϑ1.
Hence,(36)dzrdr=4r3+3ϑ2ϑ1r2+2ϑ3ϑ1r+ϑ4ϑ1.
Denote(37)r3+3ϑ2ϑ1r2+2ϑ3ϑ1r+ϑ4ϑ1=0.
Setp=r+ϑ2/4ϑ1, and then, (37) takes the form(38)p3+f1p+f2=0,where(39)f1=ϑ32ϑ1−3ϑ2216ϑ12,f2=ϑ2332ϑ13−ϑ2ϑ38ϑ12+ϑ44ϑ1.
Let(40)q1=f222+f133,q2=−1+i32.
It follows from (38) that(41)p1=−f22+q13+−f22−q13,p2=−f22+q1q23+−f22−q1q223,p3=−f22+q1q223+−f22−q1q23.
According to the discussion above, one easily obtains the expression ofcosφσ. Then, one easily gets the expression of sinφσ. Here, we suppose that(42)cosφσ=g1φ,sinφσ=g2φ.
Then,(43)g12φ+g22φ=1.
With the aid of Matlab 7.0, one can get the root (sayφ) of equation (43). Then, one obtains(44)σ1k=1φarccosg1φ+2kπ,k=0,1,2,…Case 2.
Ifsinφσ=−1−cos2φσ, applying the same way to this case, we can obtain(45)cosφσ=h1φ,sinφσ=h2φ.
Then,(46)h12φ+h22φ=1.
By Matlab 7.0, one can obtain the root (sayφ) of equation (46). So,(47)σ2k=1φarccosh1φ+2kπ,k=0,1,2,…
Set(48)σ0=minσ1k,σ2k,k=0,1,2,…,
In the sequel, we check the transversality condition of the appearance of Hopf bifurcation. We give the hypothesis as follows:ℳ2Z11Z21+Z12Z22>0, where(49)Z11=3ϱφ03ϱ−1cos3ϱ−1π2cos2φ0σ0−sin3ϱ−1π2sin2φ0σ0+2φ03ϱσ0cos3ϱπ2cos2φ0σ0−sin3ϱπ2sin2φ0σ0+2b1ϱφ02ϱ−1cos2ϱ−1π2cosφ0σ0−sin2ϱ−1π2sinφ0σ0+b1σ0φ02ϱcosϱπcosφ0σ0−sinϱπsinφ0σ0+b2ϱφ0ϱ−1cosϱ−1π2−b3σ0cosφ0σ0,Z12=3ϱφ03ϱ−1cos3ϱ−1π2sin2φ0σ0+sin3ϱ−1π2cos2φ0σ0+2φ03ϱσ0cos3ϱπ2sin2φ0σ0+sin3ϱπ2cos2φ0σ0+2b1ϱφ02ϱ−1cos2ϱ−1π2sinφ0σ0+sin2ϱ−1π2cosφ0σ0+b1σ0φ02ϱcosϱπsinφ0σ0+sinϱπcosφ0σ0+b2ϱφ0ϱ−1sinϱ−1π2+b3σ0sinφ0σ0,Z21=b3φ0sinφ0σ0−2φ03ϱ+1cos3ϱ+1π2cos2φ0σ0−sin3ϱ+1π2sin2φ0σ0−2b2ϱ02σsinσπ+b1ϱ0σsinσπ2ϱ0cos2ϱ0η0,−b1φ02ϱ+1cos2ϱ+1π2cos2φ0σ0−sin2ϱ+1π2sin2φ0σ0,Z22=b3φ0cosφ0σ0−2φ03ϱ+1cos3ϱ+1π2sin2φ0σ0+sin3ϱ+1π2cos2φ0σ0−b1φ02ϱ+1cos2ϱ+1π2sin2φ0σ0+sin2ϱ+1π2cos2φ0σ0.Lemma 5.
Suppose thatsσ=ξσ+iησ is a root of equation (26) near σ=σ0 such that ξσ0=0,ησ0=φ0; then, one has Reds/dσ|σ=σ0,φ=φ0>0.Proof.
In terms of equation (26), one knows that(50)3ϱs3ϱ−1e2sσdsdσ+2s3ϱe2sσdsdσσ+s+2b1ϱs2ϱ−1esσdsdσ+b1esσs2ϱdsdσσ+s+b2ϱsϱ−1dsdσ−b3e−sσdsdσσ+s=0,which leads to(51)3ϱs3ϱ−1e2sσ+2s3ϱe2sσσ+2b1ϱs2ϱ−1esσ+b1esσs2ϱσ+b2ϱsϱ−1−b3e−sσσdsdσ=b3se−sσ−2s3ϱ+1e2sσ−b1s2ϱ+1esσ.
Then,(52)dsdσ−1=Z1sZ2s−σs,where(53)Z1s=3ϱs3ϱ−1e2sσ+2s3ϱe2sσσ+2b1ϱs2ϱ−1esσ+b1esσs2ϱσ+b2ϱsϱ−1−b3e−sσσ,Z2s=b3se−sσ−2s3ϱ+1e2sσ−b1s2ϱ+1esσ.
So,(54)Redsdσ|σ=σ0,φ=φ0=ReZ1sZ2s|σ=σ0,φ=φ0=Z11Z21+Z12Z22Z212+Z222.
Byℳ2, one obtains(55)Redsdσ−1|σ=σ0,φ=φ0>0.
This finishes the proof of Lemma5.
Suppose that(56)ℳ3b1>0,b1b2>b3,b3>0.Lemma 6.
Provided thatσ=0 and ℳ3 is fulfilled, then system (3) is locally asymptotically stable.Proof.
Ifσ=0, then (24) takes the form(57)λ3+b1λ2+b2λ+b3=0.
Byℳ3, one knows that all roots λi of (57) satisfy argλi>ϱπ/2i=1,2,3,4. Thus, one knows that Lemma 6 is right. The proof finishes.
Based on the investigation above, we have the following conclusion.Theorem 3.
Provided thatℳ1–ℳ3 hold true, then U0u1∗,u1∗,u3∗ of system (3) is locally asymptotically stable if σ∈0,σ0, and a Hopf bifurcation exists near U0u1∗,u1∗,u3∗ if σ=σ0.


## 4. Computer Simulations

In this paper, we apply implicit Euler’s scheme which is introduced in [52] to carry our numerical simulations. Given the following fractional-order predator-prey model(58)dϱu1tdtϱ=u1t8−3u1t−σ−2u2t−σ+0.3u3t−σ,dϱu2tdtϱ=u2t7−2u1t−σ−2u2t−σ+0.3u3t−σ,dϱu3tdtϱ=u3t6+0.2u1t−σ+0.2u2t−σ−3u3t−σ,we can easily know that model (58) possesses the equilibrium point U01.0000,2.8384,2.2559. Let ϱ=0.79. Then, r1=8,γ2=7,γ3=6,α11=3,α12=2,α13=0.3,α21=2,α22=2,α23=0.3,α31=0.2,α32=0.2,andα33=3. Hence, a11=−3,a12=−2,a13=0.3,a21=−5.6788,a22=−5.6788,, a23=0.6768,a31=0.4512,a32=0.4512,a33=−6.7677,φ0=0.3127,andσ0=0.192. Then, b1=15.4465,b2=113.7636,b3=151.8836,Z11=0.3742,Z12=0.6425,Z21=0.5004,andZ22=0.7157. Thus, the hypothesis ℳ1–ℳ3 of Theorem 3 hold true. If σ∈0,0.192, the positive equilibrium point U01.0000,2.8384,2.2559 of system (58) is locally asymptotically stable. For this case, we choose σ=0.18<σ0=0.192. The simulation plots are displayed in Figures 1–3. Figures 1–3 reveal that if σ is less than σ0=0.192, then the densities of competing species u1 and u2 of system (58) will be tardily close to 1.0000,2.8384, respectively, and the density of cooperating species u3 of system (58) will be tardily close to 2.2559. When σ∈0.192,+∞, then system (58) loses its stability and Hopf bifurcation behavior emerges. For this case, we choose σ=0.33. The numerical simulation diagrams are presented in Figures 4–6. Figures 4–6 confirm that when σ is greater than σ0=0.192, then the densities of competing species u1 and u2 and the density of cooperating species u3 will remain periodic motion around the positive equilibrium point U01.0000,2.8384,2.2559, namely, Hopf bifurcation phenomenon takes place around U01.0000,2.8384,2.2559. In order to illustrate the bifurcation phenomenon intuitively, we plot the bifurcation diagrams (see Figures 7–9). From Figures 7–9, we can easily see that the bifurcation value is 0.192. In addition, we give the relationship table between φ0 and σ0 in Table 1.Figure 1
State variableu1t versus time t of system (58) when σ=0.18<σ0=0.192.Figure 2
. State variableu2t versus time t of system (58) when σ=0.18<σ0=0.192.Figure 3
State variableu3t versus time t of system (58) when σ=0.18<σ0=0.192.Figure 4
State variableu1t versus time t when σ=0.33>σ0=0.192.Figure 5
State variableu2t versus time t when σ=0.33>σ0=0.192.Figure 6
State variableu3t versus time t when σ=0.33>σ0=0.192.Figure 7
Bifurcation diagram for system (58): σ versus u1.Figure 8
Bifurcation diagram for system (58): σ versus u2.Figure 9
Bifurcation diagram for system (58): σ versus u3.Table 1
The magnitude relation ofφ0 and σ0 for system (58).
φ0σ00.31270.19200.34990.21730.49570.32180.52000.34000.68800.47320.78320.55431.01080.76511.04830.80221.15560.9123


## 5. Conclusions

Based on the previous works on predator-prey models and noticing the effect of feedback delay between the predators and preys, we have established a new fractional-order competitor-competitor-mutualist Lotka–Volterra model. We have found the conditions to ensure the existence, uniqueness, and nonnegativity of the involved fractional-order prey-predator model. Applying the Laplace transform, stability theorem, and Hopf bifurcation theory of fractional-order dynamical systems, we set up a novel sufficient criterion to guarantee the local stability and the emergence of Hopf bifurcation of the involved fractional-order predator-prey system. The study shows that the time delay has a vital effect on controlling the stability. To check the rationality of theoretical predictions, we design the program to implement simulation experiments. The established analytical conclusions have crucial guiding significance in preserving the coexistence of biological populations. In this paper, we assume that the feedback time delay of different species is same. Of course, we can also deal with the different delay cases. We will focus on this topic in near future.



---

*Source: 6402459-2021-02-22.xml*