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

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ünwaldLetnikov scheme preserves the dynamical properties of the continuous model, while the classical scheme may fail to maintain those dynamical properties.


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 predatorprey 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 where  = () and  = () denote population densities of prey and predator at time , respectively.The pa rameters , , , , , and  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][6][7][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 with initial conditions (0) > 0 and (0) > 0. The criteria of the Allee effect are as follows [7,8]: (i) If 0 <  < ℎ, then the Allee effect is weak.
(ii) If  > ℎ, 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][11][12][13][14][15][16].To describe such memory effect, we first recall the definition of fractional integral operator as well as fractional differential operator.
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.
From Definition 3, we see the -order fractional derivative at time t is not defined locally; it relies on the total effects of the commonly used m-order integer derivative on the interval [0, ].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  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: 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: The equilibrium points of the above system are solutions to the equation ⃗ ( ⃗ ()) = 0.An equilibrium point ⃗  * is locally asymptotically stable if all eigenvalues (  ) of the Jacobian matrix

Equilibria and Their Stability
Based on Theorem 4, equilibria of model ( 6) can be determined by solving the following system: 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 equilibrium  0 = (0, 0), that is, the extinction of both prey and predator point, (b) two axial equilibria, that is, the prey extinction point  1 = (0, ) and the predator extinction point  2 = ( 2 , 0) where (c) positive or coexistence equilibrium  *  = ( *  ,  *  =  *  + ), where  *  are all possible real positive solutions of cubic equation where (2) The equilibria of system (6) for the strong Allee effect case are (a) trivial equilibrium  0 = (0, 0), that is, the extinction of both prey and predator point, (b) three axial equilibria that are the prey extinction point  1 = (0, ) and two predator extinction points:  2 = ( 2 , 0) and  3 = ( 3 , 0), where (c) positive or coexistence equilibrium point  *  = ( *  ,  *  =  *  +), where  *  are also all possible real positive solutions of cubic equation (10).
Moreover, algebraic computations show that if (13) has two positive roots, then they are If ( 13) has a positive root, then it must be 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.
Using the same argument as in the proof of Theorem 6, we obtain the following stability properties of equilibria for the case of strong Allee effect.
The stability properties of positive (coexistence) equilibrium for the case of weak and strong Allee effect are stated in Theorems 8 and 9.
Theorem 8. Stability of coexistence equilibrium for weak Allee effect (0 <  < ℎ): suppose ( *  ) is the Jacobian matrix at coexistence equilibrium  *  .Equilibrium  *  is asymptotically stable if one of the following mutually exclusive conditions holds: Similarly we have Theorem 9 for the stability of coexistence equilibrium for strong Allee effect case.

Theorem 9. Stability of coexistence equilibrium for strong Allee effect (𝑚 > ℎ):
suppose ( *  ) is the Jacobian matrix evaluated at the coexistence equilibrium  *  .Equilibrium  *  is asymptotically stable if one of the following mutually exclusive conditions holds: (i) Trace(( *  )) < 0; Det(( *  )) > 0; and (ii) Det(( *  )) > 0, Δ < 0, and 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].

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 is given by where Here, Δ 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   V is positive and satisfies 0 <   +1 <    < ⋅ ⋅ ⋅ <   1 =  for  ≥ 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 is 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 The nonnegative denominator function has to satisfy (Δ) = Δ + (Δ 2 ).
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: 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.
Next, we set  = 0.02 and ℎ = 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:  * 1 = (0.0128, 0.5128) and  * 2 = (0.1373, 0.6373).Theorems 7 and 9 state that axial equilibrium (0, 0.5) is locally stable for 0 <  < 1 and interior point  * 2 = (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  = 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  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).
Finally, we compare our numerical results obtained by the NSGL scheme to those obtained by the standard GL scheme using parameters  = 0.02,  = 0.5,  = 0.2, ℎ = 0.3,  = 0.7,  = 1,  = 0.5,  = 0.1, and  = 0.95; see Figure 5.We see that both NSGL and GL schemes using Δ = 0.005 produce solutions where their difference cannot be observed in the scale of Figure 5.Using Δ = 0.01, the numerical solutions of both NSGL and GL schemes are in excellent agreement with those of both schemes using Δ = 0.005.If we take time step Δ = 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 (Δ = 2.0 and Δ = 2.25).Although the NSGL scheme has solutions which are quantitatively different from solution with Δ = 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.