# Dynamical Behaviors of a Fractional-Order Three Dimensional Prey-Predator Model

Authors: Noor S. Sh. Barhoom; Sadiq Al-Nassir
Journal: Abstract and Applied Analysis
 (2021)

---

## Abstract

In this paper, the dynamical behavior of a three-dimensional fractional-order prey-predator model is investigated with Holling type III functional response and constant rate harvesting. It is assumed that the middle predator species consumes only the prey species, and the top predator species consumes only the middle predator species. We also prove the boundedness, the non-negativity, the uniqueness, and the existence of the solutions of the proposed model. Then, all possible equilibria are determined, and the dynamical behaviors of the proposed model around the equilibrium points are investigated. Finally, numerical simulations results are presented to confirm the theoretical results and to give a better understanding of the dynamics of our proposed model.


---

## Body

## 1. Introduction

The most influential theme in ecology and mathematical modeling is the dynamic of the relationship among species. Many authors extended or modified the work of the Lotka and Voltera model [1, 2], and they have investigated thse topics widely by using ordinary differential equations or deference equations, see [3–9], and references therein. Nowadays, authors formulate their mathematical models by fractional order differential equation due to their ability to give the precise description for various linear and nonlinear phenomena [10–15]. The increasing of mathematical models that based on fractional order differential equation has recently obtained popularity in the studying the behavior of biological models. Fractional-order differential equation has been successfully used and applied to model many areas of science, engineering, and phenomena that cannot be formulated by other types of equations [10, 16]. In [12, 16–21], authors have investigated the effects of the fractional order differential equation on a prey predator model as well as they also discussed the stability analysis of equilibrium points of fractional order model with and without harvesting, as well as the existence, uniqueness, and boundedness of the solutions that are proved.There are several different types of definitions of fractional-order differential equation in the literatures [16, 22–25], for e.g., Caputo, Riemann-Litouville, Fabrizio, Marechand,Grunwald-Letnikov, wayl, and Riesz fractional-order derivatives. Throughout this work, we used the Caputo fractional-order derivatives because it is not necessary to define the initial conditions of fractional-order and its fractional-order derivative of constant function is zero, as well as the similarity of the initial conditions of fractional order differential and the integer order ones [8]. This work is organized as follows: in Section 2, some useful definitions and concept that concern to the fractional order are presented. In Section 3, a three-dimensional fractional-order prey-predator model is considered. The uniqueness and boundedness as well as the nonnegativity of its solution are proved. In Section 4, all the equilibrium points are determined, and the conditions to achieve its local stability are set. In Section 5, numerical simulations are given to confirm the theoretical results. Finally, conclusions are given in Section 6.


## 2. Preliminaries

Definition 1.
[24, 26]: Caputo’s definition of fractional derivatives is given as follows:
Dtαft=Iα−αfαt, α>0,
whereαis the least integer which is not less than α, and Iθ is the Riemann-Liouville integral operator of order θ which is given by
Iθft=1/Γθ∫0tt−τθ−1fτdτ, θ>0.
Γθ is the Euler Gamma function.Throughout this work, we need some useful results for the fractional-order derivative. We start with the following results which are proved in [13, 16, 27].Lemma 1.(1)
Suppose thatft∈Ca,b and Daαft∈Ca,b with 0<α≤1, then we have(1)ft=fa+1ΓαDaαfξt−aα,wherea≤ξ≤x, ∀x∈a,b.
(2)
Supposeft∈Ca,b and Dtαft∈Ca,b, 0<α≤1. IfDaαft≥0, ∀t∈a,b, then ft is a nondecreasing function ∀t∈a,b, and if Dt0caαft≤0, ∀t∈a,b, then ft is a nonincreasing function ∀t∈a,b(3)
Consider the Cauchy problem(2)Daαxt=λxt+ft,xa=bb∈ℝ.With0<α≤1 and λ∈ℝ, then the form of the solution is given by
(3)xt=bEαλt−aα+∫att−sα−1Eα,αλt−sαfsds.While the solution to the problem(4)Daαxt=λxt,xa=bb∈ℝ,is given by
(5)xt=bEαλt−aα.(4)
Letut be a continuous function on t0,∞ and satisfy(6)Daαut≤−λut+μ,where0<α<1,λ,μ∈ℝ2,λ≠0, and t0≥0 is the initial time. Then, its solution has the form
(7)ut≤ut0−μλEα−λt−t0α+μλ.For the existence a unique solution for the general fractional-order derivative, we have the next lemma that appeared in [28].Lemma 2.
Consider the system
Dtαxt=ft,x, t>t0
with the initial conditionxt0, where0<α≤1, f:t0,∞×Ω⟶ℝn. If ft,x satisfies the locally Lipchitz condition with respect to x, then there exists a unique solution of the above system on t0,∞×Ω.


## 3. Model Formulation

Frist, the following three dimension prey-predator model with functional response of Holling type III is considered as follows:(8)dXdT=R0X1−XK0−a1YX2B1+X2dYdT=a2YX2B1+X2−b1ZY2B2+Y2−c1Y−h2YdZdT=b2ZY2B2+Y2−c2Z,whereX, Y,andZ are the density of prey species, middle species, and the top predator species. R0, K0,,c1,c2, and b1are the intrinsic growth rates, the carrying capacity of the prey population, the death rates for middle predator population and top predator population, and capture rate of prey and middle predator, respectively. The parametera1is the conservation rate of preyX, and the parameters a2and b2 are the conservation rate of middle predator Y to the top predator Z, respectively. The constantsB1 and B2 are the predation rates and half saturation constants, respectively. The parameter h2 denotes to the constant harvesting rate.We use the following transformation to reduce the dimension of model (8). x=X/k,y=a1/k0R0Y,z=a1b1/R02k0Z, and t=R0T, and then we obtain the following system:
(9)dxdt=x1−x−ax2y1+ax2dydt=bx2y1+ax2−dy2z1+dy2−m1y−hydzdt=cy2z1+dy2−m2z,wherea=k02/B1, b=aa2/R0, c=R0d/b2, d=k02R02/a12B2, m1=c1k02R0/a1, m2=R0/c2, and h=h2k0R0/a1.Next, we introduce the fractional–order derivativeα in model (2) with the help of fractional order Caputo type derivative. Then, the system (9) becomes as follows:
(10)Dtαxt=x1−x−ax2y1+ax2Dtαyt=bx2y1+ax2−dy2z1+dy2−m1y−hyDtαzt=cy2z1+dy2−m2z.Remark 3.
Because of the biological significance, we are interesting solutions that are nonnegative and bounded only, so that we will prove the following theorem.Theorem 4.
LetΩ+=x,y,z∣x≥0,y≥0andz≥0 and denotes all nonnegative real numbers in ℝ3, and then all solutions of the system (10) with x0≥0,y0≥0 and z0≥0 are uniformly bounded and nonnegative.Proof.
we have to show thatxt≥0,∀t≥0, assuming x0>0 for t=0 considering that xt≥0,∀t≥0 is not true. Then, there exists a constant t1>0 such that
(11)xt>0,0≤t<t1,xt=0,t=t1,xt<0,t>t1.
From the first equation of system (10), we have
(12)Dt1αxtt=t1=0.
According to part 1 in Lemma 1, we havext1+=0, which contradicts the fact that xt1+<0, i.e., xt<0,t>t1. Therefore, we have xt≥0,∀t≥0. By the same arguments, yt≥0,∀t≥0 and zt≥0 can get ∀t≥0.For the uniformly bounded, Let Vt=1/ax+1/by+d/bcz, and then we have
(13)DαVt=1aDαxt+1bDαyt+dbcDαzt=1ax1−x−1bm1+hy−dbcm2z.
Now, for eachη>0, we have
(14)DαVt+ηVt=1ax1−x−1bm1+hy−dbcm2z+η1ax+η1by+ηdbcz≤−1ax2+η+11ax+η−m1+h1by+η−m2dbcz≤−1ax2−η+1x+14η+12−14η+12+η−m1+h1by+η−m2dbcz=−1ax−12η+12+η+124a+η−m1+h1by+η−m2dbcz≤14aη+12+η−m1+h1by+η−m2dbcz.
Therefore,
DαVt+ηVt≤l, by taking η<minm1+h,m2,
wherel=1/4aη+12>0.
Now by using part 4 in Lemma 1, we have(15)Vt≤V0−lηEα−ηt−0α+lη,≤V0Eα−ηtα+lη1−Eα−ηtα.
Thus,Vt⟶l/η as t⟶∞ and 0<Vt≤l/η. Therefore, all solutions of system (10) that starts from Ω+ are confined in the regionΩ=x,y,z∈ℝ3∣Vt≤l/η+ε,foranyε>0.Now, the existence and uniqueness of the solution of the considered system (10) are shown in the next theorem:Theorem 5.
Letγ be a sufficiently large, and then for each S0=x0,y0,z0∈x,y,z∈ℝ3∣maxx,y,z≤γ, there exists a unique solution S=x,y,z∈Ω of the fractional-order system (10) with initial condition S0, which is defined for all t≥0.Proof.
LetS=x,y,z, S¯=x¯,y¯,z¯, and consider a mapping HS=H1S,H2S,H3S, such that
(16)H1S=x1−x−ax2y1+ax2H2S=bx2y1+ax2−dy2z1+dy2−m1y−hyH3S=cy2z1+dy2−m2z.
For anyS,S¯∈Ω, we have
(17)HS−HS¯=H1S−H1S¯+H2S−H2S¯+H3S−H3S¯=x1−x−ax2y1+ax2−x¯1−x¯+ax¯2y¯1+ax¯2+bx2y1+ax2−dy2z1+dy2−m1y−hy−bx¯2y¯1+ax¯2+dy¯2z¯1+dy¯2+m1y¯+hy¯+cy¯2z¯1+dy¯2−m2z−cy¯2z¯1+dy¯2+m2z¯=x−x¯1−x+x¯−ax2y−x¯2y¯+ax2x¯2y−y¯1+ax21+ax¯2+bx2y−x¯2y¯+ax2x¯2y−y¯1+ax21+ax¯2−dy2z−y¯2z¯+dy2y¯2z−z¯1+dy21+dy¯2−m1y−y¯−hy−y¯+cy2z−y¯2z¯+dy2y¯2z−z¯1+dy21+dy¯2−m2z−z¯.
Since1/1+dy21+dy¯2≤1 and 1/1+ax21+ax¯2≤1 and maxx,y,z≤γ, we can have
(18)HS−HS¯≤x−x¯1−x+x¯+a+bx2y−x¯2y¯+ax2x¯2y−y¯1+ax21+ax¯2+d+cy2z−y¯2z¯+dy2y¯2z−z¯1+dy21+dy¯2+m1+hy−y¯+m2z−z¯≤1+2γ+2γ2a+bx−x¯+a+bγ2+aa+bγ4+2γ2d+c+m1+hy−y¯+d+cγ2+dd+cγ4+m2z−z¯=LS−S¯,
whereL=max1+2γ+2γ2a+b,a+bγ2+aa+bγ4+2γ2d+c+m1+h,d+cγ2+dd+cγ4+m2.Therefore,HS satisfies the Lipchitz condition. By Lemma 2, the proof is finished.


## 4. The Equilibrium Points and Local Stability Analysis

In this section, the existence and the local stability of equilibrium points of the considered system (10) are studied and investigated. Before that, we need the following theorem.Theorem 6 (see [29]).
Consider the fractional–order differential systemdαx/dtα=fx; x0=x0, with α∈0,1 and x∈ℝn. The equilibrium points of the above system are solutions to the equation fx=0. An equilibrium point is locally asymptotically stable if all the eigenvalues λj of the Jacobian matrix J=∂f/∂x evaluated at the equilibrium satisfy argλj>απ/2.Proof (see [29]).
To obtain the equilibrium points of the considered system (10), we solve the following simultaneous equations:
(19)Dtαxt=0,Dtαyt=0,Dtαzt=0.
Thus, the considered system (10) has the following possible equilibrium points:
(1)
The trivial equilibrium pointE0=0,0,0 and the equilibrium point E1=1,0,0 without any restriction on the parameters(2)
The top predator free equilibrium pointE2=x¯,y¯,0 exists if b>am1+h and 0<x¯<1, where x¯=m1+h/b−am1+h and y¯=1/ax1−x¯1+ax¯2(3)
The interior equilibrium pointE3=x∗,y∗,z∗ exists if c>m3d and bx∗2/1+ax∗2>m1+h, where y∗=m2/c−m2d, x∗, is the positive root of the equation ax3−ax2+ay∗+1x−1=0 and z∗=1+dy∗2/dy∗bx∗2/1+ax∗2−m1+h
The Jacobian matrix of system (10) at x,y,z is given by
(20)Jx,y,z=1−2x−2axyk2−ax2k02bxyk2bx2k−2dyzk12−m1−h−dy2k102cyzk12cy2k1−m2,
wherek=1+ax2 and k1=1+dy2.The next theorems give the local stability ofE0, E1, and E2, respectively.Theorem 7.
The trivial equilibrium pointE0 of the system (10) is always unstable point.Proof.
The Jacobian matrixJ at E0 is given by
(21)JE0=1000−m1+h000−m2.The characteristic equation of the Jacobian matrix JE0 is given by
(22)1−λ−m1+h−λ−m2−λ=0.
Therefore, the roots of the above characteristic equation areλ1=1, λ2=−m1+h, and λ3=−m3. Then, argλ1=0<απ/2, argλ2=π>απ/2 and also argλ3=π>απ/2. Hence, the trivial equilibrium point is always unstable point.Theorem 8.
The equilibrium pointE1 of system (10) is locally asymptotically stable if <m1+h1+a.Proof.
The Jacobian matrixJ at E1 is given by
(23)JE1=−1−a1+a00b1+a−m1+h000−m2.The characteristic equation of the Jacobian matrix JE1 is given by
(24)−1−λb1+a−m1+h−λ−m2−λ=0.
The roots of the above characteristic equation areλ1=−1, λ2=b/1+a−m1+h and λ3=−m2. Then, argλ1=argλ2=π>απ/2 and argλ3=π>απ/2. Hence, the equilibrium point E1 is locally asymptotically stable if <m1+h1+a . Otherwise, it is unstable point.Theorem 9.
The equilibrium pointE2 is locally asymptotically stable if one of the following holds:
(i)
cy2/k1<m2l1<0, l2>0, and ∆≥0(ii)
cy2/k1<m2l1<0 and ∆<0 and argλi=tan−1−∆/l1>απ/2,i=2,3 where λ2,3=l1∓l12−4l2/2, l1=1−2x−21−x−bx2/k−m1+h, l2=bx2/k1−2x−m1+h1−2x−21−x/k, and ∆=l12−4l2Proof.
The Jacobian matrixJ at E2 is given by
(25)JE2=1−2x−2axyk2−ax2k02bxyk2bx2k−m1+h−dy2k100cy2k1−m2.
The characteristic equation of the Jacobian matrixJE2 is given by
(26)cy2k1−m2−λλ2−l1λ+l2=0.
It is clear thatλ1=cy2/k1−m2 is negative if cy2/k1<m2 so that argλ1=π>απ/2. The other two roots λ2 and λ3 are λ2,3=l1∓l12−4l2/2. Hence, if ∆≥0. then the eigenvalue λ2 and the eigenvalue λ3 are taken negative value; thus, argλ2=argλ3=π>απ/2. Therefore, the point E2 is locally asymptotically stable. Also, if ∆<0 and if argλi=tan−1−∆/l1>0,i=2,3 hold, then the equilibrium point E2 is asymptotically stable. Otherwise, it is unstable point.In order to discuss the local stability of the interior point, we need the following definition and lemma.Definition 10 (see [30]).
The discriminantDf of a polynomial
(27)fx=xn+a1xn−1+a2xn−2+⋯+anis defined by Df=−1nn−1/2Rf,f′, where f′ is the derivative of f and where gx=xn+b1xl−1+b2xl−2+⋯+bl and Rf,g is an n+l⨂n+l determinant. If n=3, then Df is given by
(28)Df=18a1a2a3+a1a22−4a3a13−4a23−27a32.Lemma 11 (see [30]).
LetDp be a discriminant of the cubic polynomial equation, pλ=λ3+a1λ2+a2λ+a3=0, then
(1)
If Dp>0,a1>0,a3>0,a1a2>a3, then argλi>απ/2,i=1,2,3∀α∈0,1(2)
IfDp<0,a1≥0, a2≥0, and a3>0, then argλi>απ/2,i=1,2,3 forα<2/3(3)
IfDP<0, a1<0, a2<0, and α>2/3, then argλi<απ/2,i=1,2,3(4)
IfDP<0, a1>0, a2>0, and a1a2=a3, then argλi>απ/2,i=1,2,3. ∀α∈0,1Proof (see [30]).Theorem 12.
The interior equilibrium pointE3 is locally asymptotically stable if one of the following conditions holds:
(1)
IfDp>0,N1>0,N3>0 and N1N2−N3>0, for 0<α<1(2)
IfDp<0,N1≥0,N2≥0, N3>0, for 0<α<2/3(3)
IfDp<0,N1>0,N2>0, N1N2=N3, for 0<α<1
whereN1=−1+2x+2axy/k2−bx2/k+2dyz/k12+m1+h−cy2/k1+m2, N2=1−2x−2axy/k2bx2/k−2dyz/k12−m1+h+cy2/k1−m2+bx2/k−2dyz/k12−m1+hcy2/k1−m2+2abx3y/k3+2cdy3z/k13, and N3=−1−2x−2axy/k2bx2/k−2dyz/k12−m1+hcy2/k1−m2+2cdy3z/k13−2abx3y/k3cy2/k1−m2.Proof.
The Jacobian matrix of system (10) at E3 is given by
(29)JE3=1−2x−2axyk2−ax2k02bxyk2bx2k−2dyzk12−m1+h−dy2k102cyzk12cy2k1−m2.
Therefore, the characteristic equation ofJE3 is given by
(30)pλ=λ3+N1λ2+N2λ+N3=0.
Thus, the discriminantDp of the cubic polynomial pλ is
(31)Dp=−1N1N201N132N1N2N30N2N300032N1003N202N1N2,
whereDp=18N1N2N3+N1N22−4N13N3−4N23−27N32.Then, all results can be obtained by Lemma 6.


## 5. Numerical Simulation

All previous theoretical results in Section4 are confirmed by giving numerical simulations. At different set of values of parameters, the local stability of E1, E2, and E3 of the considered system (10) is investigated numerically.First, for the pointE1, the following set of values of the parameters is chosen :a1=0.5;c1=0.1;b1=0.5;d1=0.2;m1=0.15;m2=0.1;h=0.8,α=0.98, and (0.1,0.5,0.9) is the initial value. Therefore, the condition in Theorem 4 is satisfied. Figure 1 shows that the point E1 is locally stable.Figure 1
This figure shows the local stability of the pointE1 and the time series of the three species. (a) The point E1. (b) Time series of the three species.
(a)(b)For the pointE2, we choose the values of parameters as following E2, a1=0.6; c1=0.1; b1=0.3; d1=0.2; m1=0.03; m2=0.32; and h=0.01,α=0.9 5, and (0.5,2.5,03) is the initial value. According to Theorem 9, the local stability of E2 is satisfied. Figure 2 illustrates the local stability of the point E2.Figure 2
This figure shows the local stability of the pointE2 and the time series of the three species. (a) The point E2. (b) Time series of the three species.
(a)(b)To show the local stability of the unique positive equilibrium pointE3, these values of parameters are chosen a1=0.75; c1=0.17; b1=0.3; d1=0.3; m1=0.03; m2=0.3; h=0.02, and α=0.9 5, (0.82,1.95,0.9) is the initial value so that according to Theorem 12, the equilibrium point E3 is local stability point, and the trajectories of prey and predator are shown in Figure 3.Figure 3
This figure shows the local stability of the pointE3 and the time series of the three species. (a) The point E3. (b) Time series of the three species.
(a)(b)


## 6. Conclusions

Fractional-order differential equations have recently been successfully used and applied to model many areas of science, especially, phenomena that cannot be formulated by other types of equations, so that this work concerns a study of a three-dimensional fractional-order prey-predator model with Holling type III functional response and constant rate harvesting. It is shown that the model possesses the existence, nonnegativity, boundedness, and uniqueness of the solutions, as desired in any population dynamics. The local stability of all possible equilibrium points of the considered system is investigated analytically then confirms by numerical verifications.



---

*Source: 1366797-2021-12-15.xml*