Computers, Materials & Continua DOI:10.32604/cmc.2022.019314
Article

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

Moustafa El-Shahed1,* and Asmaa M. Al-Dububan2

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 [2–16]. 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 [19–26]. 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(1−Xk)−βXYA+X−H1X,dYdT=e1βXYA+X−γYZB+Y−dY−H2Y2,dZdT=e2γYZB+Y−δZ−H3Z, (1)

with initial values X(0)=X0≥0,Y(0)=Y0≥0,Z(0)=Z0≥0.

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

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

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

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

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+b1x−a2yz1+b2y−μ1y−ηy2,cDαz=a2yz1+b2y−μ2z, (3)

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

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

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 t≥0.

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

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

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

where

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 X∈R+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

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

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

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

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

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η(R1−1), and x2 satisfies the equation

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

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

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

E3(x3,y3,z3) exists if

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)=(γ−2x−a1y(1+b1x)2−a1x1+b1x0a1y(1+b1x)2a1x1+b1x−a2z(1+b2y)2−μ1−2ηy−a2y1+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

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

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

Choosing L1=1−R0γR0, one obtains cDαV1≤0.

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

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−ηy2−a2y21+b2y200a2y21+b2y2−μ2), (8)

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

The above equation has the following eigenvalue λ1=a2y21+b2y2−μ2=μ2(R2−1). 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.

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

Choosing L2=11+b1x2, then cDαV2≤0 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

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

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

where

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

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

Choosing L3=11+b1x3 and L4=1+b2y3, then cDαV3≤0 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

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,

then W2TDFμ1(E1,μ1∗)W1=ν1τ2(b1γ+1)a1≠0, and

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+b1x−a2yz1+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)(t−s)α−1Γ(α)ds+∫0tg(X)(t−s)α−1dσ(s)Γ(α)ds. (13)

According to [36–39], 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 [40–43]

GLDαf(t)=Limh→0h−α∑j=0[t−ah](−1)j(αj)f(t−jh), (14)

where [t−ah] means the integer part of t−ah. This formula can be reduced to

GLDαf(tn)≈h−α∑j=0nwjαf(tn−j), (15)

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

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 [44–46]

CDαf(t)=GLDαf(t)−f(0)t−αΓ(1−α)=GLDαf(t)−Limn→∞(−1)nhα(α−1n)f(0)≈1hα∑j=0nwjα(f(tn−j)−f(0)). (16)

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

xn=x0+hα(xn−1(γ−xn−1)−a1xn−1yn−11+b1xn−1+σ1xn−1hζ1n)−∑j=1nwjα(xn−j−x0)yn=y0+hα(a1xn−1yn−11+b1xn−1−a2yn−1zn−11+b2yn−1−μ1yn−1−ηyn−12+σ2yn−1hζ2n)−∑j=1nwjα(yn−j−y0)zn=z0+hα(a2yn−1zn−11+b2yn−1−μ2zn−1+σ3zn−1hζ3n)−∑j=1nwjα(zn−j−z0), (17)