# Stability Analysis of a Fractional Order Modified Leslie-Gower Model with Additive Allee Effect

Authors: Agus Suryanto; Isnani Darti; Syaiful Anam
Journal: International Journal of Mathematics and Mathematical Sciences
 (2017)

---

## Abstract

We analyze the dynamics of a fractional order modified Leslie-Gower model with Beddington-DeAngelis functional response and additive Allee effect by means of local stability. In this respect, all possible equilibria and their existence conditions are determined and their stability properties are established. We also construct nonstandard numerical schemes based on Grünwald-Letnikov approximation. The constructed scheme is explicit and maintains the positivity of solutions. Using this scheme, we perform some numerical simulations to illustrate the dynamical behavior of the model. It is noticed that the nonstandard Grünwald-Letnikov scheme preserves the dynamical properties of the continuous model, while the classical scheme may fail to maintain those dynamical properties.


---

## Body

## 1. Introduction

The dynamical interaction of predator and prey is one of important subjects in ecological science. In recent years, one of the most important species interactions is predator-prey model [1]. One of well-known mathematical models which describe the dynamics of prey-predator interaction is the modified Leslie-Gower model proposed by Aziz-Alaoui and Okiye [2]. In this model, the growth rate of predator is in the form of logistics-type where its carrying capacity is proportional to the prey number and environment protection for predator. One of important parameters describing the prey-predator interaction is the functional response which describes the predator’s rate of prey consumption per capita. Aziz-Alaoui and Okiye [2] and Yu [3] have considered a modified Leslie-Gower model with Holling type II functional response, while Yu [4] considered the same model but with Beddington-DeAngelis functional response. In the normalized variables, the modified Leslie-Gower equation with Beddington-DeAngelis functional response can be written as (1)dNdt=N1-N-ωNPa+bN+cPdPdt=sP1-PN+k,where N=Nt and P=Pt denote population densities of prey and predator at time t, respectively. The parameters ω, a, b, c, s, and k are positive constants.One of other factors that influence the interaction of predator and prey is Allee effect, referring to a decrease in per capita fertility rate at low population densities. Allee effect may occur under several mechanisms, such as difficulties in finding mates when population density is low or social dysfunction at small population sizes. When such a mechanism operates, the per capita fertility rate of the species increases with density; that is, positive interaction among species occurs [5–8]. Recently Indrajaya et al. [9] investigate a modified Leslie-Gower equation with Beddington-DeAngelis functional response and additive (both weak and strong) Allee effect on prey(2)dNdt=N1-N-mN+h-ωNPa+bN+cPdPdt=sP1-PN+k,with initial conditions N0>0 and P0>0. The criteria of the Allee effect are as follows [7, 8]:(i)
If0<m<h, then the Allee effect is weak.(ii)
Ifm>h, then the Allee effect is strong.It is shown in system (2) that the growth rates of both prey and predator depend only on the current state. In many situations, the growth rate is also dependent on the history of variable or its memory. With the rapid development of fractional calculus, fractional differential equations have been implemented in various fields including biological system. This is due the fact that fractional differential equations are naturally related to the real life phenomena with memory which exists in most of biological system [10–16]. To describe such memory effect, we first recall the definition of fractional integral operator as well as fractional differential operator.Definition 1 (see [17]).
The Riemann-Liouvilleα-order fractional integral operator of any function u∈L10,a, x∈0,a is defined by(3)Jαux=1Γα∫0ax-tα-1utdt,where Γ· is the Euler Gamma function.Definition 2 (see [17]).
Letm be an integer which satisfies m-1<α<m. The Riemann-Liouville α-order fractional derivative of function u∈L10,a is defined as(4)DRLαux≔dmdxmJm-αux=1Γm-αdmdxm∫0ax-tm-α-1utdt,where dm/dxm is the common m-order derivative.The Riemann-Liouville fractional derivative is historically the first concept of fractional derivative and theoretically well established. However, in the case of Riemann-Liouville fractional differential equation, the initial value is usually given in the form of fractional derivative, which is not practical. Consequently, one applies the Caputo fractional derivative which is defined as follows.Definition 3 (see [17]).
The Caputo fractional differential operator of orderα>0, with m-1<α<m, m∈N, is defined by(5)DCαux≔Jm-αumx=1Γm-α∫0ax-tm-α-1umtdt.For simplicity, the Caputo fractional derivative of function ut of order α is denoted by dαut/dtα.From Definition3, we see the α-order fractional derivative at timet is not defined locally; it relies on the total effects of the commonly usedm-order integer derivative on the interval 0,t. So it can be used to describe the variation of a system in which the instantaneous change rate depends on the past state, which is called the ‘‘memory effect’’ in a visualized manner [18].In this paper we reconsider system (2). By assuming that the growth rates of both prey and predator at time t do not only depend instantaneously on the current state but also depend on the past state, we replace the first order derivatives in system (2) with the fractional order Caputo type derivatives:(6)dαNdtα=N1-N-mN+h-ωNPa+bN+cPdαPdtα=sP1-PN+k,where 0<α<1. Hence we have a system of fractional differential equation. In the following we discuss the dynamical properties of system (6). To study the stability of equilibrium points, we apply the following stability theorem.Theorem 4 (see [19]).
Consider the following autonomous nonlinear fractional order system:(7)dαu→tdtα=f→u→t;u→0=u→0;0<α<1.The equilibrium points of the above system are solutions to the equation f→u→t=0. An equilibrium point u→∗ is locally asymptotically stable if all eigenvalues λj of the Jacobian matrix J=∂f→/∂u→ evaluated at equilibrium u→∗ satisfy argλj>απ/2.


## 2. Equilibria and Their Stability

Based on Theorem4, equilibria of model (6) can be determined by solving the following system:(8)N1-N-mN+h-ωNPa+bN+cP=0sP1-PN+k=0.System (8) has been solved by Indrajaya et al. [9], and the obtained equilibria are as follows:(1)
The equilibria of system (6) for the weak Allee effect case are
(a)
trivial equilibriumEw0=0,0, that is, the extinction of both prey and predator point,(b)
two axial equilibria, that is, the prey extinction pointEw1=0,k and the predator extinction point Ew2=Nw2,0 where (9)Nw2=121-h+1-h2+4h-m,(c)
positive or coexistence equilibriumEw∗=Nw∗,Pw∗=Nw∗+k, where Nw∗ are all possible real positive solutions of cubic equation(10)N3+3η1N2+3η2N+η3=0,where(11)η1=-b+c1-h-a+ck-ω3b+c,η2=-a+ck1-h+b+ch-m-ωh+k3b+c,η3=-a+ckh-m-hωkb+c.(2)
The equilibria of system (6) for the strong Allee effect case are
(a)
trivial equilibriumEs0=0,0, that is, the extinction of both prey and predator point,(b)
three axial equilibria that are the prey extinction pointEw1=0,k and two predator extinction points: Es2=Ns2,0 and Es3=Ns3,0, where (12)Ns2=121-h+1-h2+4h-m,Ns3=121-h-1-h2+4h-m,(c)
positive or coexistence equilibrium pointEs∗=Ns∗,Ps∗=Ns∗+k, where Ns∗ are also all possible real positive solutions of cubic equation (10).Using transformation z=N∗+η1, (10) can be reduced to(13)hz=z3+3pz+q=0,where p=η2-η12 and q=η3-3η1η2+2η1. Implementing Cardan’s method as performed by Cai et al. [7], we obtain the existence of the positive equilibria as follows.Lemma 5 (existence of positive equilibria).
Let(N∗,P∗) be the interior equilibrium of model (6) for both weak and strong Allee effects where N∗is a real positive solution of (10). Then the following statements hold:(a)
Ifq<0, then (13) has a single positive root z1. As a result, model (6) has a unique positive equilibrium point, that is, E∗=N∗,N∗+k=z1-η1,z1-η1+k, with z1>η1.(b)
Suppose thatq>0 and p<0, then
(b1)
Ifq2+4p3=0, then (13) has a positive root of multiplicity two. Thus, model (6) has a unique positive equilibrium point, that is, E∗=N∗,N∗+k=-p,-p+k.(b2)
Ifq2+4p3<0, then (13) has two positive roots z1 and z2. Thus, model (6) has two positive equilibrium points, namely, E1∗=N1∗,N1∗+k=z1-η1,z1-η1+k and E2∗=N2∗,N2∗+k=z2-η1,z2-η1+k, with E1∗=N1∗,N1∗+kz1,2>η1.(c)
Ifq=0 and p<0, then (13) has a unique positive root z1=-3p. As a result, model (6) has a unique positive equilibrium point, that is, E∗=N∗,N∗+k=-3p,-3p+k with qm=0=0.Moreover, algebraic computations show that if (13) has two positive roots, then they are(14)z1=-4q+44p3+q22/3-4p2-4q+44p3+q22/3,z2=-z12+z13+4q2z1.If (13) has a positive root, then it must be (15)z1=-4q+44p3+q22/3-4p2-4q+44p3+q22/3.To check the local stability of each equilibrium point, we linearize system (6) around the equilibrium and verify all eigenvalues of the Jacobian matrix evaluated at the equilibrium. The stability properties of trivial and axial equilibrium points for the case of weak and strong Allee effect are, respectively, stated in Theorems 6 and 7.Theorem 6.
Stability of trivial and axial equilibrium for weak Allee effect(0<m<h):(i)
the trivial equilibriumEw0=0,0 and the axial equilibrium Ew2=Nw2,0 are always unstable;(ii)
the axial equilibriumEw1=0,k is asymptotically stable if h-m<hωk/(a+ck).Proof.
(i) The Jacobian matrix atEw0 is JEw0=1-m/h00s, and the eigenvalues are λ1=1-m/h>0 and λ1=s>0. It is clear that arg⁡λ1,2=0<απ/2,α>0. Hence Ew0 is unstable. The Jacobian matrix at Ew2 is JEw2=1-2Nw2∗-mh/2Nw2∗+h2-ωNw2∗/(a+bNw2∗)0s where one of its eigenvalues is λ=s>0 and therefore Ew2 is unstable because arg⁡λ=0<απ/2, α>0.
(ii) The Jacobian matrix atEw1 is JEw1=1-m/h-ωk/(a+ck)00-s. The eigenvalues of JEw1 are λ1=-s<0 and λ2=1-m/h-ωk/(a+ck). Thus arg⁡λ1=π>απ/2 and arg⁡λ2=π>απ/2 whenever h-m<hωk/(a+ck). This proves part (ii).Using the same argument as in the proof of Theorem6, we obtain the following stability properties of equilibria for the case of strong Allee effect.Theorem 7.
Stability of trivial and axial equilibrium for strong Allee effect(m>h): (i)
the trivial equilibriumEs0=0,0; the axial equilibria: Es2=Ns2,0 and Es3=Ns3,0 are always unstable;(ii)
the axial equilibriumEs1=0,k is always asymptotically stable.The stability properties of positive (coexistence) equilibrium for the case of weak and strong Allee effect are stated in Theorems8 and 9.Theorem 8.
Stability of coexistence equilibrium for weak Allee effect(0<m<h):
supposeJEw∗ is the Jacobian matrix at coexistence equilibrium Ew∗. Equilibrium Ew∗ is asymptotically stable if one of the following mutually exclusive conditions holds:(i)
TraceJEw∗<0; DetJEw∗>0 and Δ=Trace2JEw∗-4DetJEw∗≥0.(ii)
DetJEw∗>0, Δ<0 and (16)ΔTraceJEw∗>tan⁡απ2.Proof.
The characteristics equation ofJEw∗ is given by (17)λ2-TraceJEw∗λ+DetJEw∗=0.(i)
IfTraceJEw∗<0; DetJEw∗>0; and Δ≥0 then λ1,2<0; hence argλ1,2=π>απ/2 and the result follows.(ii)
Ifλ is an eigenvalue of JEw∗ and Δ<0, then λ- is also an eigenvalue. Using Δ/TraceJEw∗>tanαπ/2, we have that (λ-λ-)/(λ+λ-)=Imλ/Re⁡λ=argλ>tan⁡απ/2. Therefore the stability of Ew∗ follows.Similarly we have Theorem9 for the stability of coexistence equilibrium for strong Allee effect case.Theorem 9.
Stability of coexistence equilibrium for strong Allee effect(m>h):
supposeJEs∗ is the Jacobian matrix evaluated at the coexistence equilibrium Es∗. Equilibrium Es∗ is asymptotically stable if one of the following mutually exclusive conditions holds:(i)
TraceJEs∗<0; DetJEs∗>0; and (18)Δ=Trace2JEs∗-4DetJEs∗≥0.(ii)
DetJEs∗>0, Δ<0, and(19)ΔTraceJEs∗>tan⁡απ2.Based on the above theorems it can be seen that the stability properties of both trivial and axial equilibrium points are not dependent onα (order of fractional derivative). But α may influence significantly the stability of coexistence equilibrium point. Coexistence point can be asymptotically stable although the eigenvalue of Jacobian matrix has positive real part, provided that conditions of Theorem 8(ii) or Theorem 9(ii) are met. This is in contrast to the coexistence equilibrium point of the integer-order model (2) where coexistence point is asymptotically stable only if all real parts of the eigenvalues of the Jacobian matrix are negative [9].


## 3. Numerical Simulations

To solve system (6), we implement a nonstandard Grünwald-Letnikov scheme which is a combination of the Grünwald-Letnikov approximation [20] and the nonstandard finite difference (NSFD) method [21, 22]. According to [20], the explicit (or implicit) Grünwald-Letnikov (GL) approximation for a fractional differential equation with initial value(20)dαytdtα=fyt,y0=y00<α<1is given by(21)ytn+1-∑v=1n+1cvαytn+1-v-rn+1αy0=Δtαfynor=Δtαfyn+1,where cvα=1-(α+1)/vcv-1α; c1α=α; and rn+1α=Δtαr0αtn+1=n+1-α/Γ1-α. Here, Δt represents the time step of numerical integration. The Grünwald-Letnikov approximation is proceeding iteratively but the sum in the scheme becomes longer and longer which represents the memory effects. Scherer et al. [20] have shown that the coefficient cvα is positive and satisfies 0<cn+1α<cnα<⋯<c1α=α for n≥1. Observe that in the standard Grünwald-Letnikov approximation (21), the right hand side of (20) is approximated locally. We implement a nonstandard method which is adopted from the NSFD method [23]. A numerical scheme for an initial value problem (22)du→dt=f→t,u→;u→0=u→0is called a NSFD method if at least one of the following conditions is satisfied [21, 22]:(i)
The left hand side is approximated by the generalization of forward difference scheme(23)du→ndt≈u→n+1-u→nψΔt.The nonnegative denominator function has to satisfyψΔt=Δt+OΔt2.(ii)
The approximation offt,u→ is nonlocal.By implementing the Grünwald-Letnikov approximation for the fractional derivative and the nonlocal approximation for the right hand side of system (6) we get the following scheme:(24)Nn+1=rn+1αN0+∑j=1n+1cjαNn+1-j+ΔtαNn1+ΔtαNn+m/Nn+h+ωPn/a+bNn+cPnPn+1=rn+1αP0+∑j=1n+1cjαPn+1-j+ΔtαsPn1+ΔtαsPn/Nn+k.Observe that scheme (24) is explicit and hence it is simple and easy to be implemented. Besides that, the nonstandard Grünwald-Letnikov scheme (24) also maintains the positivity solutions.To verify our stability analysis as well as the effectiveness our numerical scheme, we perform some numerical simulations. First we use hypothetic values of parameterss=1, k=0.5, m=0.2, h=0.1, ω=0.7, a=1, b=0.5, c=0.05, and Δt=0.1. Model (6) with these parameters has four equilibrium points: (0,0), (0.1298,0), (0.7701,0), and (0,0.5). According to Theorem 7, only axial equilibrium (0,0.5) is stable for any order of fractional derivative α, where 0<α<1. This stability behavior is confirmed by our numerical solutions; see Figures 1 and 2. This shows that strong Allee effect may lead to an extinction of prey population. It is shown in Figure 2 that our numerical solutions for α=0.5,0.6,0.7,0.8,0.9,1.0 are convergent to the axial equilibrium (0,0.5), indicating that equilibrium (0,0.5) is stable asymptotically for any order of fractional derivative. Detail observation shows that as the order of fractional derivative increases the convergence of solution is faster and the solution of system (6) closes to the integer-order model (2).Figure 1
Phase portrait of system (6). The values of parameters are s=1, k=0.5, m=0.2, h=0.1, ω=0.7, a=1, b=0.5, c=0.05, and Δt=0.1.
(a)
α
=
0.80
(b)
α
=
0.95Figure 2
Solutions of system (6) with initial value (1, 0.2) and parameter values: s=1, k=0.5, m=0.2, h=0.1, ω=0.7, a=1, b=0.5, c=0.05, and Δt=0.1, for various values of α. All numerical solutions are convergent to axial equilibrium (0, 0.5). Solution of system (6) with larger order of fractional derivative α has faster convergence compared to that with smaller α.
(a)
(b)Next, we sets=0.02 and h=0.3 and keep the rest of parameters as in Figure 1. This weak Allee case has a trivial equilibrium (0,0); two axial equilibria: (0.8217,0) and (0,0.5); and two interior points: Ew1∗=0.0128,0.5128 and Ew2∗=0.1373,0.6373. Theorems 7 and 9 state that axial equilibrium (0,0.5) is locally stable for 0<α<1 and interior point Ew2∗=0.1373,0.6373 is locally stable if α<α∗=0.886, and other equilibria are always unstable. Such behavior is in accordance with our numerical results depicted in Figures 3(a) and 3(b). It is clearly seen in Figure 3(a) that α=0.8 produces bistable dynamic where depending on the initial values, solutions may be convergent to the extinction of prey point (0,0.5) or to interior point. In other words, the solution of system (6) is highly sensitive to the initial conditions. An initially relatively small prey will converge to the prey extinction point. On the other hand, if the prey is initially relatively large then prey and predator will coexist. If we increase the order of derivative such that α=0.95, the axial equilibrium (0,0.5) is still locally stable but the interior point becomes unstable; see Figure 3(b). In latter case, there exists a stable limit cycle which shows that both prey and predator are fluctuating around the interior point. However, the appearance of limit cycle may be suppressed by increasing the coefficient of predator interference. For example we plot in Figure 4(a) the numerical solution using the same parameters as in Figure 3(b) but with c=0.1. We see that the interior point is now stable while the axial equilibrium point (0,0.5) is unstable. It can be said that relatively large predator interference can stabilize the interior point. On the other hand, strong Allee effect can destabilize or even remove the interior equilibrium and can cause the extinction of prey population. For example, we show numerical simulation using parameters the same as in Figure 3(b) except m=0.3. This simulation shows that all initial values converge to the axial equilibrium (0,0.5) which shows the extinction of prey population but predator species can still survive in the habitat because there is an enough environmental protection; see Figure 4(b).Figure 3
Phase portrait of system (6). The values of parameters are s=0.02, k=0.5, m=0.2, h=0.3, ω=0.7, a=1, b=0.5, c=0.05, and Δt=0.1.
(a)
α
=
0.80
(b)
α
=
0.95Figure 4
Phase portrait of system (6). The values of parameters are s=0.02,k=0.5,h=0.3, ω=0.7, a=1, b=0.5, α=0.95, and Δt=0.1.
(a)
m
=
0.2; c=0.1
(b)
m
=
0.3; c=0.05Finally, we compare our numerical results obtained by the NSGL scheme to those obtained by the standard GL scheme using parameterss=0.02, k=0.5,m=0.2, h=0.3, ω=0.7, a=1, b=0.5, c=0.1, and α=0.95; see Figure 5. We see that both NSGL and GL schemes using Δt=0.005 produce solutions where their difference cannot be observed in the scale of Figure 5. Using Δt=0.01, the numerical solutions of both NSGL and GL schemes are in excellent agreement with those of both schemes using Δt=0.005. If we take time step Δt=0.1, both schemes have comparable solutions which are initially distorted from solutions with much smaller time step; see Figure 5(a). In Figure 5(b), we plot solutions using relatively large time step (Δt=2.0 and Δt=2.25). Although the NSGL scheme has solutions which are quantitatively different from solution with Δt=0.005, nevertheless those solutions still have the same behavior as before; that is, they are always positive and convergent to the correct equilibrium point. However, the GL scheme in this case gives unrealistic negative value for prey population. If the time step is further increased, the GL scheme will be unstable and leads to blowing up solutions.Figure 5
Phase portrait of system (6) calculated using nonstandard (NSGL) and standard (GL) schemes. The values of parameters are s=0.02, k=0.5, m=0.2, h=0.3, ω=0.7, a=1, b=0.5, c=0.1, and α=0.95.
(a)
(b)


## 4. Conclusion

The dynamic of a fractional order modified Leslie-Gower model with Beddington-DeAngelis functional response and additive Allee effect has been analyzed. Our model has four types of equilibria that are the trivial (extinction of both prey and predator) equilibrium, two axial equilibria (the prey extinction point and the predator extinction point), and the interior (coexistence) point. The trivial and the predator extinction for both weak and strong Allee effects are always unstable. For the case of weak Allee effect, the prey extinction is conditionally stable while for that of strong Allee effect, the prey extinction is always stable. Our analysis also shows that the interior point for both weak and Allee effects is conditionally stable. The order of fractional derivative may influence the stability of interior point. Here, when the orderα is larger than critical order α∗, then the interior point may be destabilized. These dynamical properties are confirmed by our NSGL schemes which shows the effectiveness of NSGL scheme. It is also shown that the NSGL scheme preserves the positivity of numerical solutions. Furthermore, our numerical results show that the NSGL scheme produces numerical solutions which satisfy the dynamical behavior of our model. However, the standard GL scheme may fail to preserve such properties; for example, it can produce nonrealistic negative solutions.



---

*Source: 8273430-2017-05-28.xml*