[BACK]

Deterministic and Stochastic Fractional-Order Hastings-Powell Food Chain Model

1Department of Mathematics, Unaizah College of Sciences and Arts, Qassim University, Unaizah, 51911, Saudi Arabia
2Department of Mathematics, Buraydah College of Sciences and Arts, Qassim University, Buraydah, 51431, Saudi Arabia
*Corresponding Author: Moustafa El-Shahed. Email: elshahedm@yahoo.com
Received: 09 April 2021; Accepted: 11 June 2021

Abstract: In this paper, a deterministic and stochastic fractional-order model of the tri-trophic food chain model incorporating harvesting is proposed and analysed. The interaction between prey, middle predator and top predator population is investigated. In order to clarify the characteristics of the proposed model, the analysis of existence, uniqueness, non-negativity and boundedness of the solutions of the proposed model are examined. Some sufficient conditions that ensure the local and global stability of equilibrium points are obtained. By using stability analysis of the fractional-order system, it is proved that if the basic reproduction number R0<1, the predator free equilibrium point E1 is globally asymptotically stable. The occurrence of local bifurcation near the equilibrium points is investigated with the help of Sotomayor’s theorem. Some numerical examples are given to illustrate the theoretical findings. The impact of harvesting on prey and the middle predator is studied. We conclude that harvesting parameters can control the dynamics of the middle predator. A numerical approximation method is developed for the proposed stochastic fractional-order model.

Keywords: Food chain model; global dynamics; stability; stochastic; Hastings-Powell model; harvesting

1  Introduction

Mathematical analysis is one of the important tools for understanding and interpreting different interactions in the environment around us. The food chain model system is attractive to researchers in theoretical ecology because it helps to understand the relationships between populations and describe the behavior of the ecosystem. Hastings et al. [1] investigated a three-species food chain model with Holling type II functional responses. Hastings-Powell model, has been revisited recently by many authors [216]. In the real world, food chain models are always affected by environmental noise. Thus, the stochastic models may be a more appropriate way of modeling the Hastings-Powell food chain model in many circumstances [17,18]. Recently, fractional calculus has been applied to describe different mathematical models, and it has been shown to be more accurate in some cases compared to the classical models [1926]. The main objective of this paper is to study the deterministic and stochastic fractional-order Hastings-Powell model incorporating harvesting. This model is an extension of the classical Hastings-Powell model (1) by introducing fractional derivative, stochastic and harvesting. The qualitative behavior of the proposed model is investigated. Also, a numerical approximation method is developed for the proposed stochastic fractional-order model.

The paper is arranged as follows: In Section 2, the mathematical model is described. Some preliminary results, such as existence, uniqueness, nonnegativity and boundedness are presented in Section 3. The local and global stability of equilibrium points of the fractional-order food chain model is analyzed in Section 4. With the help of Sotomayor’s theorem, the transcritical bifurcation of the proposed model is investigated in Section 5. Section 6 extends the deterministic fractional- order food chain model to the stochastic fractional-order model. In Section 7, some numerical simulations are presented to verify the obtained theoretical results. Finally, the conclusions are given in Section 8.

2  Mathematical Model

Recently, Nath et al. [27] studied the following system of integer order differential equations to understand the underlying dynamics of the food chain model:

dXdT=rX(1Xk)βXYA+XH1X,dYdT=e1βXYA+XγYZB+YdYH2Y2,dZdT=e2γYZB+YδZH3Z, (1)

with initial values X(0)=X00,Y(0)=Y00,Z(0)=Z00.

In the above model X(T), Y(T) and Z(T) represent the population density of prey, middle predator and top predator at time T, respectively. The prey grows with an intrinsic growth rate r and carrying capacity k in the absence of predation. It is assumed that the middle predator consumes prey following Holling type II functional response, and β is the maximum predation rate of the middle predator in prey. The middle predator is attacked by the top predator at rate γ using Holling type II functional response. A and B are the half saturation constants. e1 and e2 are conversion rate for middle predator and top predator, respectively. d and δ are the mortality rate of the middle predator and top predator, respectively. Here H1 and H3 are the harvesting rate of the prey population and top predator, respectively. The term H2Y2 represents the quadratic harvesting of the middle predator.

Introducing dimensionless

x=Xk,y=Yke1,z=Zke1e2,t=rT.

Making the above substitutions in the model (1). Then the system yields the following form

dxdt=x(1x)a1xy1+b1xhx,dydt=a1xy1+b1xa2yz1+b2yμ1yηy2,A=πr2dzdt=a2yz1+b2yμ2z. (2)

The dimensionless parameters in the food chain model (2) are defined as

a1=kβe1rA,a2=kγe1e2rB,b1=kA,b2=ke1B,h=H1r,μ1=dr,μ2=δ+H3r,η=ke1H2r

Following [13,15,16], by replacing the integer derivative by a Caputo fractional derivative in (2), one can obtain the following fractional-order system

cDαx=x(γx)a1xy1+b1x,cDαy=a1xy1+b1xa2yz1+b2yμ1yηy2,cDαz=a2yz1+b2yμ2z, (3)

where, 0<α<1 and γ=1h.

3  Some Preliminary Results

3.1 Existence and Uniqueness

In this section, the existence and uniqueness of the solutions of the fractional-order system (3) are investigate in the region Ω×(0,T] where

Ω={(x,y,z)R+3:max(|x|,|y|,|z|)φ},

for sufficiently large φ.

Theorem 1. For each X0=(x0,y0,z0)Ω, there exists a unique solution X(t)Ω of the fractional-order system (3), which is defined for all t0.

Proof. Define a mapping F(X)=(F1(X),F2(X),F3(X)), in which

F1(X)=x(γx)a1xy1+b1x,F2(X)=a1xy1+b1xa2yz1+b2yμ1yηy2,F3(X)=a2yz1+b2yμ2z. (4)

For any X,X¯Ω, it follows from (3) that

F(X)F(X¯)=|F1(X)F1(X¯)|+|F2(X)F2(X¯)|+|F3(X)F3(X¯)|=|x(γx)a1xy1+b1xx¯(γx¯)+a1x¯y¯1+b1x¯|+|a1xy1+b1xa2yz1+b2yμ1yηy2a1x¯y¯1+b1x¯+a2y¯z¯1+b2y¯+μ1y¯+ηy¯2|+|a2yz1+b2yμ2za2y¯z¯1+b2y¯+μ2z¯|(γ+2φ+2a1φ)|xx¯|+(μ1+2ηφ+2a1φ+2a1b1+2a2φ)|yy¯|+(μ2+4a2b2)|zz¯|M1XX¯,

where

M1=max{γ+2φ(1+a1),μ1+2ηφ+2a1φ+2a1b1+2a2φ,μ2+4a2b2}.

Hence, F(X) satisfies the Lipschitz condition with respect to X. According to Cresson et al. [28], as F(X) locally Lipschitz. Then there exists a unique local solution to the fractional-order system (3).

3.2 Non-Negativity and Boundedness

The following results show the non-negativity of the solutions of the fractional-order system (3). According to [28], a model of the form dXdt=F(X) satisfies the positivity property if and only if for all i=1,2,3, Fi(X)0 for all XR+3 such that Xi=0. Thus, the solution of the integer-order model (1), with nonnegative initial conditions remains nonnegative. Also, the solution satisfies the Lipschitz condition, as stated in Theorem 1. By Theorem 5 and Theorem 6 in Cresson et al. [28], the solution of the fractional-order model (3) also satisfies the non-negativity. The boundedness of the solutions of model (3) are given in the following theorem.

Theorem 2. All the solutions of the fractional-order food chain Model (3) starting in R+3 are uniformly bounded.

Proof. The approach of [29,30] is utilized. Let (x(t),y(t),z(t)) to be any solution of the system (3) with non-negative initial conditions. Let w(t)=x(t)+y(t)+z(t), then

cDαw(t)x2+γxμ1yμ2zx2+2γxν(x+y+z),

where, ν<min{γ,μ1,μ2}, thus, cDαw(t)+νwγ2. In accordance with Lemma 9 in Choi et al. [31], it follows that

0w(t)w(0)Eα(νtα)+γ2tαEα,α+1(νtα),

where Eα is the Mittag-Leffler function. According to Lemma 5 and Corollary 6 in Choi et al. [31], it follows

0w(t)γ2ν,as t.

Hence all the solutions of fractional-order a tri-trophic food chain model (3) that start in R+3 are uniformly bounded in the region

H={(x,y,z)R+3:w(t)γ2ν+ξ,for anyξ>0}.

One can also prove that xγ.

4  Equilibria and Stability

The fractional-order system (3) has the following equilibrium points:

1) E0=(0,0,0), that is, the extinction of prey, middle predator, and top predator. The trivial equilibrium E0 always exists.

2) The middle predator and top predator free equilibrium point E1=(γ,0,0). Therefore E1 exists if h<1. Using the next generation method, one can obtain the basic reproduction number

R0=a1γμ1(1+b1γ).

3) Following [32], R0 determine the local and global stability of E1. Here γa1(1+γb1) and 1μ1 are the birth rate and the mean lifespan of middle predator at E1, respectively. Subsequently their product gives the basic reproduction number at E1. In the following, two critical parameters R1 and R2, can be used to classify the dynamics of the fractional-order model (3). The threshold parameter R1 defined by R1=a1x2μ1(1+b1x2), while the threshold parameter R2 defined by R2=a2y2μ2(1+b2y2).

4) The top predator extinction equilibrium point E2=(x2,y2,0), where y2=μ1η(R11), and x2 satisfies the equation

b1x22+(1γb1)x2+a1μ1η(R11)γ=0.

5) The top predator extinction equilibrium point E2 exists if b1γ>1 and 1+γηa1μ1<R1<1+η(1+b1γ)24a1b1μ1.

6) The coexistence positive equilibrium point E3=(x3,y3,z3), where

y3=μ2a2μ2b2,z3=(1+b2y3)a2(a1x31+b1x3μ1ηy3),

7) and x3 is the positive root of the equation:

b1x32+(1γb1)x3+(a1y3γ)=0. (5)

E3(x3,y3,z3) exists if

a2>μ2b2,a1x31+b1x3>μ1+ηy3,b1γ>1,a1y3>γ,(1γb1)2>4b1(a1y3γ).

The locally and globally asymptotically stable of equilibrium points of fractional-order food chain model (3) are now investigated. The Jacobian matrix is given as follows:

J(x,y,z)=(γ2xa1y(1+b1x)2a1x1+b1x0a1y(1+b1x)2a1x1+b1xa2z(1+b2y)2μ12ηya2y1+b2y0a2z(1+b2y)2a2y1+b2yμ2). (6)

The eigenvalues of J around E0 are γ,μ1 and μ2, therefore E0 is unstable saddle point. The stability of middle predator extinction equilibrium point E1 is investigated as follows

Theorem 3. If R0<1, then E1 is locally asymptotically stable.

Proof. The Jacobian matrix of model (3) at E1 is

J(E1)=(γa1γ1+b1γ00a1γ1+b1γμ1000μ2). (7)

The eigenvalues of J(E1) are γ,μ2 and a1γ1+b1γμ1. Thus, E1 is locally asymptotically stable if a1γ1+b1γ<μ1, which is equivalent to R0<1.

Theorem 4. If R0<1, then E1 is globally asymptotically stable.

Proof. Let us consider the following positive definite Lyapunov function

V1=L12(xγ)2+12y2+y+z.

The time derivative of V1 along the solution of fractional-order tri-trophic food chain model (3), one obtains

cDαV1L1(xγ)dxdt+(1+y)dydt+dzdtL1(xγ)(x(γx)a1xy1+b1x)+(1+y)[a1xy1+b1xa2yz1+b2yμ1yηy2]+(a2yz1+b2yμ2z)L1x(xγ)2+(a1x1+b1xμ1)y2+(L1a1γx1+b1x+a1x1+b1xμ1)y.

According to Theorem 1, xγ, and consequently a1x1+b1xa1γ1+b1γ. The above inequality can be written as

cDαV1(a1γ1+b1γμ1)y2+(L1a1γ21+b1γ+a1γ1+b1γμ1)y.

Choosing L1=1R0γR0, one obtains cDαV10.

According to generalized Lyapunov–Lasalle’s invariance principle [33], E1 is globally asymptotically stable.

Hence the equilibrium point E1 is globally asymptotically stable if R0<1.

The global stability of the equilibrium point E1 can be further demonstrated from the second equation of the system (3) as follows

cDαy(a1x1+b1xμ1)y(a1γ1+b1γμ1)yμ1(R01)y,

when R0<1, y(t). The extinction of middle predator y(t) implies the extinction of top predator z(t). Furthermore, the equilibrium point E1 is globally asymptotically stable. Following [32], R0<1, implies that the conversion of prey into middle predator during lifespan of middle predator is less than unity. As a result, the biomass flow from prey to middle predator and top predator will stop. In this case, the middle predator and top predator will be extinct.

The stability of top predator extinction equilibrium point E2 is investigated as follows

Theorem 5. If R2<1 and a1b1y2(1+b1x2)2<1 then the top predator extinction equilibrium point E2 is locally asymptotically stable.

Proof. The Jacobian matrix of model (3) at E2 is

J(E2)=(Υa1x21+b1x20a1y2(1+b1x2)2ηy2a2y21+b2y200a2y21+b2y2μ2), (8)

where Υ=γ2x2a1y2(1+b1x2)2. The eigenvalues of J(E2) are the roots of equation

(a2y21+b2y2μ2λ)(λ2+(ηy2Υ)λ+a12x2y2(1+b1x2)3ηy2Υ)=0.

The above equation has the following eigenvalue λ1=a2y21+b2y2μ2=μ2(R21). The eigenvalues λ2,3 have negative real parts if Υ<0. Thus, if R2<1 and Υ<0, then the top predator extinction equilibrium point E2 is locally asymptotically stable.

Theorem 6. If a1b1y2(1+b1x2)<1 and a2y2<μ2, then the top predator extinction equilibrium point E2 is globally asymptotically stable.

Proof. The following positive definite Lyapunov function is considered.

V2=L2(xx2x2lnxx2)+yy2y2lnyy2+z.

By calculating the time derivative of V2 along the solution of system (3), one obtains,

cDαV2L2(xx2)[γxa1y1+b1x]+(yy2)[a1x1+b1xa2z1+b2yμ1ηy]+z(a2y1+b2yμ2)L2(xx2)2(a1b1y2(1+b1x)(1+b1x2)1)+[a1a1L2a1b1L2x2](1+b1x)(1+b1x1)(xx2)(yy2)+z(a2y2μ2).

Choosing L2=11+b1x2, then cDαV20 when a1b1y2(1+b1x2)<1 and a2y2<μ2, and hence, according to generalized Lyapunov–Lasalle’s invariance principle [33], the top predator extinction equilibrium point E2 is globally asymptotically stable when a1b1y2(1+b1x2)<1 and a2y2<μ2.

The stability of coexistence equilibrium point E3 is investigated as follows.

Theorem 7. If a1b1y3(1+b1x3)2<1 and a2b2z3(1+b2y3)2<η, then the coexistence equilibrium point E3 is locally asymptotically stable.

Proof. The Jacobian matrix of model (3) at E3 is

J(E3)=(C1a1x31+b1x30a1y3(1+b1x3)2C2a2y31+b2y30a2z3(1+b2y3)20),

where C1=a1b1x3y3(1+b1x3)2x3 and C2=a2b2y3z3(1+b2y3)2ηy3. The characteristic polynomial of J(E3) is

λ3+ϕ1λ2+ϕ2λ+ϕ3=0, (9)

where

ϕ1=(C1+C2),ϕ2=[a12x3y3(1+b1x3)3+a22y3z3(1+b2y3)3+C1C2],ϕ3=a22C1y3z3(1+b2y3)3.

If we choose C1<0 and C2<0, then ϕ1>0, ϕ2>0, ϕ3>0 and ϕ1ϕ2>ϕ3. According to Routh-Hurtwitz criteria, the system (3) is locally asymptotically stable around the coexistence equilibrium point E3, when a1b1y3(1+b1x3)2<1 and a2b2z3(1+b2y3)2<η.

Following [34], it is interested to note that the function

Φi(α,θ)=απ2|arg(λi(θ))|,i=1,2,3, (10)

has a similar effect as the real part of the eigenvalue in the integer order system. If Φi(α,θ)<0 for all i=1,2,3, then E3 is locally asymptotically stable. If there exist i such that Φi(α,θ)>0, then E3 is unstable.

Theorem 8. If a1b1y31+b1x3<1 and a2b2z31+b1y3<η, then the coexistence equilibrium point E3 is globally asymptotically stable.

Proof. Let us consider the function

V3=L3(xx3x3lnxx3)+(yy3y3lnyy3)+L4(zz3z3lnzz3).

By calculating the time derivative of V3 along along with the solution of system (3), one obtains,

cDαV3L3(xx3)[γxa1y1+b1x]+(yy3)[a1x1+b1xa2z1+b2yμ1ηy]+L4(zz3)[a2y1+b2yμ2]=L3(xx3)[(x3x)+a1y31+b1x3a1y1+b1x]+L4(zz3)[a2y1+b2ya2y31+b2y3]+(yy3)[a1x1+b1xa1x31+b1x3+a2z31+b2y3a2z1+b2y+η(y3y)]=L3(a1b1y3(1+b1x)(1+b1x3)1)(xx3)2+(xx3)(yy3)(1+b1x)(1+b1x3)[a1a1L3a1b1L3x3]+(yy3)(zz3)(1+b2y)(1+b2y3)[a2L4a2a2b2y3]+(a2b2z3(1+b2y)(1+b2y3)η)(yy3)2.

Choosing L3=11+b1x3 and L4=1+b2y3, then cDαV30 when a1b1y31+b1x3<1 and a2b2z31+b2y3<η. Thus, according to generalized Lyapunov–Lasalle’s invariance principle [33], the coexistence equilibrium point E3 is globally asymptotically stable.

5  Bifurcation Analysis

In this section we will investigate the local bifurcation near the free predator equilibrium point of the food chain model (2) with the help of Sotomayor’s theorem [35] to discuss the bifurcation analysis of the underlying system. The food chain model (2) can be rewritten in a vector form dXdt=F(X), where X=(x,y,z)T and F=(F1,F2,F3)T with Fi,i=1,2,3 are given in (4).

Theorem 9. The food chain model (2) undergoes a transcritical bifurcation with respect to the bifurcation parameter μ1 around the free predator equilibrium point E1=(γ,0,0) if R0=1.

Proof. The Jacobian matrix of the food chain model (2) at the free predator equilibrium point E1 with μ1=μ1=γa11+b1γ has zero eigenvalue takes the form

J(E1)=(γγa11+b1γ000000μ2).

The eigenvector corresponding to J(E1)W1=0 is W1=(ν1,(1+b1γ)a1,0)T where ν1 is any non zero real number. Similarly, the eigenvector corresponding to J(E1)TW2=0 is given by W2=(0,τ2,0)T, where τ2 is any non zero real number. Consider, Fμ1=Fμ1(X,μ1)=(0,y,0)T, thus, W2TFμ1(E1,μ1)=0. Therefore, according to Sotomayor’s theorem for local bifurcation, the food chain model (2) has no saddle-node bifurcation near E1 at μ1=γa11+b1γ. Now,

DFμ1(E1,μ1)=(000010000),

then W2TDFμ1(E1,μ1)W1=ν1τ2(b1γ+1)a10, and

W2TD2F(X,μ1)(W1,W1)=τ2(2a1ν1ν2(b1γ+1)22a2ν3ν22ην22)0.

Thus, according to Sotomayor’s theorem, the food chain model (2) has a transcritical bifurcation at μ1=γa11+b1γ as the parameter μ1 passes through the value μ1, thus the proof is complete. It is interested to note that R0=1 is equivalent to μ1=γa11+b1γ.

6  Stochastic Fractional-order Model

This section extends the deterministic fractional-order food chain model (3) to the following stochastic fractional-order model.

CDαx=x(γx)a1xy1+b1x+σ1xdσ1dt,CDαy=a1xy1+b1xa2yz1+b2yμ1yηy2+σ2ydσ2dt,CDαz=a2yz1+b2yμ2z+σ3zdσ3dt, (11)

where σi(i=1,2,3) are independent standard Brownian motions with σi(0)=0 and σi>0 denote the intensities of the white noise. The stochastic fractional-order food chain model (11) can be written in the general form:

CDαX(t)=F(X)+g(X)dσdt, (12)

where F(X) is given in (4), g(x)=(σ1x,σ2y,σ3z) and dσdt=(dσ1dt,dσ2dt,dσ3dt)T. Applying Riemann–Liouville integral to both sides of (12), one can obtain the following stochastic Volterra integral equation.

X(t)=X0+0tF(X)(ts)α1Γ(α)ds+0tg(X)(ts)α1dσ(s)Γ(α)ds. (13)

According to [3639], under some conditions on the coefficient functions, the global existence and uniqueness of solutions for the stochastic fractional-order system (13) can be investigated. Because Grunwald-Letnikov’s definition is the most straightforward from the point of view of numerical implementation, so one can use it to solve the Hastings-Pwoell system of fractional-order stochastic differential equations. Grunwald-Letnikov (GLDα) fractional derivative of order α defined by [4043]

GLDαf(t)=Limh0hαj=0[tah](1)j(αj)f(tjh), (14)

where [tah] means the integer part of tah. This formula can be reduced to

GLDαf(tn)hαj=0nwjαf(tnj), (15)

where h is the time step, tn=nh and wjα are the Grunwald-Letnikov coefficients satisfy the following recurrence relationship

w0α=1,wjα=(11+αj)wj1α,j=1,2,3,

If f(t) is continuous function and f(t) is integrable function in the interval [0,T], then the relation between Caputo and Grunwald-Letnikov fractional derivative takes the form [4446]

CDαf(t)=GLDαf(t)f(0)tαΓ(1α)=GLDαf(t)Limn(1)nhα(α1n)f(0)1hαj=0nwjα(f(tnj)f(0)). (16)

Now, the fractional-order stochastic Hastings-Pwoell model (11) in Grunwald-Letnikov sense can be written as

xn=x0+hα(xn1(γxn1)a1xn1yn11+b1xn1+σ1xn1hζ1n)j=1nwjα(xnjx0)yn=y0+hα(a1xn1yn11+b1xn1a2yn1zn11+b2yn1μ1yn1ηyn12+σ2yn1hζ2n)j=1nwjα(ynjy0)zn=z0+hα(a2yn1zn11+b2yn1μ2zn1+σ3zn1hζ3n)j=1nwjα(znjz0), (17)

where, σi and ζi(n) represent real constants and a 3D Gaussian white noise processes, respectively, i=1,2,3.  ζi satisfy the follows:

ζj(t)=0,(j=0,1,2,3)andζi(t1)ζj(t2)=δijδ(t1tj),

δij is Kronecker delta and δ(t1tj) is the Dirac delta function.

7  Numerical Simulations

In this part, the numerical simulations of the mathematical model (3) will be examined, and the focus will be on the effect of harvesting parameters. The numerical results will be compared with the theorems formulated in the previous sections. The interactions between prey, middle predator and the top predator will be simulated by the following parameters: a1=2,a2=0.2,b1=2,b2=1,h=0.1,μ1=0.1,μ2=0.02.

In order to show the effects of fractional derivative α on the dynamics of system (3) one can draw the bifurcation diagram considering fractional-order α as a bifurcation parameter. The supercritical Hopf bifurcation value centralizes at α=0.912047 as indicated in Fig. 1. When α<α, the coexistence equilibrium point E3(0.8155,0.1111,2.88787) is locally asymptotically stable as shown in Fig. 1 and coincide with Fig. 2 when α=0.9. For α>α, the system undergoes a limit cycle oscillations as exhibits in Fig. 1 and coincide with Fig. 2 when α=0.92.

Figure 1: Bifurcation diagram of the fractional-order system (3) with respect to α

For better understand the effect of the quadratic harvesting of the middle predator η around the positive equilibrium point E3, one can draw the bifurcation diagram with respect to η as a bifurcation parameter. It can be seen that the supercritical Hopf bifurcation value localized at η=0.229688 as shown in Fig. 3. It can also be observed that when η<0.229688 the fractional-order system (3) undergoes limit cycle behaviour as shown in Fig. 3 and coincide with Fig. 4 when η=0.15. For η>0.229688 the interior equilibrium point E3=(0.815539,0.111111,2.64157418) is locally asymptotically stable as indicated in Fig. 3 and coincide with Fig. 4 when η=0.4. It can also be observed that the conditions of local stability that are proven in Theorem 7 are verified because a1b1y3(1+b1x3)2=0.0642022<1 and a2b2z3(1+b2y3)2=0.027935<η.

Figure 2: Time series and phase diagram of the equilibrium point E3 of model (3) with different values of α

Figure 3: Bifurcation diagram of the fractional-order system (3) with respect to η

The effect of prey harvesting rate h is shown in Fig. 5. From the bifurcation diagram with respect to h as a bifurcation parameter, it can be seen that the supercritical Hopf bifurcation value localized at h=0.3583 as shown in Fig. 5. For h<0.3583 the fractional-order system (3) undergoes limit cycle behaviour as indicated in Fig. 6 when h=0.2. It can also be observed that when h>0.3583 the interior equilibrium point is locally asymptotically stable as shown in Fig. 5 and coincide with Fig. 6 when h=0.4.

Figure 4: Time series and phase diagram of the equilibrium point E3 of model (3) with different values of η

Figure 5: Bifurcation diagram of fractional-order system (3) with respect to h

Figure 6: Time series and phase diagram of the equilibrium point E3 of model (3) with h = 0.02 and h = 0.4

Figure 7: Bifurcation diagram of fractional-order system (3) with respect to μ1

If we increase the value of middle predator death rate μ1 and keeping all other parameters value fixes, it can be seen that a transcritical bifurcation occurs at μ1=0.642857 as shown in Fig. 7 and stated in Theorem 9. For μ1<0.642857 the interior equilibrium point E3=(0.911269,0.111111,1.36443) is locally asymptotically stable as indicated in Fig. 7 and coincide with Fig. 8 when μ1=0.4. For μ1>0.642857 the free predator equilibrium point E3=(0.99,0,0) is locally asymptotically stable as shown in Fig. 7 and coincide with Fig. 8 when μ1=0.7.

Figure 8: Time series of the food chain model (3) with μ1 = 0.4 and μ1 = 0.7

8  Conclusion

In this paper, a deterministic and stochastic fractional-order model of the tri-trophic food chain model incorporating harvesting has been proposed. It is shown that the proposed model has bounded and non-negative solution as desired in any population dynamics. By using stability analysis of fractional-order system, we have proved that if the basic reproduction number R0<1, the predator free equilibrium point E1 is globally asymptotically stable. The interaction between prey, the middle predator, and the top predator was investigated. The impact of harvesting on prey and middle predator was studied. Some sufficient conditions that ensure the local and global stability of equilibrium points have been obtained. We conclude that harvesting parameters can control the dynamics of the middle predator. A numerical approximation method is developed for the proposed stochastic fractional-order model.

Funding Statement: The authors gratefully acknowledge Qassim University, represented by the Deanship of Scientific Research, on the financial support under the number (cosao-bs-2019-2-2-I-5469) during the academic year 1440 AH/2019 AD.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

## References

1. A. Hastings and T. Powell, “Chaos in a three-species food chain,” Ecology, vol. 72, no. 3, pp. 896–903, 199
2. B. Buonomo, F. Giannino, S. Saussure and E. Venturino, “Effects of limited volatiles release by plants in tritrophic interactions,” Mathematical Biosciences and Engineering, vol. 16, no. 5, pp. 3331–3344, 2019.
3. P. Ghosh, P. Das and D. Mukherjee, “Chaos to order—Effect of random predation in a Holling type IV tri-trophic food chain system with closure terms,” International Journal of Biomathematics, vol. 9, no. 5, pp. 1650073, 2016.
4. K. Cheng, H. You and T. Yang, “Dynamics of three species food chain model with Holling type II functional response,” arXiv preprint arXiv: 20012237v1, 2020.
5. B. Sahoo and S. Poria, “The chaos and control of a food chain model supplying additional food to top-predator,” Chaos, Solitons & Fractals, vol. 58, pp. 52–64, 2014.
6. S. G. Mortoja, P. Panja, A. Paul, S. Bhattacharya and S. K. Mondal, “Is the intermediate predator a key regulator of a tri-trophic food chain model?: An illustration through a new functional response,” Chaos, Solitons & Fractals, vol. 132, no. 4, pp. 109613, 2020.
7. N. Pal, S. Samanta and J. Chattopadhyay, “Revisited Hastings and Powell model with omnivory and predator switching,” Chaos, Solitons & Fractals, vol. 66, no. 12, pp. 58–73, 2014.
8. N. Pal, S. Samanta and S. Rana, “The impact of constant immigration on a tri-trophic food chain model,” International Journal of Applied and Computational Mathematics, vol. 3, no. 4, pp. 3615–3644, 2017.
9. S. Raw, B. Tiwari and P. Mishra, “Analysis of a plankton-fish model with external toxicity and nonlinear harvesting,” Ricerche di Matematica, vol. 69, no. 2, pp. 653–681, 2020.
10. M. Haque, N. Ali and S. Chakravarty, “Study of a tri-trophic prey-dependent food chain model of interacting populations,” Mathematical Biosciences, vol. 246, no. 1, pp. 55–71, 2013.
11. B. Ghosh, D. Pal, T. Legovi´c and T. Kar, “Harvesting induced stability and instability in a tri-trophic food chain,” Mathematical Biosciences, vol. 304, no. 5, pp. 89–99, 2018.
12.  B. Nath, N. Kumari, V. Kumar and K. P. Das, “Refugia and Allee effect in prey species stabilize chaos in a tri-trophic food chain model,” Differential Equations and Dynamical Systems, vol. 3, no. 12, pp. 591, 2019.
13. A. Matouk, A. Elsadany, E. Ahmed and H. Agiza, “Dynamical behavior of fractional-order Hastings-Powell food chain model and its discretization,” Communications in Nonlinear Science and Numerical Simulation, vol. 27, no. 1–3, pp. 153–167, 2015.
14. M. N. Huda, T. Trisilowati and A. Suryanto, “Dynamical analysis of fractional-order Hastings-Powell food chain model with alternative food,” Journal of Experimental Life Science, vol. 7, no. 1, pp. 39–44, 2017.
15. P. Panja, “Stability and dynamics of a fractional-order three-species predator-prey model,” Theory in Biosciences, vol. 138, no. 2, pp. 251–259, 2019.
16. L. Wang, H. Chang and Y. Li, “Dynamics analysis and chaotic control of a fractional-order three- species food-chain system,” Mathematics, vol. 8, no. 3, pp. 409, 2020.
17. L. Addison, “Analysis of a predator-prey model: A deterministic and stochastic approach,” J. Biom. Biostat., vol. 8, no. 4, pp. 359–368, 20
18. R. Ali, M. S. Arif, M. Rafiq, M. Bibi and R. Fayyaz, “Numerical analysis of stochastic vector borne plant disease model,” Computers, Materials & Continua, vol. 63, no. 1, pp. 65–83, 2020.
19. B. Ahmad, S. K. Ntouyas, A. Alsaedi and A. F. Albideewi, “A study of a coupled system of Hadamard fractional differential equations with nonlocal coupled initial-multipoint conditions,” Advances in Difference Equations, vol. 2021, no. 1, pp. 1–16, 2021.
20. S. Ahmad, R. Ullah and D. Baleanu, “Mathematical analysis of tuberculosis control model using nonsingular kernel type Caputo derivative,” Advances in Difference Equations, vol. 2021, no. 1, pp. 1–18, 2021.
21. R. Almeida, A. M. B. da Cruz, N. Martins and M. T. T. Monteiro, “An epidemiological MSEIR model described by the Caputo fractional derivative,” International Journal of Dynamics and Control, vol. 7, no. 2, pp. 776–784, 2019.
22. N. Sweilam, S. Al-Mekhlafi, A. Albalawi and D. Baleanu, “On the optimal control of coronavirus (2019-nCov) mathematical model; a numerical approach,” Advances in Difference Equations, vol. 2020, no. 1, pp. 1–13, 2020.
23. Q. Liu, D. Jiang, T. Hayat, A. Alsaedi and B. Ahmad, “Dynamical behavior of a higher order stochastically perturbed SIRI epidemic model with relapse and media coverage,” Chaos, Solitons & Fractals, vol. 139, no. Suppl. 3, pp. 110013, 2020.
24. A. I. Aliyu, M. Inc, A. Yusuf and D. Baleanu, “A fractional model of vertical transmission and cure of vector-borne diseases pertaining to the Atangana-Baleanu fractional derivatives,” Chaos, Solitons & Fractals, vol. 116, pp. 268–277, 2018.
25. I. Tejado, E. P´erez and D. Val´erio, “Fractional derivatives for economic growth modelling of the group of twenty: Application to prediction,” Mathematics, vol. 8, no. 1, pp. 50, 2020.
26. K. Shah, F. Jarad and T. Abdeljawad, “On a nonlinear fractional order model of dengue fever disease under Caputo-Fabrizio derivative,” Alexandria Engineering Journal, vol. 59, no. 4, pp. 2305–2313, 2020.
27. B. Nath and K. P. Das, “Harvesting in tri-trophic food chain stabilises the chaotic dynamics-conclusion drawn from Hastings and Powell model,” International Journal of Dynamical Systems and Differential Equations, vol. 10, no. 2, pp. 95–115, 2020.
28. J. Cresson and A. Szafran´ska, “Discrete and continuous fractional persistence problems-the positivity property and applications,” Communications in Nonlinear Science and Numerical Simulation, vol. 44, no. 3, pp. 424–448, 2017.
29. H. Li, L. Zhang, C. Hu, Y. Jiang and Z. Teng, “Dynamical analysis of a fractional-order predator- prey model incorporating a prey refuge,” Journal of Applied Mathematics and Computing, vol. 54, no. 1–2, pp. 435–449, 2017.
30. M. Sambath, P. Ramesh and K. Balachandran, “Asymptotic behavior of the fractional order three species prey-predator model,” International Journal of Nonlinear Sciences and Numerical Simulation, vol. 19, no. 7–8, pp. 721–733, 2018.
31. S. K. Choi, B. Kang and N. Koo, “Stability for Caputo fractional differential systems,” Abstract and Applied Analysis, vol. 2014, no. 3–4, pp. 1–6, 2014.
32. K. pada Das, “A mathematical study of a predator-prey dynamics with disease in predator,” International Scholarly Research Network, vol. 807486, pp. 1–16, 2011.
33. J. Huo, H. Zhao and L. Zhu, “The effect of vaccines on backward bifurcation in a fractional order HIV model,” Nonlinear Analysis: Real World Applications, vol. 26, no. 3, pp. 289–305, 2015.
34. M. S. Abdelouahab, N. E. Hamri and J. Wang, “Hopf bifurcation and chaos in fractional-order modified hybrid optical system,” Nonlinear Dynamics, vol. 69, no. 1–2, pp. 275–284, 2012.
35. L. Perko, Differential Equations and Dynamical Systems, vol. 7. New York, NY, USA: Springer, 2013.
36. W. Wang, S. Cheng, Z. Guo and X. Yan, “A note on the continuity for Caputo fractional stochastic differential equations,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 30, no. 7, pp. 73106, 2020.
37. T. Doan, P. Huong, P. Kloeden and A. Vu, “Euler-maruyama scheme for Caputo stochastic fractional differential equations,” Journal of Computational and Applied Mathematics, vol. 380, no. 4, pp. 112989, 2020.
38. Z. Guo, J. Hu and W. Wang, “Caratheodory’s approximation for a type of Caputo fractional stochastic differential equations,” Advances in Difference Equations, vol. 2021, no. 1, pp. 68, 2021.
39. D. T. Son, P. T. Huong, P. E. Kloeden and H. T. Tuan, “Asymptotic separation between solutions of Caputo fractional stochastic differential equations,” Stochastic Analysis and Applications, vol. 36, no. 4, pp. 654–664, 2018.
40. H. Aminikhah, A. H. R. Sheikhani, T. Houlari and H. Rezazadeh, “Numerical solution of the distributed-order fractional Bagley-Torvik equation,” IEEE/CAA Journal of Automatica Sinica, vol. 6, no. 3, pp. 760–765, 2019.
41. S. S. Ray and A. Patra, “Numerical solution of fractional stochastic neutron point kinetic equation for nuclear reactor dynamics,” Annals of Nuclear Energy, vol. 54, no. 4, pp. 154–161, 2013.
42. M. Hou, J. Yang, S. Shi and H. Liu, “Logical stochastic resonance in a nonlinear fractional-order system,” European Physical Journal Plus, vol. 135, no. 9, pp. L453, 2020.
43. G. Zou and B. Wang, “On the study of stochastic fractional-order differential equation systems,” arXiv preprint arXiv: 1611.07618, 2016.
44. I. Podlubny, Fractional Differential Equations. New York, NY, USA: Academic Press, 1999.
45. C. Huang, X. Yu, C. Wang, Z. Li and N. An, “A numerical method based on fully discrete direct discontinuous Galerkin method for the time fractional diffusion equation,” Applied Mathematics and Computation, vol. 264, no. 1, pp. 483–492, 2015.
46. C. Huang and M. Stynes, “Error analysis of a finite element method with GMMP temporal discretisation for a time-fractional diffusion equation,” Computers & Mathematics with Applications, vol. 79, no. 9, pp. 2784–2794, 2020.