# Fractional-Order Model of Two-Prey One-Predator System

Authors: Mohammed Fathy Elettreby; Ahlam Abdullah Al-Raezah; Tamer Nabil
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

We propose a fractional-order model of the interaction within two-prey and one-predator system. We prove the existence and the uniqueness of the solutions of this model. We investigate in detail the local asymptotic stability of the equilibrium solutions of this model. Also, we illustrate the analytical results by some numerical simulations. Finally, we give an example of an equilibrium solution that is centre for the integer order system, while it is locally asymptotically stable for its fractional-order counterpart.


---

## Body

## 1. Introduction

The branch of mathematics that concerned biology is called mathematical biology. Mathematical biology tries to model, study, analyze, and interpret biological phenomenon such as the interactions, coexistence, and evolution of different species [1–4]. These interactions may be among the individuals of the same species and among the individuals of different species or interactions against the environment, disease, and food supply. The most important one is the interaction among the individuals of different species which can be predatory, competitive, or mutualistic. A large number of simple mathematical models have been suggested to understand the predator-prey interaction. The work done independently by Lotka [5] and Volterra [6], which is known as Lotka-Volterra model, was the first stone in this field. Later, the model was extended to include density-dependent prey growth and a functional response [7]. After that, a huge number of variants of this model were suggested. Initially, the authors proposed a system of two-dimensional coupled differential equations [8, 9]. Then, the scenario changed to discrete mathematical models which included a lot of complex dynamical behaviours [10, 11]. In the last two decades, the fractional-order differential equations appeared and began to study the predator-prey models in the fractional-order form [12–19].Fractional calculus generalizes the concept of ordinary differentiation and integration to noninteger order. Fractional calculus is a fertile field for researchers to study very important real phenomena in many fields like physics, engineering, biology, and so forth [20–30]. The fractional differential equations are naturally related to systems with memory. Also, they are closely related to fractals which are numerous in biological system. The definition of fractional derivative involves an integration which is nonlocal operator. This means that the fractional derivative is a nonlocal operator. Studies and results of solutions obtained using the fractional differential equations are more general and are as stable as their integer order counterpart.There are a lot of approaches to define the fractional differential operator such as Grunwald-Letnikov, Riemann-Liouville, Caputo, and Hadamard. The Riemann-Liouville and Caputo approaches are the most widely used in applications [21, 22, 29, 31]. The fractional derivative is defined via the fractional integral operator. So, we will start by the definition of the fractional integral operator. The fractional integral operator of order α (>0) of a function f(t) such that t>0 is given by(1)Iαft=∫0tt-sα-1Γαfsds.The Riemann-Liouville derivative of order α (>0) is given by(2)Dαft=DmIm-αft,where m=α and D=d/dt.In this paper, we used Caputo approach to define the fractional derivative. It is a modification to the Riemann-Liouville definition. The Caputo fractional derivative of orderα (>0) is denoted by D∗α and is given in the following form:(3)D∗αft=In-αDnft,where n=α and t>0.The following are some important properties of the fractional derivatives and integrals [18]. Let β,γ∈R+ and α∈0,1; then we have the following:(1)
IfIαa:L1→l1 and fx∈L1, then IaβIaγfx=Iaβ+γfx.(2)
limβ→n⁡Iaβfx=Ianfx uniformly, where n=1,2,3,… and Ia1fx=∫axds.(3)
limβ→0⁡Iaβfx=fx weakly.(4)
Ifk is constant and fx=k≠0, then D∗αk=0.(5)
Iffx is absolutely continuous on a,b, then limα→1⁡D∗αfx=dfx/dx.In this paper, we proposed a fractional-order model to study the interaction of a system consists of two-prey and one-predator species. In Section2, we proved the existence and the uniqueness of the solutions of our model. In Section 3, we studied the local asymptotic stability of the equilibrium points of the system. The numerical solution of the fractional-order two-prey one-predator model is given in Section 4.


## 2. The Fractional-Order Prey-Predator Model

Letx1t represent the first prey herds (gazelles) and x2t be the second prey herds (buffalos) densities, respectively. Suppose x3t represent the predator (lions) density. It is known that the logistic scenario is the most appropriate to describe the growth of the preys. So, the terms ax1t1-x1t and bx2t1-x2t are the growth of the two preys, where the positive parameters a and b are the intrinsic growth rate of them. The nature requires cooperations between preys against the predator to avoid the predation and to facilitate getting food. This cooperation is represented by the term x1tx2tx3t. Also, the predation to the preys is represented by the terms x1tx3t and x2(t)x3(t). Considering these assumptions, we get the following fractional-order two-prey one-predator system:(4)D∗αx1t=f1x1,x2,x3=ax1t1-x1t-x1tx3t+x1tx2tx3t,t∈0,T,D∗αx2t=f2x1,x2,x3=bx2t1-x2t-x2tx3t+x1tx2tx3t,t∈0,T,D∗αx3t=f3x1,x2,x3=-cx32t+dx1tx3t+ex2tx3t,t∈0,T,with the initial values (5)x1tt=0=x10,x2tt=0=x20,x3tt=0=x30,where c is the death rate of the predator, 0<α≤1, and x1(t)≥0, x2(t)≥0, x3(t)≥0. The constants a,b,c,d, and e are all positive. In the following, we studied the above model to understand the long time behaviour prey-predator interaction.Lemma 1.
The initial value problem (1, 2) can be written as the following matrix form:(6)D∗αXt=AXt+x1tBXt+x2tCXt+x3tDXt+x1tx2tEXt,X0=X0,where(7)Xt=x1tx2tx3t,X0=x10x20x30,A=a000b0000,B=-a00000000,C=0000-b0000,D=-1000-10de-c,E=001001000.Definition 2.
LetC∗0,T be the class of continuous column vector Xt=x1tx2tx3t where C0,T is the class of continuous functions defined on the interval 0,T and xit∈C0,T, i=1,2,3.Lemma 3.
The initial value problem (1, 2) has a unique solution in the region0,T×η, where η=x1,x2,x3∈R3:max⁡x1,x2,x3≤M and T<+∞.Proof.
LetFXt=AXt+x1tBXt+x2tCXt+x3tDXt+x1tx2tEXt be a continuous function where Xt is a continuous column vector. Suppose that Xt and Yt are two distinct continuous column vectors solutions of the initial value problem (1, 2) such that Xt=x1tx2tx3t, Yt=y1ty2ty3t. Then, (8)FX-FY=AXt+x1tBXt+x2tCXt+x3tDXt+x1tx2tEXt-AYt+y1tBYt+y2tCYt+y3tDYt+y1ty2tEyt≤AXt-Yt+x1tBXt-Yt+x1t-y1tBYt+x2tCXt-Yt+x2t-y2tCYt+x3tDXt-Yt+x3t-y3tDYt+x1tx2tEXt-Yt+x1tx2t-y1ty2tEYt.Since x1tx2t-y1ty2t=x1tx2t-x1ty2t+x1ty2t-y1ty2t≤x1tx2t-y2t+y2tx1t-y1t for each i=1,2,3.
Then, we have(9)FX-FY≤A+Bx1t+Yt+Cx2t+Yt+Dx3t+Yt+Ex1tx2t+x1t+y2tXt-Yt.Then (10)FX-FY≤LXt-Yt,where (11)L=A+4MB+4MC+4MD+M2+2ME.So, the continuous function FXt satisfies Lipschitz condition and has a unique solution [32].


## 3. Equilibrium Solutions and Stability Analysis

So far there is no known method for solving nonlinear equations. Therefore, it is difficult and even impossible to find an analytical solution to system (4). So, we will need the qualitative study. Finding the equilibrium points and studying their stability are the most important. To evaluate the equilibrium points of system (4), let (12)D∗αx-1=f1x-1,x-2,x-3=0,D∗αx-2=f2x-1,x-2,x-3=0,D∗αx-3=f3x-1,x-2,x-3=0,where x-1,x-2,x-3 is the equilibrium point of system (4).Solving the equations of system (12), we get the following equilibrium points. The dormant state E00,0,0, the boundary states E11,0,0, E20,1,0, E31,1,0, E40,bc/bc+e,be/bc+e, E5ac/ac+d,0,ad/ac+d, the two coexistence states E61,1,d+e/c, and E7(bca+e(a-b)/db+ea,acb-d(a-b)/db+ea,ab) which exists under the following conditions:(13)eba≤bc+e,dab≤ac+d.To study the stability of these equilibrium points, we have to linearize system (4) and compute its Jacobian matrix:(14)J=a1-2x1-x31-x2x1x3-x11-x2x2x3b1-2x2-x31-x1-x21-x1dx3ex3-2cx3+dx1+ex2.Substituting by the equilibrium pointE00,0,0 in the above Jacobian matrix, we get(15)J0=a000b0000,which has the eigenvalues λ=0,a,b. Since we have two positive eigenvalues, then the equilibrium point E0 is unstable. This means that the state of extinction of all species is not acceptable. There will be no life. The Jacobian matrix for the E11,0,0 is (16)J1=-a0-10b000d,and their eigenvalues are λ=-a,d,b. Then the equilibrium point E1 is unstable. Similarly, the eigenvalues of the equilibrium point E20,1,0 are λ=-b,a,e. So, it is unstable. Similarly, the eigenvalues for the equilibrium point E31,1,0 are λ=-a,-b,d+e. So, point E3 is unstable. This means that the state of existence of one prey alone cannot be continuous forever.For equilibrium pointE40,bc/bc+e,be/bc+e, the Jacobian matrix is(17)J4=a-be2bc+e200b2cebc+e2-b2cbc+e-bcbc+ebdebc+ebe2bc+e-bcebc+e,which has the following characteristic polynomial: (18)a-be2bc+e2-λλ2+a1λ+a2=0,where a1=bc(b+e)/bc+e and a2=b2ce/bc+e.A sufficient condition to say that an equilibrium point is a locally asymptotically stable is that all eigenvaluesλ satisfy argλ>απ/2 [33]. This condition implies that the characteristic polynomial of that point should satisfy Routh-Hurwitz conditions [17]. Since a1 and a2 are both positive, then the stability of point E4 depends on the first bracket in (18). Thus E4 is locally asymptotically stable if λ=a-be2/bc+e2<0.Similarly, equilibrium pointE5(ac/ac+d,0,ad/ac+d) has the Jacobian matrix: (19)J5=-a2cac+da2cdac+d2-acac+d0b-ad2ac+e20ad2ac+dadeac+d-acdac+d,which has the characteristic polynomial: (20)b-ad2ac+d2-λλ2+b1λ+b2=0,where b1=ac(a+d)/ac+d and b2=a2cd/ac+d. Since b1 and b2 are both positive, then E5 is locally asymptotically stable if λ=b-ad2/ac+d2<0.The study of the coexistence points is more important. This study gives us the conditions that lead to coexistence between the prey and the predator. The Jacobian matrix for the first coexistence equilibrium pointE6(1,1,d+e/c) has the characteristic polynomial: (21)d+e+λλ2+c1λ+c2=0,where c1=a+b and c2=d+e2/c2+ab. Since c1 and c2 are both positive and λ=-d+e<0. Thus E6 is locally asymptotically stable.For the equilibrium pointE7, the Jacobian matrix is given by(22)J7=-ax-1x-1x-3-x-11-x-2x-2x-3-bx-2-x-21-x-1dx-3ex-3-cx-3,where x-1,x-2,x-3=bca+ea-b/db+ea,acb-da-b/db+ea,ab. The above matrix has the following characteristic polynomial: (23)Pλ=λ3+d1λ2+d2λ+d3=0,whered1=ax-1+bx-2+cx-3,d2=acx-1x-3+bcx-2x-3+ex-2x-31-x-1+dx-1x-31-x-2,d3=x-1x-2x-3ae1-x-1+bd1-x-2+cx-3+d+ex-3.Since 0<x-1<1, 0<x-2<1, and x-3>0, then, from the above, d1>0, d2>0, and d3>0.Definition 4.
The discriminantDPλ of the polynomial Pλ=λ3+d1λ2+d2λ+d3 is defined by DPλ=-13RP,P′, where P′ is the derivative of the function P and RP,P′ is the determinant of the corresponding Sylvester matrix; then,(24)DPλ=-1d1d2d3001d1d2d332d1d200032d1d200032d1d2,DPλ=18d1d2d3+d1d22-4d13d3-4d23-27d32.Proposition 5.
If the discriminantDPλ is positive, then the Routh-Hurwitz conditions (d1>0, d3>0 and d1d2>d3) are the necessary and sufficient conditions for that all eigenvalues satisfy argλ>απ/2.Proposition 6.
From the above definition of the discriminant and the coefficients of the characteristic equation, we find that (with long calculations) all the four necessary and sufficient conditions are satisfied. Then equilibrium pointE7 is locally asymptotically stable.


## 4. Numerical Results

The Adams-type predictor-corrector method for the numerical solution of the fractional differential equations was discussed in [10]. This method can be used for both linear and nonlinear problems. It may be extended to multiterm equations (involving more than one differential operator) too [10].In this paper, we used Adams-type predictor-corrector method for the numerical solution of our system. First, we will give the Adams-type predictor-corrector method for solving general initial value problem with Caputo derivative:(25)D∗αyt=ft,yt,with the initial condition y0=y0 and t∈0,T. We assumed that a set of points tj,yj, where yj=ytj, are the points used for our approximation and tj=jh, j=0,1,…,Ninteger, h=T/N. The general formula for Adams-type predictor-corrector method is(26)yn+1=∑k=0α-1tn+1kk!y0k+hαΓα+2∑j=0nσj,n+1ftj,yj+hαΓα+2σn+1,n+1ftn+1,yn+1P,where(27)σj,n+1=nα+1n-αn+1α,ifj=0n-j+2α+1+n-jα+1-2n-j+1α+1,if1≤j≤n1,ifj=n+1,yn+1P=∑k=0α-1tn+1kk!y0k+1Γα∑j=0nρj,n+1ftj,yj,ρj,n+1=hααn+1-jα-n-jα.Applying the above algorithm for system (4), we have the following:(28)x1,n+1=x1,0+hαΓα+2∑j=0nσ1,j,n+1ax1,j1-x1,j-x1,jx3,j+x1,jx2,jx3,j+hαΓα+2σ1,n+1,n+1ax1,n+1P1-x1,n+1P-x1,n+1Px3,n+1P+x1,n+1Px2,n+1Px3,n+1P,x2n+1=x2,0+hαΓα+2∑j=0nσ2,j,n+1bx2,j1-x2,j-x2,jx3,j+x1,jx2,jx3,j+hαΓα+2σ2,n+1,n+1bx2,n+1P1-x2,n+1P-x2,n+1Px3,n+1P+x1,n+1Px2,n+1Px3,n+1P,x3,n+1=x3,0+hαΓα+2∑j=0nσ3,j,n+1-cx3,j2+dx1,jx3,j+ex1,jx3,j+hαΓα+2σ3,n+1,n+1-cx3,n+12P+dx1,n+1Px3,n+1P+ex2,n+1Px3,n+1P,where (29)x1,n+1P=x1,0+1Γα∑j=0nρ1,j,n+1ax1,j1-x1,j-x1,jx3,j+x1,jx2,jx3,j,x2,n+1P=x2,0+1Γα∑j=0nρ2,j,n+1bx2,j1-x2,j-x2,jx3,j+x1,jx2,jx3,j,x3,n+1P=x3,0+1Γα∑j=0nρ3,j,n+1-cx3,j2+dx1,jx3,j+ex2,jx3,j.Therefore, for i=1,2,3, (30)σi,j,n+1=nα+1-n-αn+1α,ifj=0n-j+2α+1+n-jα+1-2n-j+1α+1,if1≤j≤n,1,ifj=n+1,ρi,j,n+1=hααn+1-jα-n-jα.We used the step size to be0.05 in Figures 1, 2, and 3. In Figure 1, we used the initial values x1(0)=0.1, x2(0)=0.1, and x3(0)=0.2 and the parameter values α=0.8, a=1.5, b=2.0, c=2.0, d=0.9, and e=0.9, such that equilibrium point E6 is locally asymptotically stable. We have the approximate solutions between x1 and x2 in Figure 1(a), between x1 and x3 in Figure 1(b), and between x2 and x3 in Figure 1(c) and we showed that the solutions are stable.Figure 1
An example for the locally asymptotically stable equilibrium pointE6 where α=0.8, a=1.5, b=2.0, c=2.0, d=0.9, e=0.9.
(a)
The approximate solutions betweenx1t and x3t
(b)
The approximate solutions betweenx1t and x3t
(c)
The approximate solutions betweenx2t and x3t
(d)
The stability of approximate solutionsx1, x2 and x3Figure 2
An example for the locally asymptotically stable of the equilibrium pointE6 where α=0.8, a=1.0, b=2.0, c=2.0, d=0.7, e=0.7.
(a)
The approximate solutions betweenx1t and x2t
(b)
The approximate solutions betweenx1t and x3t
(c)
The approximate solutions betweenx2t and x3t
(d)
The stability of approximate solutionsx1, x2 and x3Figure 3
Example for the locally asymptotically stable equilibrium pointE7, where α=0.9, a=2.5, b=2.0, c=2.0, d=0.8, and e=0.9.
(a)
The approximate solutions betweenx1t and x2t
(b)
The approximate solutions betweenx1t and x3t
(c)
The approximate solutions betweenx2t and x3t
(d)
The stability of approximate solutionsx1, x2, and x3In Figures2 and 3, we used the initial values x1(0)=0.1, x2(0)=0.125, and x3(0)=1.75. In Figure 2, the constant values are α=0.8, a=1.0, b=2.0, c=2.0, d=0.7, and e=0.7. In Figure 3, the constant values are α=0.9, a=2.5, b=2.0, c=2.0, d=0.8, and e=0.9 which satisfied the existence conditions of the locally asymptotically stable equilibrium point E7. We plot the approximate solutions between x1 and x2 in Figures 2(a) and 3(a), between x1 and x3 in Figures 2(b) and 3(b), and between x2 and x3 in Figures 2(c) and 3(c) and we showed that the solutions are stable.We used the trapezoidal method [34] to compare the results and plot some new figures; see Figures 4, 5, and 6 to show this comparison.Figure 4
An example for the locally asymptotically stable equilibrium pointE6, where α=0.8, a=1.5, b=2.0, c=2.0, d=0.9, and e=0.9.
(a)
The approximate solutions betweenx1t and x3t
(b)
The approximate solutions betweenx1t and x3t
(c)
The approximate solutions betweenx2t and x3t
(d)
The stability of approximate solutionsx1, x2, and x3Figure 5
An example for the locally asymptotically stable equilibrium pointE6, where α=0.8, a=1.0, b=2.0, c=2.0, d=0.7, and e=0.7.
(a)
The approximate solutions betweenx1t and x2t
(b)
The approximate solutions betweenx1t and x3t
(c)
The approximate solutions betweenx2t and x3t
(d)
The stability of approximate solutionsx1, x2, and x3Figure 6
Example for the locally asymptotically stable equilibrium pointE7, where α=0.9, a=2.5, b=2.0, c=2.0, d=0.8, and e=0.9.
(a)
The approximate solutions betweenx1t and x2t
(b)
The approximate solutions betweenx1t and x3t
(c)
The approximate solutions betweenx2t and x3t
(d)
The stability of approximate solutionsx1, x2, and x3We used the step size to be0.05 in Figures 4, 5, and 6. In Figure 4, we used the initial values x1(0)=0.1, x2(0)=0.1, and x3(0)=0.2 and the parameter values α=0.8, a=1.5, b=2.0, c=2.0, d=0.9, and e=0.9, such that equilibrium point E6 is locally asymptotically stable. We have the approximate solutions between x1 and x2 in Figure 4(a), between x1 and x3 in Figure 4(b), and between x2 and x3 in Figure 4(c) and we showed that the solutions are stable.In Figures5 and 6, we used the initial values x1(0)=0.1, x2(0)=0.125, and x3(0)=1.75. In Figure 5, the constant values are α=0.8, a=1.0, b=2.0, c=2.0, d=0.7, and e=0.7. In Figure 6, the constant values are α=0.9, a=2.5, b=2.0, c=2.0, d=0.8, and e=0.9 which satisfied the existence conditions of the locally asymptotically stable equilibrium point E7. We plot the approximate solutions between x1 and x2 in Figures 5(a) and 6(a), between x1 and x3 in Figures 5(b) and 6(b), and between x2 and x3 in Figures 5(c) and 6(c) and we showed that the solutions are stable.


## 5. Conclusion

In this paper, we studied the existence, the uniqueness, the stability of the equilibrium points, and numerical solutions of a fractional-order two-prey one-predator model. The coexistence equilibrium pointsE6 and E7 were stable equilibrium points under some conditions at the ordinary differential equation form of the model. But, in our fractional form, we found that the same points are stable without any conditions on the parameters. This means that the predators (lions) can live stably together with the two preys (gazelles and buffalos) without extinction of any of them. This is an example of the equilibrium point which is centre for the integer order system but locally asymptotically stable for its fractional-order counterpart. This means that the fractional-order differential equations are, at least, as stable as their integer order counterpart. So, we recommend to restudy most of the complex models (biological, medical, etc.) in the fractional-order form.



---

*Source: 6714538-2017-08-30.xml*