# Stability and Bifurcation Analysis of a Fractional-Order Food Chain Model with Two Time Delays

Authors: Hao Qi; Wencai Zhao
Journal: Journal of Mathematics
 (2022)

---

## Abstract

In this study, the stability and bifurcation problems of a fractional food chain system with two kinds of delays are studied. Firstly, the nonnegative, bounded, and unique properties of the solutions of the system are proved. The asymptotic stability of the equilibrium points of the system is discussed. Furthermore, the global asymptotic stability of the positive equilibrium point is deduced by using Lyapunov function method. Secondly, the system takes two kinds of time delays as bifurcation parameters and calculates the critical values of Hopf bifurcation accurately. The results show that Hopf bifurcation can advance with increasing fractional order and another delay. In conclusion, numerical simulation verifies and illustrates the theoretical results.


---

## Body

## 1. Introduction

In the ecosystem, no species exists in isolation. The different populations are all related to each other. Predator relationship, competition relationship, reciprocity relationship, and parasitic relationship are the main population relationships. In these major relationships, predator-prey relationship is universal in nature and is of great significance to complex ecosystems. It is precisely because of the important application background and practical value of the predation system that the food chain system has been researched extensively by many scholars [1–5].In nature, the phenomenon of time delay is exhibited universally in biological population. The phenomenon of time delay is mainly caused by many factors such as gestation, maturation, and food digestion of population. The phenomenon of time delay signifies that the related properties of the system are related to not only the present state but also the previous period. Aiello and Freedman studied a single population system with a time delay and stage structure [6]. Beddington et al. [7] proved that time delay could affect the stability of the dynamical model. Gazi et al. [8] researched the influence of harvest and discrete time delay on the prey-predator populations and obtained the discrete time-delay length required to remain the stability of the system. Jana et al. [9] analyzed the time-delay predator-prey system including prey shelter and demonstrated the global asymptotic stability of the system. Yan et al. [10] considered the predator-prey model with delayed reaction diffusion and analyzed the global asymptotic stability of the positive equilibrium point of the model. Vinoth et al. [11] put forward a delayed prey-predator system with additive Allee effect, and the local asymptotic stability of the model at equilibrium point was studied. Numerous studies have shown that a population system with time delay could exhibit more complex nonlinear dynamic behaviors. Therefore, time delay has a profound impact on the stability behavior of biological systems.Differential equation theory has been widely used in automation system, aerospace technology, information engineering, and so on. In these practical applications, the system usually has some parameters. If the parameters of the system change, the topological structure of the phase diagram in phase space also changes; then, the phenomenon is called bifurcation [12–14]. Hopf bifurcation theory has become a classical tool to research the generation and extinction periodic solutions of small amplitude differential equations. When a parameter passes a marginal value, the equilibrium point will lose stability and a periodic solution will appear [15–19]. Deka et al. [15] proposed and analyzed a one-predator and two-prey system with a general Gauss type, and the stability and direction of the Hopf bifurcation were proved by regarding the mortality of the predator as the bifurcation argument. In [16], a predator-prey model with discrete time delay of habitat complexity and sanctuary for prey was proposed and the occurrence criterion of Hopf bifurcation was obtained by taking the time lag as argument. In [20], Guo et al. established a food chain system with a couple of time lags and Holling II type functions:(1)x′t=xtr1−a1xt−yt−τ11+mxt,y′t=yta2xt−τ21+mxt−r2−a3zt−τ1,z′t=zta4yt−τ2−r3.Among them, the biological significance of each parameter of system (1) is well illustrated in Table 1.Table 1
Biological significance of symbols.
SymbolsBiological significancextThe density of prey population at timetytThe density of the primary predator population at timetztThe density of the top predator population at timetr1The intrinsic growth rate of prey populationr2The death rate of the primary predator populationr3The death rate of the top predator populationa1The internal competition rate of the prey populationa2The nutrient conversion rate from prey to primary predatora3The rate of capture by top predators on primary predatorsa4The digestibility of top predators to primary predatorsmThe semisaturation of predatorτ1The capture timeτ2The maturity timeAt the same time, the existence of the positive equilibrium point was proved, and the occurrence criterion of Hopf bifurcation was obtained by taking the time lag as the parameter.Fractional-order calculus is a method that rises recently. It is a method that extends the ordinary integral calculus to nonintegral calculus [21–27]. So far, fractional calculus has been applied to many domains, such as neural network [28], medicine [29], finance system [30], and safety communication [31]. A great deal of studies have proved that the fractional dynamical system is to a higher degree suitable to biological systems because the fractional differential is connected with the entire time zone, while the integer differential is only related to a particular moment. Because biological systems generally have the characteristics of heredity and memory, so more and more scholars believe that the method of fractional calculus can better characterize the behavior of biological system. At present, some scholars have spread the classical integer-order differential systems to the fractional-order differential systems [32–37]. Rihan et al. [32] studied a fractional-order food chain model with time delay as well as infection in predators; sufficient criterion for asymptotic stability of the stable condition of the model was established. Huang et al. [36] discovered that the bifurcation dynamics of the model could be resultfully controlled as long as other parameters of the system are determined, and the extended feedback delay or fractional order is carefully adjusted.Based on the above discussion, model (1) is extended in this study to obtain the following fractional-order food chain model:(2)Dθxt=xtr1−a1xt−yt−τ11+mxt,Dθyt=yta2xt−τ21+mxt−r2−a3zt−τ1,Dθzt=zta4yt−τ2−r3,where 0<θ≤1; the biological significance of each variable and parameter of model (2) is the same as that of model (1). The initial conditions are xt=ζ1t≥0,yt=ζ2t≥0,zt=ζ3t≥0,andt∈−maxτ1,τ2,0. The model is established on the sense of Caputo derivative.The rest of the study is organized as follows. Several definitions as well as lemmas are addressed in Section2. In Section 3, the corresponding nondelay system of (2) is discussed. The Hopf bifurcation of system (2) is studied in Section 4. Some numerical simulations are presented in Section 5. Conclusions are drawn in the end.


## 2. Preliminaries

For the theoretical derivation, we first give the relevant definitions and lemmas of Caputo calculus.Definition 1 (see [21]).
The fractional integral of orderθ for a function ft is defined as(3)Iθft=1Γθ∫t0tt−sθ−1fsds,where t≥t0, θ>0, and Γ⋅ is the Gamma function, Γθ=∫0∞tθ−1e−tdt.Definition 2 (see [21]).
The Caputo fractional derivative of orderθ for a function ft is defined as(4)Dθft=1Γn−θ∫t0tfnst−sθ−n+1ds,where n is a positive integer, n−1<θ≤n, and t≥t0. So, specifically, if 0<θ<1,(5)Dθft=1Γ1−θ∫t0tf′st−sθds.Lemma 1 (see [22]).
Definew∈Cθt0,T,R. Suppose that there exist t1∈t0,T, such that wt1=0 and wt>0 for t0≤t<t1; then,(6)Dθwt1<0.Lemma 2 (see [21]).
Defineθ>0, n−1<θ≤n. Suppose mt is n times continuous differentiable function and Dθmt is piecewise continuous on t0,∞; we have(7)LDθmt=sθϒs−∑k=0n−1sθ−k−1mkt0,where ϒs=Lmt.Lemma 3 (see [38]).
AssumeM represents the complex plane, for ∀a>0 and b>0, and Q∈M; then,(8)Ltb−1Ea,bQta=sa−bsa−Q,for ℝs>Q1/a, ℝs signifies the real part of the complex number s, and Ea,b is the following Mittag–Leffler function described by(9)Ea,bz=∑n=0∞znΓan+b.Lemma 4 (see [24]).
Consider the system:(10)DθXt=ϱt,X,t0>0,with initial condition Xt0=Xt0, where θ∈0,1 and ϱ:t0,∞×η⟶Rn, η⊆Rn, if ϱt,X meets the local Lipschitz criteria with respect to X∈Rn:(11)ϱt,X−ϱt,X¯≤δX−X¯,so it has a unique solution of (10) on t0,∞, where(12)Xu1,u2,…,un−X¯v1,v2,…,vn=∑i=1nui−vi,ui,vi∈R.Lemma 5 (see [23]).
Define0<θ<1, Xt∈Rn, and f:Rn⟶Rn; think about the following nonlinear fractional system of the same order:(13)DθXt=fXt,X0=X0.If the eigenvaluesλii=1,…,n of the Jacobian matrix corresponding to the equilibrium point of the system meet the following criterion argλi>θπ/2,i=1,2,…n, then system (13) is asymptotically stable.Lemma 6 (see [26]).
Assumelt∈R+ is a continuous and differentiable function. For ∀t≥t0, Dθlt−l∗−l∗lnlt/l∗≤1−l∗/ltDθlt, l∗∈R+, ∀θ∈0,1.Lemma 7 (see [23]).
Think about the undern-dimensional linear fractional-order time-delay system:(14)Dθ1x1t=b11x1t−τ11+b12x2t−τ12+⋯+b1nxnt−τ1n,Dθ2x2t=b21x1t−τ21+b22x2t−τ22+⋯+b2nxnt−τ2n,⋮Dθnxnt=bn1x1t−τn1+bn2x2t−τn2+⋯+bnnxnt−τnn.Among them,∀θi∈0,1, and the initial conditions xit=ξit are provided for −max1≤i,j≤nτij=−τmax≤t≤0, i=1,2,…,n. It is defined as(15)Λs=sθ1−b11e−sτ11−b12e−sτ12⋯−b1ne−sτ1n−b21e−sτ21sθ2−b22e−sτ22⋯−b2ne−sτ2n⋮⋮⋱⋮−bn1e−sτn1−bn2e−sτn2⋯sθn−bnne−sτnn.If all roots ofdetΛs=0 have negative real parts, so the zero solution of system (14) is Lyapunov globally asymptotically stable.


## 3. Analysis of the Nondelayed Model

First, we research the delay-free system of (2):(16)Dθxt=xtr1−a1xt−yt1+mxt,Dθyt=yta2xt1+mxt−r2−a3zt,Dθzt=zta4yt−r3.The nonnegativity and boundedness, existence, and uniqueness of solutions about systems (2) and (16) are discussed in Sections 3.1 and 3.2. The local stability of the equilibrium points of system (2) is discussed, and the global asymptotic stability of the positive equilibrium point of system (2) is demonstrated in Section 3.3.


### 3.1. Nonnegativity and Boundedness of Solutions

Think about the biological implications of reality, it is significant to analyze the nonnegativity of the system. To prove the following theorem, letR+ denote the collection of entire positive real numbers containing 0, η+=x,y,z∈η:x,y,z∈R+.Theorem 1.
The solutions about system (16) from η+ are nonnegative and uniformly bounded.Proof.
Whent=0, then x0>0; we desire to obtain the solution xt from η+ is nonnegative, i.e., xt≥0, for ∀t≥0. Suppose it exhibits a constant t′>0, xt′<0; according to xt which is a continuous function, there exists t″∈0,t′ and xt″=0. Define t1=mint″∈0,t′|xt″=0; then, when t=t1>0, from system (16), one obtains Dθxt|t=t1=xtr1−a1xt−yt/1+mxt=0. However, according to the definition of t1, x0>0 and xt1=0; moreover, xt>0,t∈0,t1; by Lemma 1, we have Dθxt1<0. Hence, we derive a contradiction; therefore, xt≥0,∀t∈0,∞. Likewise, we can demonstrate yt,zt≥0,∀t∈0,∞.
For boundedness, we think about the following function:(17)Wxt,yt,zt=xt+1a2yt+a3a2a4zt.
According to system (16), one has(18)DθWt=r1xt−a1x2t−1a2r2yt−a3r3a2a4zt,DθWt+k^Wt=−a1xt−r1+k^2a12+r1+k^24a1−r2−k^a2yt−a3a2a4r3−k^zt≤−a1xt−r1+k^2a12+r1+k^24a1≤r1+k^24a1=A,where k^=minr2,r3>0. Therefore,(19)DθWt+k^Wt≤A.
By Lemma2, making Laplace transform of both sides of (19), we obtain(20)sθϒs−sθ−1W0+k^ϒs≤As,where ϒs=LWt. From this, we can obtain(21)ϒs≤sθ−1W0sθ+k^+Assθ+k^.
Making inverse Laplace transform of (21), then(22)Wt≤W0L−1sθ−1sθ+k^+AL−1sθ−θ+1sθ+k^.
By Lemma3, one has(23)Wt≤W0Eθ,1−k^tθ+AtθEθ,θ+1−k^tθ.
According to(24)Eϑ,ιz=zEϑ,ϑ+ιz+1Γι,so we have(25)Eθ,1−k^tθ=−k^tθEθ,θ+1−k^tθ+1Γ1,tθEθ,θ+1−k^tθ=−1k^Eθ,1−k^tθ−1.
Hence,(26)Wt≤W0−Ak^Eθ,1−k^tθ+Ak^,where, if t⟶∞, we have Eθ,1−k^tθ⟶0.
Furthermore, the setD attracts all the solutions of system (16), where(27)D=x,y,z∈η+|xt+1a2yt+a3a2a4zt≤1k^r1+k^24a1+ϵ,ϵ>0.


### 3.2. Existence and Uniqueness of Solutions

Theorem 2.
System (16) only exhibits a solution Xt=xt,yt,zt∈η+ for any given initial value Xt0=xt0,yt0,zt0∈η+.Proof.
According to Theorem1, the solutions of system (16) from η+ are nonnegative and uniformly bounded; then, there exists a constant P, such that maxxt,yt,zt≤P. Define a mapping QX=Q1X,Q2X,Q3X, in which(28)Q1X=xtr1−a1xt−yt1+mxt,Q2X=yta2xt1+mxt−r2−a3zt,Q3X=zta4yt−r3.
LetXandX¯ be any two solutions to system (16); we can derive(29)QX−QX¯=Q1X−Q1X¯+Q2X−Q2X¯+Q3X−Q3X¯=xr1−a1x−y1+mx−x¯r1−a1x¯−y¯1+mx¯+ya2x1+mx−r2−a3z−y¯a2x¯1+mx¯−r2−a3z¯+za4y−r3−z¯a4y¯−r3≤r1x−x¯+a1x2−x¯2+xy−x¯y¯+mxyx¯−mx¯y¯x1+mx1+mx¯+r2y−y¯+a3yz−y¯z¯+a2xy−x¯y¯+mxyx¯−mx¯y¯x1+mx1+mx¯+a4yz−y¯z¯+r3z−z¯≤r1x−x¯+2a1Px−x¯+Py−y¯+Px−x¯+mP2y−y¯+mP2x−x¯+mP2x−x¯+r2y−y¯+a3Pz−z¯+a3Py−y¯+a2Py−y¯+Px−x¯+mP2y−y¯+mP2x−x¯+mP2x−x¯+a4Pz−z¯+a4Py−y¯+r3z−z¯=r1+2a1P+P+2mP2+a2P+2a2mP2x−x¯+P+mP2+r2+a3P+a2P+a2mP2+a4Py−y¯+a3P+a4P+r3z−z¯≤δX−X¯,where δ=maxM1,M2,M3 and M1=r1+2a1P+P+2mP2+a2P+2a2mP2, M2=P+mP2+r2+a3P+a2P+a2mP2+a4P, and M3=a3P+a4P+r3. Hence, QX meets the Lipschitz criteria about X. By Lemma 4, it has only a solution Xt of system (16) with initial value Xt0=xt0,yt0,zt0.


### 3.3. Stability Analysis of Balance Point

For analyzing the possible equilibrium of system (16), we first present the following assumptions:(i)
ℍ1: r1a2−mr2−a1r2>0.(ii)
ℍ2: a4r1−r3>0 and a2x∗−r21+mx∗>0, in which x∗ is the positive root of the following equation:(30)−ma1a4x2+mr1a4−a1a4x+a4r1−r3=0.We can find the following four biologically feasible equilibrium points:(1)
E0=0,0,0(2)
E1=r1/a1,0,0(3)
E2=xˇ,yˇ,0; it exhibits if the condition ℍ1 is true, where xˇ=r2/a2−mr2 and yˇ=a2r1a2−mr2−a1r2/a2−mr22(4)
E3=x∗,y∗,z∗; it exhibits if the condition ℍ2 is true, where y∗=r3/a4 and z∗=a2x∗−r21+mx∗/1+mx∗a3The Jacobian matrix about system (16) at arbitrary point x,y,z is as follows(31)J=r1−2a1x−y1+mx2−x1+mx0a2y1+mx2a2x1+mx−r2−a3z−a3y0a4za4y−r3.Theorem 3.
The trivial equilibriumE0 of system (16) is unstable.Proof.
The Jacobian matrix atE0=0,0,0 is as follows:(32)JE0=r1000−r2000−r3.
Letλ=sθ; the characteristic equation of (32) is(33)λ−r1λ+r2λ+r3=0.
So, the roots of the characteristic equation areλ1=r1,λ2=−r2,andλ3=−r3. Therefore,(34)argλ1=0<θπ2,argλ2=argλ3=π>θπ2.
Consequently, by Lemma5, the trivial equilibrium E0 of system (16) is unstable.Theorem 4.
Ifa2r1/a1+mr1−r2<0 is met, the boundary equilibrium E1 of system (16) is locally asymptotically stable.Proof.
The corresponding Jacobian matrix atE1=r1/a1,0,0 is shown below:(35)JE1=−r1−r1a1+mr100a2r1a1+mr1−r2000−r3.
At this point, the characteristic equation corresponding to (35) is(36)λ+r1λ+r3λ−a2r1a1+mr1−r2=0.
So, the roots of the characteristic equation areλ1=−r1,λ2=a2r1/a1+mr1−r2,λ3=−r3. Owing to the assumptions,(37)argλ1=argλ2=argλ3=π>θπ2.
Consequently,E1 of system (16) is locally asymptotically stable by Lemma 5.Theorem 5.
In the case ofℍ1, if one of the following conditions is met,(1)
r1a2−mr2−a1r2=a1a2/m<r3a2−mr22/a2a4(2)
r1a2−mr2−a1r2<mina1a2/m,r3a2−mr22/a2a4(3)
a1a2/m<r1a2−mr2−a1r2<r3a2−mr22/a2a4 and p1/2p2<cosθπ/2, where p1 and p2 are given in the following proof; then, the boundary equilibrium E2=xˇ,yˇ,0 of system (16) is locally asymptotically stableProof.
In the case ofℍ1, the boundary equilibrium E2=xˇ,yˇ,0 of system exhibits. The Jacobian matrix at E2=xˇ,yˇ,0 is as follows:(38)JE2=mxˇyˇ1+mxˇ2−a1xˇ−xˇ1+mxˇ0a2yˇ1+mxˇ20−a3yˇ00a4yˇ−r3.
The characteristic equation of (38) is(39)λ2−p1λ+p2λ−p3=0,where p1=mxˇyˇ/1+mxˇ2−a1xˇ, p2=a2xˇyˇ/1+mxˇ3, and p3=a4yˇ−r3.(1)
Ifr1a2−mr2−a1r2=a1a2/m<r3a2−mr22/a2a4, we can find that p3<0 and p1=0; all characteristic roots of equation (39) are λ1=p2i,λ2=−p2i,andλ3=p3; therefore,(40)argλ1=argλ2=π2>θπ2,argλ3=π>θπ2.Hence, by Lemma5, the equilibrium E2 about system (16) is locally asymptotically stable.(2)
Ifr1a2−mr2−a1r2<mina1a2/m,r3a2−mr22/a2a4, we can find that p3<0 and p1<0. Obviously,(41)argλ1>θπ2,argλ2>θπ2,argλ3=π>θπ2.Accordingly, by Lemma5, the equilibrium E2 about system (16) is locally asymptotically stable.(3)
Using the given conditions, we can obtain all characteristic roots of equation (39) are λ1=p1/2+4p2−p12/2i,λ2=p1/2−4p2−p12/2i,andλ3=p3. Owing to p1/2p2<cosθπ/2, therefore,(42)argλ1>θπ2,argλ2>θπ2,argλ3=π>θπ2.
As a result, by Lemma5, the equilibrium E2 about system (16) is locally asymptotically stable.Theorem 6.
Ifℍ3:a1>my∗/1+mx∗ is satisfied, then the positive equilibrium E3=x∗,y∗,z∗ about system (16) is globally asymptotically stable.Proof.
The Jacobian matrix atE3=x∗,y∗,z∗ is as follows:(43)JE3=mx∗y∗1+mx∗2−a1x∗−x∗1+mx∗0a2y∗1+mx∗20−a3y∗0a4z∗0.
The characteristic equation of (43) is(44)λ3+A1λ2+A2λ+A3=0,where A1=−mx∗y∗/1+mx∗2−a1x∗, A2=a2x∗y∗/1+mx∗3+a3a4y∗z∗, and A3=−a3a4y∗z∗mx∗y∗/1+mx∗2−a1x∗. Based on our assumptions, we have(45)A1>0,A1A31A2>0,A3>0.
By the Routh–Hurwitz criterion, all roots of (44) are negative real parts; therefore,(46)argλ1>θπ2,argλ2>θπ2,argλ3>θπ2.
As a result, by Lemma5, the equilibrium E3 about system (16) is locally asymptotically stable.
Let us consider the Lyapunov function:(47)Vx,y,z=x−x∗−x∗lnxx∗+1+mx∗a2y−y∗−y∗lnyy∗+a31+mx∗a2a4z−z∗−z∗lnzz∗.
Obviously,Vx,y,z>0 for any x,y,z>0, except for the positive equilibrium E3=x∗,y∗,z∗.
By Lemma6, we have(48)DθVx,y,z≤1−x∗xDθxt+1+mx∗a21−y∗yDθyt+a31+mx∗a2a41−z∗zDθzt=x−x∗r1−a1x−y1+mx+1+mx∗a2y−y∗a2x1+mx−r2−a3z+a31+mx∗a2a4z−z∗a4y−r3=x−x∗a1x∗−a1x−y1+mx+y∗1+mx+1+mx∗a2y−y∗a2x1+mx−a2x∗1+mx∗−a3z+a3z∗+a31+mx∗a2a4y−y∗a4y−y∗=−a1x−x∗2+my∗x−x∗21+mx1+mx∗.
Sincea1>my∗/1+mx∗, then we have DθVx,y,z≤0. Thus, E3 is globally asymptotically stable.


## 3.1. Nonnegativity and Boundedness of Solutions

Think about the biological implications of reality, it is significant to analyze the nonnegativity of the system. To prove the following theorem, letR+ denote the collection of entire positive real numbers containing 0, η+=x,y,z∈η:x,y,z∈R+.Theorem 1.
The solutions about system (16) from η+ are nonnegative and uniformly bounded.Proof.
Whent=0, then x0>0; we desire to obtain the solution xt from η+ is nonnegative, i.e., xt≥0, for ∀t≥0. Suppose it exhibits a constant t′>0, xt′<0; according to xt which is a continuous function, there exists t″∈0,t′ and xt″=0. Define t1=mint″∈0,t′|xt″=0; then, when t=t1>0, from system (16), one obtains Dθxt|t=t1=xtr1−a1xt−yt/1+mxt=0. However, according to the definition of t1, x0>0 and xt1=0; moreover, xt>0,t∈0,t1; by Lemma 1, we have Dθxt1<0. Hence, we derive a contradiction; therefore, xt≥0,∀t∈0,∞. Likewise, we can demonstrate yt,zt≥0,∀t∈0,∞.
For boundedness, we think about the following function:(17)Wxt,yt,zt=xt+1a2yt+a3a2a4zt.
According to system (16), one has(18)DθWt=r1xt−a1x2t−1a2r2yt−a3r3a2a4zt,DθWt+k^Wt=−a1xt−r1+k^2a12+r1+k^24a1−r2−k^a2yt−a3a2a4r3−k^zt≤−a1xt−r1+k^2a12+r1+k^24a1≤r1+k^24a1=A,where k^=minr2,r3>0. Therefore,(19)DθWt+k^Wt≤A.
By Lemma2, making Laplace transform of both sides of (19), we obtain(20)sθϒs−sθ−1W0+k^ϒs≤As,where ϒs=LWt. From this, we can obtain(21)ϒs≤sθ−1W0sθ+k^+Assθ+k^.
Making inverse Laplace transform of (21), then(22)Wt≤W0L−1sθ−1sθ+k^+AL−1sθ−θ+1sθ+k^.
By Lemma3, one has(23)Wt≤W0Eθ,1−k^tθ+AtθEθ,θ+1−k^tθ.
According to(24)Eϑ,ιz=zEϑ,ϑ+ιz+1Γι,so we have(25)Eθ,1−k^tθ=−k^tθEθ,θ+1−k^tθ+1Γ1,tθEθ,θ+1−k^tθ=−1k^Eθ,1−k^tθ−1.
Hence,(26)Wt≤W0−Ak^Eθ,1−k^tθ+Ak^,where, if t⟶∞, we have Eθ,1−k^tθ⟶0.
Furthermore, the setD attracts all the solutions of system (16), where(27)D=x,y,z∈η+|xt+1a2yt+a3a2a4zt≤1k^r1+k^24a1+ϵ,ϵ>0.


## 3.2. Existence and Uniqueness of Solutions

Theorem 2.
System (16) only exhibits a solution Xt=xt,yt,zt∈η+ for any given initial value Xt0=xt0,yt0,zt0∈η+.Proof.
According to Theorem1, the solutions of system (16) from η+ are nonnegative and uniformly bounded; then, there exists a constant P, such that maxxt,yt,zt≤P. Define a mapping QX=Q1X,Q2X,Q3X, in which(28)Q1X=xtr1−a1xt−yt1+mxt,Q2X=yta2xt1+mxt−r2−a3zt,Q3X=zta4yt−r3.
LetXandX¯ be any two solutions to system (16); we can derive(29)QX−QX¯=Q1X−Q1X¯+Q2X−Q2X¯+Q3X−Q3X¯=xr1−a1x−y1+mx−x¯r1−a1x¯−y¯1+mx¯+ya2x1+mx−r2−a3z−y¯a2x¯1+mx¯−r2−a3z¯+za4y−r3−z¯a4y¯−r3≤r1x−x¯+a1x2−x¯2+xy−x¯y¯+mxyx¯−mx¯y¯x1+mx1+mx¯+r2y−y¯+a3yz−y¯z¯+a2xy−x¯y¯+mxyx¯−mx¯y¯x1+mx1+mx¯+a4yz−y¯z¯+r3z−z¯≤r1x−x¯+2a1Px−x¯+Py−y¯+Px−x¯+mP2y−y¯+mP2x−x¯+mP2x−x¯+r2y−y¯+a3Pz−z¯+a3Py−y¯+a2Py−y¯+Px−x¯+mP2y−y¯+mP2x−x¯+mP2x−x¯+a4Pz−z¯+a4Py−y¯+r3z−z¯=r1+2a1P+P+2mP2+a2P+2a2mP2x−x¯+P+mP2+r2+a3P+a2P+a2mP2+a4Py−y¯+a3P+a4P+r3z−z¯≤δX−X¯,where δ=maxM1,M2,M3 and M1=r1+2a1P+P+2mP2+a2P+2a2mP2, M2=P+mP2+r2+a3P+a2P+a2mP2+a4P, and M3=a3P+a4P+r3. Hence, QX meets the Lipschitz criteria about X. By Lemma 4, it has only a solution Xt of system (16) with initial value Xt0=xt0,yt0,zt0.


## 3.3. Stability Analysis of Balance Point

For analyzing the possible equilibrium of system (16), we first present the following assumptions:(i)
ℍ1: r1a2−mr2−a1r2>0.(ii)
ℍ2: a4r1−r3>0 and a2x∗−r21+mx∗>0, in which x∗ is the positive root of the following equation:(30)−ma1a4x2+mr1a4−a1a4x+a4r1−r3=0.We can find the following four biologically feasible equilibrium points:(1)
E0=0,0,0(2)
E1=r1/a1,0,0(3)
E2=xˇ,yˇ,0; it exhibits if the condition ℍ1 is true, where xˇ=r2/a2−mr2 and yˇ=a2r1a2−mr2−a1r2/a2−mr22(4)
E3=x∗,y∗,z∗; it exhibits if the condition ℍ2 is true, where y∗=r3/a4 and z∗=a2x∗−r21+mx∗/1+mx∗a3The Jacobian matrix about system (16) at arbitrary point x,y,z is as follows(31)J=r1−2a1x−y1+mx2−x1+mx0a2y1+mx2a2x1+mx−r2−a3z−a3y0a4za4y−r3.Theorem 3.
The trivial equilibriumE0 of system (16) is unstable.Proof.
The Jacobian matrix atE0=0,0,0 is as follows:(32)JE0=r1000−r2000−r3.
Letλ=sθ; the characteristic equation of (32) is(33)λ−r1λ+r2λ+r3=0.
So, the roots of the characteristic equation areλ1=r1,λ2=−r2,andλ3=−r3. Therefore,(34)argλ1=0<θπ2,argλ2=argλ3=π>θπ2.
Consequently, by Lemma5, the trivial equilibrium E0 of system (16) is unstable.Theorem 4.
Ifa2r1/a1+mr1−r2<0 is met, the boundary equilibrium E1 of system (16) is locally asymptotically stable.Proof.
The corresponding Jacobian matrix atE1=r1/a1,0,0 is shown below:(35)JE1=−r1−r1a1+mr100a2r1a1+mr1−r2000−r3.
At this point, the characteristic equation corresponding to (35) is(36)λ+r1λ+r3λ−a2r1a1+mr1−r2=0.
So, the roots of the characteristic equation areλ1=−r1,λ2=a2r1/a1+mr1−r2,λ3=−r3. Owing to the assumptions,(37)argλ1=argλ2=argλ3=π>θπ2.
Consequently,E1 of system (16) is locally asymptotically stable by Lemma 5.Theorem 5.
In the case ofℍ1, if one of the following conditions is met,(1)
r1a2−mr2−a1r2=a1a2/m<r3a2−mr22/a2a4(2)
r1a2−mr2−a1r2<mina1a2/m,r3a2−mr22/a2a4(3)
a1a2/m<r1a2−mr2−a1r2<r3a2−mr22/a2a4 and p1/2p2<cosθπ/2, where p1 and p2 are given in the following proof; then, the boundary equilibrium E2=xˇ,yˇ,0 of system (16) is locally asymptotically stableProof.
In the case ofℍ1, the boundary equilibrium E2=xˇ,yˇ,0 of system exhibits. The Jacobian matrix at E2=xˇ,yˇ,0 is as follows:(38)JE2=mxˇyˇ1+mxˇ2−a1xˇ−xˇ1+mxˇ0a2yˇ1+mxˇ20−a3yˇ00a4yˇ−r3.
The characteristic equation of (38) is(39)λ2−p1λ+p2λ−p3=0,where p1=mxˇyˇ/1+mxˇ2−a1xˇ, p2=a2xˇyˇ/1+mxˇ3, and p3=a4yˇ−r3.(1)
Ifr1a2−mr2−a1r2=a1a2/m<r3a2−mr22/a2a4, we can find that p3<0 and p1=0; all characteristic roots of equation (39) are λ1=p2i,λ2=−p2i,andλ3=p3; therefore,(40)argλ1=argλ2=π2>θπ2,argλ3=π>θπ2.Hence, by Lemma5, the equilibrium E2 about system (16) is locally asymptotically stable.(2)
Ifr1a2−mr2−a1r2<mina1a2/m,r3a2−mr22/a2a4, we can find that p3<0 and p1<0. Obviously,(41)argλ1>θπ2,argλ2>θπ2,argλ3=π>θπ2.Accordingly, by Lemma5, the equilibrium E2 about system (16) is locally asymptotically stable.(3)
Using the given conditions, we can obtain all characteristic roots of equation (39) are λ1=p1/2+4p2−p12/2i,λ2=p1/2−4p2−p12/2i,andλ3=p3. Owing to p1/2p2<cosθπ/2, therefore,(42)argλ1>θπ2,argλ2>θπ2,argλ3=π>θπ2.
As a result, by Lemma5, the equilibrium E2 about system (16) is locally asymptotically stable.Theorem 6.
Ifℍ3:a1>my∗/1+mx∗ is satisfied, then the positive equilibrium E3=x∗,y∗,z∗ about system (16) is globally asymptotically stable.Proof.
The Jacobian matrix atE3=x∗,y∗,z∗ is as follows:(43)JE3=mx∗y∗1+mx∗2−a1x∗−x∗1+mx∗0a2y∗1+mx∗20−a3y∗0a4z∗0.
The characteristic equation of (43) is(44)λ3+A1λ2+A2λ+A3=0,where A1=−mx∗y∗/1+mx∗2−a1x∗, A2=a2x∗y∗/1+mx∗3+a3a4y∗z∗, and A3=−a3a4y∗z∗mx∗y∗/1+mx∗2−a1x∗. Based on our assumptions, we have(45)A1>0,A1A31A2>0,A3>0.
By the Routh–Hurwitz criterion, all roots of (44) are negative real parts; therefore,(46)argλ1>θπ2,argλ2>θπ2,argλ3>θπ2.
As a result, by Lemma5, the equilibrium E3 about system (16) is locally asymptotically stable.
Let us consider the Lyapunov function:(47)Vx,y,z=x−x∗−x∗lnxx∗+1+mx∗a2y−y∗−y∗lnyy∗+a31+mx∗a2a4z−z∗−z∗lnzz∗.
Obviously,Vx,y,z>0 for any x,y,z>0, except for the positive equilibrium E3=x∗,y∗,z∗.
By Lemma6, we have(48)DθVx,y,z≤1−x∗xDθxt+1+mx∗a21−y∗yDθyt+a31+mx∗a2a41−z∗zDθzt=x−x∗r1−a1x−y1+mx+1+mx∗a2y−y∗a2x1+mx−r2−a3z+a31+mx∗a2a4z−z∗a4y−r3=x−x∗a1x∗−a1x−y1+mx+y∗1+mx+1+mx∗a2y−y∗a2x1+mx−a2x∗1+mx∗−a3z+a3z∗+a31+mx∗a2a4y−y∗a4y−y∗=−a1x−x∗2+my∗x−x∗21+mx1+mx∗.
Sincea1>my∗/1+mx∗, then we have DθVx,y,z≤0. Thus, E3 is globally asymptotically stable.


## 4. Analysis of the Delayed Model

The conditions for nonnegativity boundedness, existence, and uniqueness derived for system (16) also apply to system (2). Systems (2) and (16) have identical equilibrium points. Due to the impact of time lags τ1 and τ2, the stability of system (2) needs to be rediscussed. Next, the stability and branch of system (2) are studied by selecting τ1 and τ2 as key parameters, and the critical bifurcation value is discussed precisely.


### 4.1. The Bifurcation of System (2) Caused by Delay τ1

In the following analysis, we focus on time delayτ1 as the bifurcation parameter of system (2) and obtain the critical value of Hopf bifurcation of the system.Making transformation,P1t=xt−x∗, P2t=yt−y∗, and P3t=zt−z∗. In consequence, system (2) is able to be transformed into(49)DθP1t=P1t+x∗r1−a1P1t+x∗−P2t−τ1+y∗1+mP1t+x∗,DθP2t=P2t+y∗a2P1t−τ2+x∗1+mP1t+x∗−r2−a3P3t−τ1+z∗,DθP3t=P3t+z∗a4P2t−τ2+y∗−r3.The linearized scheme from system (49) results in(50)DθP1t=b11P1t+b12P2t−τ1,DθP2t=b21P1t+b22P2t+b23P1t−τ2+b24P3t−τ1,DθP3t=b31P3t+b32P2t−τ2,where(51)b11=r1−2a1x∗−y∗1+mx∗2,b12=−x∗1+mx∗,b21=−a2mx∗y∗1+mx∗2,b22=a2x∗1+mx∗−r2−a3z∗,b23=a2y∗1+mx∗,b24=−a3y∗,b31=a4y∗−r3,b32=a4z∗.The characteristic equation of system (50) is as shown below:(52)U1s+U2se−sτ1=0,where(53)U1s=s3θ+−b31−b22−b11s2θ+b22b31+b11b31+b11b22sθ−b11b22b31,U2s=b12b21b31−b21b12sθ+b12b23b31+b11b32b24−b12b23sθ−b32b24sθe−sτ2.The real and imaginary parts ofUksk=1,2 are represented by Ukr and Uki. Suppose s is a purely imaginary root of (52), where s=ω1cosπ/2+isinπ/2ω1>0; it follows from (52) that(54)U2rcosω1τ1+U2isinω1τ1=−U1r,U2icosω1τ1−U2rsinω1τ1=−U1i.In view of (54), we derive that(55)cosω1τ1=−h1ω1h3ω1,sinω1τ1=−h2ω1h3ω1,where h1ω1=U1rU2r+U1iU2i, h2ω1=U1rU2i−U2rU1i, and h3ω1=U2r2+U2i2. It is apparent from (55) that(56)h12ω1+h22ω1−h32ω1=0.In terms ofcosω1τ1=−h1ω1/h3ω1, we obtain(57)τ10k=1ω1arccos−h1ω1h3ω1+2kπ,k=0,1,2….Suppose the equation of (56) has a positive real root ω10; we make(58)τ10=minτ10k,k=0,1,2…,where τ10k is provided by (57).Ifτ1=0, then (52) becomes(59)ς1s+ς2se−sτ2=0,where(60)ς1s=s3θ+−b31−b22−b11s2θ+b22b31+b11b31+b11b22−b12b21sθ+b12b21b31−b11b22b31,ς2s=−b12b23−b32b24sθ+b12b23b31+b11b32b24.Suppose thatςkr and ςki represent the real and imaginary parts of ςksk=1,2, s is a purely imaginary root of (59), and s=ω¯1cosπ/2+isinπ/2ω¯1>0; we can get that(61)ς2rcosω¯1τ2+ς2isinω¯1τ2=−ς1r,ς2icosω¯1τ2−ς2rsinω¯1τ2=−ς1i.Based on (61), we have(62)cosω¯1τ2=−k1ω¯1k3ω¯1,sinω¯1τ2=−k2ω¯1k3ω¯1,where k1ω¯1=ς1rς2r+ς1iς2i, k2ω¯1=ς1rς2i−ς2rς1i, and k3ω¯1=ς2r2+ς2i2. It is apparent from (62) that(63)k12ω¯1+k22ω¯1−k32ω¯1=0.In the light ofcosω¯1τ2=−k1ω¯1/k3ω¯1, we obtain(64)τ¯20k=1ω¯1arccos−k1ω¯1k3ω¯1+2kπ,k=0,1,2….Suppose the equation of (63) has a positive real root, we make(65)τ¯20=minτ¯20k,k=0,1,2…,where τ¯20k is provided by (64).Remark 1.
If equation (56) has no positive roots, then the system does not have bifurcation points. On the contrary, if equation (56) has more than one positive root, we take the minimum of all the roots. As mentioned above, τ10=minτ10k,k=0,1,2…. Similarly, τ¯20k is obtained this way.In order to better search for the criterion of the occurrence for bifurcation, the following hypotheses are helpful and essential:ℍ4:E^1F^1+E^2F^2/F^12+F^22>0, where E^1, E^2, F^1, and F^2 are described in the following.Lemma 8.
Letsτ1=ξτ1+iω1τ1 be the root of (17) near τ1=τ1j meeting ξτ1j=0 and ω1τ1j=ω10, so the following transversality criteria are true:(66)Redsdτ1|ω1=ω10,τ1=τ10>0,where ω10andτ10 are the critical frequency and the bifurcation point individually.Proof.
After differentiating equation (52) about τ1, we have(67)U1′sdsdτ1+U2′sdsdτ1e−sτ1+U2se−sτ1−τ1dsdτ1−s=0.
So, we can obtain(68)dsdτ1=E^sF^s,where(69)E^s=sU2se−sτ1,F^s=U1′s+U2′s−τ1U2se−sτ1.
LetE^1 and E^2 be the real and imaginary parts of E^s individually. F^1 and F^2 be the real and imaginary parts of F^s severally. After several algebraic calculation, we get from (68) that(70)Redsdτ1|ω1=ω10,τ1=τ10=E^1F^1+E^2F^2F^12+F^22,where(71)E^1=ω10U2rsinω10τ10−U2icosω10τ10,E^2=ω10U2rcosω10τ10+U2isinω10τ10,F^1=U1′r+U2′r−τ10U2rcosω10τ10+U2′i−τ10U2isinω10τ10,F^2=U1′i+U2′i−τ10U2icosω10τ10−U2′r−τ10U2rsinω10τ10.
As a result, supposeℍ4 implies that the transversality criteria are true. That is the proof of Lemma 8.With the support of Lemmas7 and 8, the under theorem can be derived.Theorem 7.
In the case ofℍ2,ℍ4, and τ2∈0,τ¯20, we have the following results:(1)
The positive equilibriumE3 of system (2) is asymptotically stable when τ1∈0,τ10.(2)
System (2) exhibits a Hopf bifurcation at E3 when τ1=τ10, i.e., it has a branch of periodic solution bifurcating from E3 near τ1=τ10


### 4.2. The Bifurcation of System (2) Caused by Delay τ2

In the following discussion, time delayτ2 is taken as the bifurcation parameter of system (2), and the Hopf bifurcation criterion of the system is obtained through theoretical analysis.The characteristic equation about system (2) is available:(72)V1s+V2se−sτ2=0,where(73)V1s=s3θ+−b31−b22−b11s2θ+b22b31+b11b31+b11b22sθ−b11b22b31+b12b21b31−b21b12sθe−sτ1,V2s=b12b23b31+b11b32b24−b12b23sθ−b32b24sθe−sτ1.The real and imaginary parts ofVksk=1,2 are represented by Vkr and Vki. Suppose s is a purely imaginary root of (72), where s=ω2cosπ/2+isinπ/2ω2>0, and we have(74)V2rcosω2τ2+V2isinω2τ2=−V1r,V2icosω2τ2−V2rsinω2τ2=−V1i.With the help of (74), we have(75)cosω2τ2=−l1ω2l3ω2,sinω2τ2=−l2ω2l3ω2,where l1ω2=V1rV2r+V1iV2i, l2ω2=V1rV2i−V2rV1i, and l3ω2=V2r2+V2i2. From (75), one has(76)l12ω2+l22ω2−l32ω2=0.In terms ofcosω2τ2=−l1ω2/l3ω2, we obtain(77)τ20k=1ω2arccos−l1ω2l3ω2+2kπ,k=0,1,2….Support the equation of (76) has a positive real root ω20, and we make(78)τ20=minτ20k,k=0,1,2…,where τ20k is provided by (77).Once eliminatingτ2 from (72), then(79)v1s+v2se−sτ1=0,where(80)v1s=s3θ+−b31−b22−b11s2θ+b22b31+b11b31+b11b22sθ−b11b22b31,v2s=b12b21b31+b12b23b31+b11b32b24+−b21b12−b12b23−b32b24sθ.Labelsvkr and vki represent the real and imaginary parts of vksk=1,2. Suppose s is a purely imaginary root of (79), where s=ω¯2cosπ/2+isinπ/2ω¯2>0, and we can derive from (79) that(81)v2rcosω¯2τ1+v2isinω¯2τ1=−v1r,v2icosω¯2τ1−v2rsinω¯2τ1=−v1i.By means of (81), we have(82)cosω¯2τ1=−q1ω¯2q3ω¯2,sinω¯2τ1=−q2ω¯2q3ω¯2,where q1ω¯2=v1rv2r+v1iv2i, q2ω¯2=v1rv2i−v2rv1i, and q3ω¯2=v2r2+v2i2. It is obtained from (82) that(83)q12ω¯2+q22ω¯2−q32ω¯2=0.Because ofcosω¯2τ1=−q1ω¯2/q3ω¯2, we obtain(84)τ¯10k=1ω¯2arccos−q1ω¯2q3ω¯2+2kπ,k=0,1,2….Suppose (83) has a positive real root; we make(85)τ¯10=minτ¯10k,k=0,1,2…,k=0,1,2…,where τ¯10k is provided by (84).In order to better search for the criterion of the occurrence for bifurcation, the following hypotheses are available and essential:ℍ5:ℵ1ℑ1+ℵ2ℑ2/ℑ12+ℑ22>0, where ℵ1, ℵ2, ℑ1, and ℑ2 are described in the following.Lemma 9.
Letsτ2=ξτ2+iω2τ2 be the root of (29) near τ2=τ2j meeting ξτ2j=0 and ω2τ2j=ω20, so the following transversality criteria are true:(86)Redsdτ2|ω2=ω20,τ2=τ20>0,where ω20andτ20 are the critical frequency and the bifurcation point individually.Proof.
Differentiating equation (72) with respect to τ2, we have(87)V1′sdsdτ2+V2′sdsdτ2e−sτ2+V2se−sτ2−τ2dsdτ2−s=0.
So, we can obtain(88)dsdτ2=ℵsℑs,where(89)ℵs=sV2se−sτ2,ℑs=V1′s+V2′s−τ2V2se−sτ2.
Defineℵ1 and ℵ2 be the real and imaginary parts of ℵs individually. ℑ1 and ℑ2 be the real and imaginary parts of ℑs individually. After several algebraic calculations, we receive from (88) that(90)Redsdτ2|ω2=ω20,τ2=τ20=ℵ1ℑ1+ℵ2ℑ2ℑ12+ℑ22,where(91)ℵ1=ω20V2rsinω20τ20−V2icosω20τ20,ℵ2=ω20V2rcosω20τ20+V2isinω20τ20,ℑ1=V1′r+V2′r−τ20V2rcosω20τ20+V2′i−τ20V2isinω20τ20,ℑ2=V1′i+V2′i−τ20V2icosω20τ20−V2′r−τ20V2rsinω20τ20.
As a result, supposeℍ5 implies the transversality criteria are true. That is the proof of Lemma 9.With the support of Lemmas7 and 8, the under theorem can be derived.Theorem 8.
In the case ofℍ2,ℍ5, and τ1∈0,τ¯10, we have the following results:(1)
The equilibriumE3 of system (2) is asymptotically stable when τ2∈0,τ20(2)
System (2) occurs a Hopf bifurcation at E3 when τ2=τ20, i.e., it has a branch of periodic solution bifurcating from E3 near τ2=τ20Remark 2.
In the previous work, many authors discussed Hopf bifurcations for fractional-order systems with single delay [39, 40], but in this study, we study Hopf bifurcations for fractional-order systems with two delays, which is of great significance for the discussion of Hopf bifurcations for systems with multiple delays.Remark 3.
In fact, the fractional-order system has a wider stability region than the integer order system. In other words, the fractional-order number will affect the stability of the system, taking the fractional-order number as the bifurcation parameter will also cause Hopf bifurcation.


## 4.1. The Bifurcation of System (2) Caused by Delay τ1

In the following analysis, we focus on time delayτ1 as the bifurcation parameter of system (2) and obtain the critical value of Hopf bifurcation of the system.Making transformation,P1t=xt−x∗, P2t=yt−y∗, and P3t=zt−z∗. In consequence, system (2) is able to be transformed into(49)DθP1t=P1t+x∗r1−a1P1t+x∗−P2t−τ1+y∗1+mP1t+x∗,DθP2t=P2t+y∗a2P1t−τ2+x∗1+mP1t+x∗−r2−a3P3t−τ1+z∗,DθP3t=P3t+z∗a4P2t−τ2+y∗−r3.The linearized scheme from system (49) results in(50)DθP1t=b11P1t+b12P2t−τ1,DθP2t=b21P1t+b22P2t+b23P1t−τ2+b24P3t−τ1,DθP3t=b31P3t+b32P2t−τ2,where(51)b11=r1−2a1x∗−y∗1+mx∗2,b12=−x∗1+mx∗,b21=−a2mx∗y∗1+mx∗2,b22=a2x∗1+mx∗−r2−a3z∗,b23=a2y∗1+mx∗,b24=−a3y∗,b31=a4y∗−r3,b32=a4z∗.The characteristic equation of system (50) is as shown below:(52)U1s+U2se−sτ1=0,where(53)U1s=s3θ+−b31−b22−b11s2θ+b22b31+b11b31+b11b22sθ−b11b22b31,U2s=b12b21b31−b21b12sθ+b12b23b31+b11b32b24−b12b23sθ−b32b24sθe−sτ2.The real and imaginary parts ofUksk=1,2 are represented by Ukr and Uki. Suppose s is a purely imaginary root of (52), where s=ω1cosπ/2+isinπ/2ω1>0; it follows from (52) that(54)U2rcosω1τ1+U2isinω1τ1=−U1r,U2icosω1τ1−U2rsinω1τ1=−U1i.In view of (54), we derive that(55)cosω1τ1=−h1ω1h3ω1,sinω1τ1=−h2ω1h3ω1,where h1ω1=U1rU2r+U1iU2i, h2ω1=U1rU2i−U2rU1i, and h3ω1=U2r2+U2i2. It is apparent from (55) that(56)h12ω1+h22ω1−h32ω1=0.In terms ofcosω1τ1=−h1ω1/h3ω1, we obtain(57)τ10k=1ω1arccos−h1ω1h3ω1+2kπ,k=0,1,2….Suppose the equation of (56) has a positive real root ω10; we make(58)τ10=minτ10k,k=0,1,2…,where τ10k is provided by (57).Ifτ1=0, then (52) becomes(59)ς1s+ς2se−sτ2=0,where(60)ς1s=s3θ+−b31−b22−b11s2θ+b22b31+b11b31+b11b22−b12b21sθ+b12b21b31−b11b22b31,ς2s=−b12b23−b32b24sθ+b12b23b31+b11b32b24.Suppose thatςkr and ςki represent the real and imaginary parts of ςksk=1,2, s is a purely imaginary root of (59), and s=ω¯1cosπ/2+isinπ/2ω¯1>0; we can get that(61)ς2rcosω¯1τ2+ς2isinω¯1τ2=−ς1r,ς2icosω¯1τ2−ς2rsinω¯1τ2=−ς1i.Based on (61), we have(62)cosω¯1τ2=−k1ω¯1k3ω¯1,sinω¯1τ2=−k2ω¯1k3ω¯1,where k1ω¯1=ς1rς2r+ς1iς2i, k2ω¯1=ς1rς2i−ς2rς1i, and k3ω¯1=ς2r2+ς2i2. It is apparent from (62) that(63)k12ω¯1+k22ω¯1−k32ω¯1=0.In the light ofcosω¯1τ2=−k1ω¯1/k3ω¯1, we obtain(64)τ¯20k=1ω¯1arccos−k1ω¯1k3ω¯1+2kπ,k=0,1,2….Suppose the equation of (63) has a positive real root, we make(65)τ¯20=minτ¯20k,k=0,1,2…,where τ¯20k is provided by (64).Remark 1.
If equation (56) has no positive roots, then the system does not have bifurcation points. On the contrary, if equation (56) has more than one positive root, we take the minimum of all the roots. As mentioned above, τ10=minτ10k,k=0,1,2…. Similarly, τ¯20k is obtained this way.In order to better search for the criterion of the occurrence for bifurcation, the following hypotheses are helpful and essential:ℍ4:E^1F^1+E^2F^2/F^12+F^22>0, where E^1, E^2, F^1, and F^2 are described in the following.Lemma 8.
Letsτ1=ξτ1+iω1τ1 be the root of (17) near τ1=τ1j meeting ξτ1j=0 and ω1τ1j=ω10, so the following transversality criteria are true:(66)Redsdτ1|ω1=ω10,τ1=τ10>0,where ω10andτ10 are the critical frequency and the bifurcation point individually.Proof.
After differentiating equation (52) about τ1, we have(67)U1′sdsdτ1+U2′sdsdτ1e−sτ1+U2se−sτ1−τ1dsdτ1−s=0.
So, we can obtain(68)dsdτ1=E^sF^s,where(69)E^s=sU2se−sτ1,F^s=U1′s+U2′s−τ1U2se−sτ1.
LetE^1 and E^2 be the real and imaginary parts of E^s individually. F^1 and F^2 be the real and imaginary parts of F^s severally. After several algebraic calculation, we get from (68) that(70)Redsdτ1|ω1=ω10,τ1=τ10=E^1F^1+E^2F^2F^12+F^22,where(71)E^1=ω10U2rsinω10τ10−U2icosω10τ10,E^2=ω10U2rcosω10τ10+U2isinω10τ10,F^1=U1′r+U2′r−τ10U2rcosω10τ10+U2′i−τ10U2isinω10τ10,F^2=U1′i+U2′i−τ10U2icosω10τ10−U2′r−τ10U2rsinω10τ10.
As a result, supposeℍ4 implies that the transversality criteria are true. That is the proof of Lemma 8.With the support of Lemmas7 and 8, the under theorem can be derived.Theorem 7.
In the case ofℍ2,ℍ4, and τ2∈0,τ¯20, we have the following results:(1)
The positive equilibriumE3 of system (2) is asymptotically stable when τ1∈0,τ10.(2)
System (2) exhibits a Hopf bifurcation at E3 when τ1=τ10, i.e., it has a branch of periodic solution bifurcating from E3 near τ1=τ10


## 4.2. The Bifurcation of System (2) Caused by Delay τ2

In the following discussion, time delayτ2 is taken as the bifurcation parameter of system (2), and the Hopf bifurcation criterion of the system is obtained through theoretical analysis.The characteristic equation about system (2) is available:(72)V1s+V2se−sτ2=0,where(73)V1s=s3θ+−b31−b22−b11s2θ+b22b31+b11b31+b11b22sθ−b11b22b31+b12b21b31−b21b12sθe−sτ1,V2s=b12b23b31+b11b32b24−b12b23sθ−b32b24sθe−sτ1.The real and imaginary parts ofVksk=1,2 are represented by Vkr and Vki. Suppose s is a purely imaginary root of (72), where s=ω2cosπ/2+isinπ/2ω2>0, and we have(74)V2rcosω2τ2+V2isinω2τ2=−V1r,V2icosω2τ2−V2rsinω2τ2=−V1i.With the help of (74), we have(75)cosω2τ2=−l1ω2l3ω2,sinω2τ2=−l2ω2l3ω2,where l1ω2=V1rV2r+V1iV2i, l2ω2=V1rV2i−V2rV1i, and l3ω2=V2r2+V2i2. From (75), one has(76)l12ω2+l22ω2−l32ω2=0.In terms ofcosω2τ2=−l1ω2/l3ω2, we obtain(77)τ20k=1ω2arccos−l1ω2l3ω2+2kπ,k=0,1,2….Support the equation of (76) has a positive real root ω20, and we make(78)τ20=minτ20k,k=0,1,2…,where τ20k is provided by (77).Once eliminatingτ2 from (72), then(79)v1s+v2se−sτ1=0,where(80)v1s=s3θ+−b31−b22−b11s2θ+b22b31+b11b31+b11b22sθ−b11b22b31,v2s=b12b21b31+b12b23b31+b11b32b24+−b21b12−b12b23−b32b24sθ.Labelsvkr and vki represent the real and imaginary parts of vksk=1,2. Suppose s is a purely imaginary root of (79), where s=ω¯2cosπ/2+isinπ/2ω¯2>0, and we can derive from (79) that(81)v2rcosω¯2τ1+v2isinω¯2τ1=−v1r,v2icosω¯2τ1−v2rsinω¯2τ1=−v1i.By means of (81), we have(82)cosω¯2τ1=−q1ω¯2q3ω¯2,sinω¯2τ1=−q2ω¯2q3ω¯2,where q1ω¯2=v1rv2r+v1iv2i, q2ω¯2=v1rv2i−v2rv1i, and q3ω¯2=v2r2+v2i2. It is obtained from (82) that(83)q12ω¯2+q22ω¯2−q32ω¯2=0.Because ofcosω¯2τ1=−q1ω¯2/q3ω¯2, we obtain(84)τ¯10k=1ω¯2arccos−q1ω¯2q3ω¯2+2kπ,k=0,1,2….Suppose (83) has a positive real root; we make(85)τ¯10=minτ¯10k,k=0,1,2…,k=0,1,2…,where τ¯10k is provided by (84).In order to better search for the criterion of the occurrence for bifurcation, the following hypotheses are available and essential:ℍ5:ℵ1ℑ1+ℵ2ℑ2/ℑ12+ℑ22>0, where ℵ1, ℵ2, ℑ1, and ℑ2 are described in the following.Lemma 9.
Letsτ2=ξτ2+iω2τ2 be the root of (29) near τ2=τ2j meeting ξτ2j=0 and ω2τ2j=ω20, so the following transversality criteria are true:(86)Redsdτ2|ω2=ω20,τ2=τ20>0,where ω20andτ20 are the critical frequency and the bifurcation point individually.Proof.
Differentiating equation (72) with respect to τ2, we have(87)V1′sdsdτ2+V2′sdsdτ2e−sτ2+V2se−sτ2−τ2dsdτ2−s=0.
So, we can obtain(88)dsdτ2=ℵsℑs,where(89)ℵs=sV2se−sτ2,ℑs=V1′s+V2′s−τ2V2se−sτ2.
Defineℵ1 and ℵ2 be the real and imaginary parts of ℵs individually. ℑ1 and ℑ2 be the real and imaginary parts of ℑs individually. After several algebraic calculations, we receive from (88) that(90)Redsdτ2|ω2=ω20,τ2=τ20=ℵ1ℑ1+ℵ2ℑ2ℑ12+ℑ22,where(91)ℵ1=ω20V2rsinω20τ20−V2icosω20τ20,ℵ2=ω20V2rcosω20τ20+V2isinω20τ20,ℑ1=V1′r+V2′r−τ20V2rcosω20τ20+V2′i−τ20V2isinω20τ20,ℑ2=V1′i+V2′i−τ20V2icosω20τ20−V2′r−τ20V2rsinω20τ20.
As a result, supposeℍ5 implies the transversality criteria are true. That is the proof of Lemma 9.With the support of Lemmas7 and 8, the under theorem can be derived.Theorem 8.
In the case ofℍ2,ℍ5, and τ1∈0,τ¯10, we have the following results:(1)
The equilibriumE3 of system (2) is asymptotically stable when τ2∈0,τ20(2)
System (2) occurs a Hopf bifurcation at E3 when τ2=τ20, i.e., it has a branch of periodic solution bifurcating from E3 near τ2=τ20Remark 2.
In the previous work, many authors discussed Hopf bifurcations for fractional-order systems with single delay [39, 40], but in this study, we study Hopf bifurcations for fractional-order systems with two delays, which is of great significance for the discussion of Hopf bifurcations for systems with multiple delays.Remark 3.
In fact, the fractional-order system has a wider stability region than the integer order system. In other words, the fractional-order number will affect the stability of the system, taking the fractional-order number as the bifurcation parameter will also cause Hopf bifurcation.


## 5. Numerical Results

Now, we verify our theoretical consequences by numerical simulations. For practical reasons, we only analyze the positive equilibriumE3, rather than E0, E1, or E2. We set r1=1.5, r2=0.125, r3=0.25, a1=1.2, a2=0.4, a3=0.6, a4=0.4, and m=0.5; then, system (2) can be transformed into(92)Dθxt=xt1.5−1.2xt−yt−τ11+0.5xt,Dθyt=yt0.4xt−τ21+0.5xt−0.125−0.6zt−τ1,Dθzt=zt0.4yt−τ2−0.25.We can easily verify that the system only has a positive equilibrium,E3=0.8895,0.6250,0.2021.Case 1.
We setτ1=τ2=0,θ=0.8,0.9,1, and then, we can verify that the system meets the condition of Theorem 6; Figure 1 indicates the positive equilibrium E3 of system (92) is stable. Particularly, if θ=0.9, Figure 2 shows that E3 of system (92) is globally asymptotically stable.
To discuss the bifurcation points about system (92), let us define θ=0.9.Figure 1
Case1: the trajectories and phase diagrams of system (92) over time at τ1=τ2=0.Figure 2
Case1: the phase diagrams for system (92) at τ1=τ2=0 and θ=0.9.Case 2.
We fixτ2 and choose τ1 as the branch parameter to consider the bifurcation of system (92). If τ1=0, we can figure out that τ¯20=3.6581. In Figure 3, it implies system (92) is asymptotically stable when τ1=0 and τ2=3, and we can clearly find that the system is not stable and Hopf bifurcation appears when τ1=0 and τ2=4. Then, we can calculate that τ10=0.8649 when τ2=2.5. Based on Theorem 7, system (92) is asymptotically stable if τ1=0.3<0.8649 and τ2=2.5. However, if we increase τ1 from 0.3 to 0.9>0.8649, system (92) is unstable and Hopf bifurcation occurs, see Figure 4.Figure 3
Case2: the trajectories and phase diagrams of system (92) over time at τ1=0.Figure 4
Case2: the trajectories and phase diagrams of system (92) over time at τ2=2.5.Case 3.
We fixτ1 and choose τ2 as the branch parameter to consider the bifurcation of system (92). When τ2=0, we can calculate that τ¯10=4.0646. In Figure 5, it implies system (92) is asymptotically stable when τ1=3.2 and τ2=0, and we can clearly find that system (92) is not stable and Hopf bifurcation appears when τ1=4.2 and τ2=0. Then, we can get that τ20=3.5566 when τ1=2. By Theorem 8, system (92) is asymptotically stable when τ1=2 and τ2=1<3.5566. However, if we increase τ2 from 1 to 3.6>3.5566, system (92) is unstable and Hopf bifurcation occurs, see Figure 6.
In order to study the impact of fractional order on bifurcation points, we make the following simulation results.Figure 5
Case3: the trajectories and phase diagrams of system (92) over time at τ2=0.Figure 6
Case3: the trajectories and phase diagrams of system (92) over time at τ1=2.Case 4.
We chooseτ2=0.4. When θ=1, we can calculate τ10=1.5156. At this time, we take τ1=2>1.5156, and Hopf bifurcation appears in system (92). When θ=0.9, we can calculate τ10=3.5790; then, we take τ1=2<3.5790; system (92) is asymptotically stable. Look at Figure 7. Figure 8 shows the impact of fractional order on τ1.Figure 7
Case4: the trajectories and phase diagrams of system (92) at τ1=2andτ2=0.4.Figure 8
Case4: effect of θ on bifurcation point τ10.Case 5.
Similarly, we takeτ1=0.2. When θ=1, we can obtain τ20=3.2996. At this time, we take τ2=3.3>3.2996, and Hopf bifurcation appears in system (92). When θ=0.9, we can get τ20=5.3302, and we take τ2=3.3<5.3302; system (92) is asymptotically stable (see Figure 9). Figure 10 shows the impact of fractional order on τ2.Figure 9
Case5: the trajectories and phase diagrams of system (92) at τ1=0.2andτ2=3.3.Figure 10
Case5: effect of θ on bifurcation point τ20.


## 6. Conclusion

In this study, a fractional-order food chain system involving two time delays has been presented. Nonnegative, bounded, existence, and uniqueness about the solution of the system have been proved. For nondelay system, we have discussed the local stability of the system equilibrium point and proved the globally asymptotically stability of the positive equilibrium point by constructing Lyapunov functions. By using time delays as parameters to discuss the Hopf bifurcation, which has showed that when the delay exceeds the critical value, the Hopf bifurcation will appear in the system, that is to say, the system will change from stable to unstable and a periodic solution will appear. In particular, the periodic oscillation behavior of the system could be suppressed by fractional order, which has indicated that the fractional-order system has a larger range of stability region than the integer-order system.



---

*Source: 5313931-2022-06-15.xml*