A Modiﬁed Leslie–Gower Model Incorporating Beddington–DeAngelis Functional Response, Double Allee Effect and Memory Effect

: In this paper, a modiﬁed Leslie–Gower predator-prey model with Beddington–DeAngelis functional response and double Allee effect in the growth rate of a predator population is proposed. In order to consider memory effect on the proposed model, we employ the Caputo fractional-order derivative. We investigate the dynamic behaviors of the proposed model for both strong and weak Allee effect cases. The existence, uniqueness, non-negativity, and boundedness of the solution are discussed. Then, we determine the existing condition and local stability analysis of all possible equilibrium points. Necessary conditions for the existence of the Hopf bifurcation driven by the order of the fractional derivative are also determined analytically. Furthermore, by choosing a suitable Lyapunov function, we derive the sufﬁcient conditions to ensure the global asymptotic stability for the predator extinction point for the strong Allee effect case as well as for the prey extinction point and the interior point for the weak Allee effect case. Finally, numerical simulations are shown to conﬁrm the theoretical results and can explore more dynamical behaviors of the system, such as the bi-stability and forward bifurcation.


Introduction
Modeling interaction between prey and its predator has become a dominant topic in mathematical biology due to its ubiquitous existence and fundamentality in many biological systems. Study of the dynamics of the predator-prey model could provide qualitative explanations of numerous phenomena that can occur in predator and prey interaction. One of the crucial phenomenon in ecology that influences the per capita growth rate either in the predator or prey population is the Allee effect, which describes a condition where, at low population densities, the per capita growth rate of the population has a positive dependence with its density. There are two kinds of Allee effects, namely the strong Allee effect and the weak Allee effect. In the strong Allee effect, there is a population threshold value named the Allee threshold, below which the per capita growth rate of the population is negative [1,2]. In terms of conservation biology, if the Allee threshold is larger, then it places a population at higher risk of extinction in a low-density population. Meanwhile, in the weak Allee effect, the per capita growth rate of the population always remains positive but is still reduced at low population densities [3][4][5]. If two or more mechanisms that generate the Allee effect works simultaneously on a single population, then it is known as the double (or multiple) Allee effect [6,7]. Biological evidence of such phenomenon from both terrestrial and aquatic habitat is given in Table 2 of [7,8] and the references cited therein.
From the mathematical point of view, some scholars have been investigating the dynamics of predator-prey models with the Allee effect in the prey population [9][10][11][12][13] or predator population [14][15][16]. The main focus of all the mentioned research studies is to investigate whether the Allee effect has a tremendous impact on the occurrence of various dynamics in the predator-prey model. For the double Allee effect phenomenon, there are some papers that study the double Allee effect in the prey, see for example [17,18]. However, most of the studies just focused on the double Allee effect in the growth of prey population, although observations are showing that the double Allee effect could be discovered in the growth of the predator population. A typical example comes from the endangered species African wild dog (Lycaon pictus). Their social system requires a high population density to survive and reproduce. Being a predator, the African wild dog is a generalist species with the Thomson's gazelle (Eudorcas thomsonii) as their common prey but it also hunts other animals such as the impala (Aepyceros melampus), warthog (Phacochoerus aethiopicus), hares, etc. They live in permanent packs of about 27 adults and pups and have to share food after killing their prey. There is also interference among predators in their hunting behavior. We refer the readers to [19] for details.
In the natural world, the presence of memory must exist in prey and predator interaction since the growth rates of prey and predator at any point depend on the history of the variables at all previous times and not only on the current state which is local to that point [20][21][22][23][24]. Recently, fractional calculus through the fractional derivatives has been known to provide an excellent instrument for describing the memory and hereditary properties of various materials and processes [25], such as in biology, finance, engineering, and physics (see, for example, [26][27][28][29][30] and the references therein). Some interesting papers regarding the Allee effect in the fractional-order predator-prey models are provided in [31,32]. In [31], Suryanto et al. have studied the local stability of the modified Leslie-Gower model with the Beddington-DeAngelis functional response and additive Allee effect in the prey population. They construct the numerical scheme that preserves the dynamics of its first-order system provided by [33]. Later, in [32], Baisad and Moonchai considered the Gause predator-prey model that includes the Allee effect in the prey population and Holling type-III functional response. They also studied the local stability and sufficient conditions of a Hopf bifurcation at the positive equilibrium point. In both papers, the dynamical behaviors are influenced by the order of the derivative.
Motivated by the above mentioned points, this paper aims to study the fractional-order Leslie-Gower predator-prey model incorporating the Beddington-DeAngelis functional response and double Allee effect. The proposed model includes the Caputo fractional-order derivative to capture the effect of memory in the growth rates of both prey and predator. From what we know, the dynamic of our proposed model that incorporates the double Allee effect in the growth of predator and memory effect under the Caputo fractional-order derivative has not been proposed and investigated by other scholars. This work may reveal an interesting ecological point of view to the importance of the double Allee effect than the single Allee effect towards the management of exploited or threatened predator population.
The remaining part of this paper is organized as follows. In Section 2, the model formulation is given. The existence, uniqueness, non-negativity, and boundedness of solutions of our model are discussed in Section 3. Then, we investigate the dynamic behaviors of the model for both weak and strong Allee effects. In Section 4, the existence of non-negative equilibrium points and the local stability of non-negative equilibrium points along with Hopf bifurcation analysis are presented. Next, the sufficient conditions for the global stability of the equilibrium points are carried out in Section 5. Numerical simulations are shown in Section 6 to verify our analytical findings as well as to numerically explore the impact of capturing rate, the Allee threshold, and the order of the fractional-order system on the dynamics of our model. Finally, we draw conclusions in Section 7.

Model Formulation
One of the primary directions in modelling the interaction of prey-predator populations is based on the mass conservation principle, which says that the predators can grow only as a function of what they have consumed. Under this principle, the general model of the predator-prey dynamics takes the following model [34]: where x(t) and y(t) are, respectively, the prey and predator population densities at time t, f (x) is the per capita growth rate of prey, g(x, y) is the functional response, k is the predation efficiency, and µ is the predator per capita death rate. An alternative to the conservative predator-prey model (1) is to abandon the mass conservation principle. This type of model does not explicitly describe the relationship between predation rate and the reproduction rate of predator. A foremost predator-prey model in this direction is the Leslie-Gower model [35]. The Leslie-Gower model maintains the prey equation as in system (1) but applies a logistic type of model for the predator equation. By considering that the per capita growth rate of prey obeys the logistic growth and that predation follows the Beddington-DeAngelis functional response, we have the following Leslie-Gower model.
Notice that the logistic type forms in both prey and predator equations are written in the form as suggested by [36]. This typical logistic form is for avoiding paradoxes in the logistic equation [37,38]. All parametersr,β 1 ,b,ĉ,q,ŝ,σ in the system (2) are real and positive. The ecological meaning of the parameters are as follows:r andβ 1 are the intrinsic growth rate of the prey and the prey intraspecific competition coefficient in the absence of predation, respectively;ŝ andσ are the intrinsic growth rate of the predator and the predator intraspecific competition coefficient, respectively;b andĉ measure the effect of capturing rate and handling time by the predator to the predation rate, respectively;q is the strength of interference among predators. The coefficient of predator intraspecific competition is assumed to depend on the prey density, i.e.,σ =β 2 x , whereβ 2 is the constant of predator intraspecific competition. Such assumption makes sense because when the prey is available in abundant (x → ∞), then there will no intraspecific competition (σ → 0) and the predator can attain its maximum per capita growth rateŝ. On the contrary, if the prey is rare (x → 0), then the intraspecific competition becomes maximum (σ → ∞) and, hence, the predator will become extinct as the per capita growth rate of predator becomes −∞. When a severe scarcity of prey occurs, the predator can switch to alternative populations, which causes the reduction in intraspecific competition for hunting the favourite food x.
To account for such phenomenon, Aziz-Alaoui and Okiye [39] proposed a modified Leslie-Gower model by introducing an inhibition coefficient (l) in the intraspecific competition due to the availability of alternative food for the predator. The intraspecific competition coefficient now becomesσ =β 2 l+x . The modified Leslie-Gower model with the Beddington-DeAngelis functional response has been studied in [40,41]. Today, the Leslie-Gower type model still attracts many scholars (see [42][43][44] and the references therein).
In this paper, we reconsider the modified Leslie-Gower predator-prey model (2) but we assume that the intrinsic growth rate of the predator population is affected by the double Allee effect. Then the predator-prey model (2) takes the following form: wherem is the Allee threshold,n is the auxiliary Allee effect constant withn > 0 and n > −m. In the second equation of model (3), the intrinsic growth rate of the predatorŝ is affected by double Allee effects. Without the intraspecific competition for the predator, the per capita growth rate of the predator is reduced fromŝ toŝ y−m y+n due to the Allee effect [45]. Therefore, the Allee effect is strong ifm > 0 and weak if −n <m ≤ 0 [9,18].
In order to seize the entire time population growth condition, we consider a fractionalorder derivative to the left-hand side of the classical derivative system (3) as follows: with initial conditions x(0) > 0 and y(0) > 0. D α * represents the Caputo fractional-order derivative of a real valued function f , which is defined by the following: where Γ(·) is the Gamma function and α ∈ (n − 1, n], n ∈ Z + [25]. In order to overcome the inconsistency of time dimension between the left-hand side of system (4) with its right-hand side, we follow [46][47][48] to modify all of the biological parameters in the right-hand side that have time dimension (time −1 ). Thus, we have a new system.
For simplification, we replace system (5) with the redefined parameters as follows: wherê Notice that, the authors [49] have studied the local stability and have shown numerically a Hopf bifurcation circumstance around the interior point for the case of m = 0. In this paper, we focus on the local and global dynamics of system (6) for both m > 0 and m < 0 cases.

Preliminaries Results
In this section, we bring out the fact that system (6) is biologically well-behaved by showing that the solution of system (6) exists and is unique as well being non-negative and bounded.
For any Z = (x, y),Z = (x,ȳ), Z,Z ∈ Ω, it follows from (7) that the following is the case: Since F(Z) satisfies the Lipschitz condition with respect to Z, it follows from Theorem 3.7 in [50] that there exists a unique solution Z(t) of system (6) with initial condition Z(0) = (x(0), y(0)). Consequently, we have the following theorem. Theorem 1. For each Z(0) = (x(0), y(0)) ∈ Ω, then initial value problem of system (6) has a unique solution Z(t) ∈ Ω which is defined for all t ≥ 0.

Non-Negativity and Boundedness
In order to prove that all solutions of system (6) are non-negative and bounded, let be the non-negative quadrant on the xy− plane. In the case of the biological significance, we must ensure that when the initial condition starts in R 2 + , then the solution of system (6) remains in R 2 + for all t ≥ t 0 .
Theorem 2. If x(t 0 ) ≥ 0 and y(t 0 ) ≥ 0, then all solutions of the system (6) are non-negative and uniformly bounded.
for t ≥ t 0 be the solutions of system (6). Suppose that assumption is false, then there exists t * > t 0 such that (6), one has the following.
Next, we prove the boundedness of all solutions of system (6). From the first equation of system (6), we have the following.
By Lemma 3 in [51], we have the following: where E α is the Mittag-Leffler function. Therefore, x(t) with initial condition x(t 0 ) are confined to the region Ω 1 where the following is the case.
From the second equation of system (6), we have the following.
Again by using Lemma 3 in [51], for the strong Allee effect (m > 0), we have the following.
However, for the weak Allee effect (m < 0), we have the following.
Therefore, the solution of y(t) with initial condition y(t 0 ) are confined to region and where the following is the case.

Equilibrium Points and Their Local Stability
In this section, the equilibrium points and existence conditions are obtained and their local stability is analyzed by using the Matignon condition [25] for the weak (m < 0) and the strong (m > 0) Allee effect, respectively.
Denote W 2 as always existing.
positive roots of the quartic equation (14): where (2) The equilibrium points of system (6) for the strong Allee effect (m > 0) are as follows: (a) Both prey and predator extinction point S 0 = (0, 0), which always exists; (b) The predator extinction point S 1 = ( r β 1 , 0), which always exists; (c) The prey extinction point S 2,3 = (0,ȳ s ) whereȳ s is the positive solution of the quadratic equation y 2 − ls β 2 − n y + mls β 2 = 0. The existence of S 2,3 is described as follows: , then there exist two prey extinction points, i.e., the following is the case.
andŷ s are also all positive roots of the quartic Equation (14).
In order to study the local stability of system (6) around an equilibrium point (x * , y * ), we consider the following Jacobian matrix J of system (6), which is given by the following.
Theorem 3. The stability properties of trivial and axial equilibrium points of system (6) for the weak Allee effect (m < 0) are as follows: Proof.
(a) By substituting W 0 to (15), we obtain the Jacobian matrix.
The eigenvalues of (18) are as follows.
Theorem 4. Suppose that m < 0 (weak Allee effect) and the following is the case.
Proof. By evaluating (15) at the interior equilibrium pointŴ = (x w ,ŷ w ), we obtain the following.

Proof.
(a) By substituting S 0 to (15), we obtain the following.
The eigenvalues of (22) are as follows.
Proof. The Jacobian matrix (15) at interior equilibrium pointŜ = (x s ,ŷ s ) is given by the following.
Hopf bifurcation on a fractional-order system occurs when the Jacobian matrix evaluated at an equilibrium point has two complex conjugate eigenvalues and there is a limit-cycle when the stability of that system changes. Here, we use the conditions for the existence of a Hopf bifurcation which was introduced by [53]. According to Theorems 4 and 6, the stability of the interior equilibrium point for both weak and strong Allee effects is influenced by the order of the fractional derivative (α). Thus, we can establish the condition for the existence of a Hopf bifurcation at the interior point as α passes through the critical value α * in the following theorem.

System with Weak Allee Effect
We know from the previous analysis that in the case of the weak Allee effect, the prey extinction point W 2 = (0,ȳ w ) and the interior pointŴ = (x w ,ŷ w ) are conditionally locally asymptotically stable. In the following, we study the global asymptotic stability of those equilibrium points.
Proof. We consider the following positive definite Lyapunov function.
Calculating the α-order derivative of V 1 (x, y) along the solution of system (6) and applying Lemma 3.1 in [54], we obtain the following.
, we obtain the following.
In this case, D α * V 1 (x, y) ≤ 0 for all (x, y) ∈ R 2 + and D α * V 1 (x, y) = 0 implies that y =ȳ w . Substituting y =ȳ w to the second equation of system (6) we obtain the following.
The unique solution of (25) is x = 0, which shows that singleton {W 2 } is the only invariant set on which D α * V 1 (x, y) = 0. By Lemma 4.6 in [55], it is proven that W 2 is globally asymptotically stable. (6) is globally asymptotically stable if the following is the case:

Theorem 9. The interior pointŴ of model
Proof. Consider the following positive definite Lyapunov function.
By taking the α-order derivative of V 2 (x, y) along the solution of system (6) and applying Lemma 3.1 in [54], one has the following.

System with Strong Allee Effect
Previous analysis shows that for the system with a strong Allee effect, the predator extinction point S 1 = ( r β 1 , 0) is always locally asymptotically stable. Hence, in the following, we study the global asymptotic stability only for the equilibrium point S 1 = ( r β 1 , 0). Theorem 10. If m ≥ br β 1 s + 1 (γ 2 + n), then the predator extinction point S 1 = r β 1 , 0 of system (6) is globally asymptotically stable.
Proof. Define a positive definite Lyapunov function at S 1 .
According to Lemma 3.1. in [54], the α-order derivative of V 3 (x, y) along the solution of system (6) satisfies the following.
Based on Theorem 2, we have y(t) ≤ γ 2 and, thus, we have the following.
We note that if m ≥ br , 0). Hence, the only invariant set on which D α * V 3 (x, y) = 0 is the singleton {S 1 }. By Lemma 4.6 in [55], the predator extinction point S 1 is globally asymptotically stable.

Numerical Simulations
In this section, some numerical simulations of system (6) are presented to verify the analytical results such as the stability of equilibrium points and a Hopf bifurcation and to explore the dynamical behavior of the system (6) for both weak and strong Allee effects, respectively. The recent development of the numerical schemes for the fractionalorder system such as the Grünwald-Letnikov method [31,56], the predictor-corrector method [57][58][59], the homotopy perturbation method [60,61], and the Laplace adomian decomposition method [62] have been used to solve the fractional-order differential equations. Here, we apply the fractional-order predictor-corrector method provided by Diethelm et al. [63] to obtain numerical solutions of the system (6). Based on that, we divide this section into three subsections that demonstrate the effects of the capturing rate to the predation rate (b), the Allee threshold (m), and the order of fractional system (α) on the stability of equilibrium points. Since parameter values based on real-life observations are not available, we use hypothetical parameters that correspond to the analytical results.

The Influence of the Capturing Rate (b)
To understand how the capturing rate (b) could influence the dynamics of system (6) in the weak Allee effect case (m < 0), we use the following hypothetical parameter values. r = 0.5, β 1 = 0.05, c = 1, q = 0.1, s = 0.1, l = 1, β 2 = 0.05, m = −1, and n = 3. (26) The bifurcation diagrams for both predator and prey populations controlled by b ∈ [0.2, 0.6] and α = 0.98 are depicted in Figure 1a. When b < b * w1 ≈ 0.24103, the interior pointŴ is asymptotically stable. For example, we plot a phase portrait of the system (6) in Figure 1b for b = 0.20. In this case, one can easily compute that χ 2 1w − 4χ 2w = −0.06299 < 0 and χ 1w = −0.08727 < 0. According to Theorem 4, the interior point W = (4.30171, 8.80734) is asymptotically stable. The stability of the interior pointŴ is also achieved in the interval b ∈ (b * w2 = 0.53846, b * w3 = 0.55487). To show such typical behavior, we plot a phase portrait of the system (6) in Figure 1d for b = 0.54. Furthermore, our numerical simulations also show the existence of limit-cycles solution enclosing the interior pointŴ for b ∈ (b * w1 , b * w2 ) as provided in the green area in Figure 1a. As an example, we plot some numerical solutions of the system (6) with b = 0.27 in Figure 1c. It is observed that all solutions of system (6) converge to a limit-cycle and the interior pointŴ = (2.76631, 5.82563) losses its stability. This situation shows that the system (6) undergoes a Hopf bifurcation with respect to b. When b > b * w3 , the interior pointŴ does not exist anymore and the prey extinction point W 2 = (0, 1) becomes stable via forward bifurcation. The numerical solutions of system (6), which show the stability of W 2 = (0, 1), are shown in Figure 1e.
Next, we numerically investigate the effect of the capturing rate (b) under the strong Allee effect case (m > 0) using the following hypothetical parameters. Here, we also take α = 0.98. We notice that the predator extinction point S 1 = (5, 0) is always asymptotically stable, see Theorem 5b. Then, by varying b from 0.2 to 0.6, we plot the bifurcation diagrams for both predator and prey population in Figure 2a. It is found that the system (6) in this case has similar qualitative behavior as in the previous simulation (see Figures 1 and 2), with the exception of the stability properties of the predator extinction point S 1 = (5, 0). In the first simulation, the predator extinction point is always unstable while, in the latter case, it is always locally asymptotically stable. The local stability properties of the predator extinction point can be clearly observed from the phase portraits depicted in Figure 2b-e. Figures 1 and 2, the capturing rate (b) has great influence to both prey and predator population densities. For parameter values as in (26) and (27), for both weak and strong Allee effects, the densities of both prey and predator populations decrease with increasing capturing rates. In particular, increasing the capturing rate such that b > b * w3 (for the weak Allee effect case) or b > b * s3 (for the strong Allee effect case) may stabilize the prey extinction point (0, 1) as stated in Theorems 3c and 5c. This means that the prey population will become extinct while the predator still survives. We also notice that for the case of strong Allee effect (m > 0), the system (6) may exhibit a bistability phenomenon, see Figure 2b,d,e. Hence, the system (6) is highly sensitive to the initial condition when the Allee effect is strong. From the ecological point of view, this result is quite interesting. We can still maintain the existence of both prey and predator even under conditions of a strong Allee effect, that is, as long as we have a sufficiently large initial predator density. Furthermore, Figure 1b,d,e confirm the global asymptotic stability behavior of system (6) for the weak Allee effect case.

The Impacts of the Allee Threshold (m)
We next study numerically the impacts of the Allee threshold (m) to the dynamics of system (6) using the following hypothetic values of parameters. r = 0.5, β 1 = 0.1, b = 0.4, c = 1, q = 0.1, s = 0.1, l = 1, β 2 = 0.05, and n = 3. (28) In Figure 3a we show the bifurcation diagrams of both predator and prey populations for parameter m ∈ [−2, 1] and α = 0.98. It is shown that when m < m * 1 ≈ −1.73333, the prey extinction point W 2 is asymptotically stable. If we increase m such that m > m * 1 − 1.73333, then the prey extinction point W 2 loses its stability when an asymptotically stable interior pointŴ appears. This shows that the system (6) exhibits a forward bifurcation. The interior pointŴ remains asymptotically stable for m ∈ (m * 1 , m * 2 ). Moreover, Figure 3a shows the occurence of a Hopf bifurcation when m passes through m * 2 ≈ −1.50769. Indeed, in the interval m ∈ (m * 2 , m * 3 = 0), there is no stable equilibrium point and there exists a limit cycle around the interior pointŴ; see the green area in Figure 3a. Our further observation shows that the system (6) with m ∈ (m * 3 , m * 4 ≈ 0.82356) has a bistability phenomenon where both interior pointŜ and the predator extinction point S 1 = (5, 0) are locally asymptotically stable. If we increase m such that m > m * 4 , the interior pointŜ loses its existence and the predator extinction point S 1 becomes a unique stable equilibrium point of the system (6).
In order to provide a better view of the dynamics of the system (6) with parameter values (28) and α = 0.98, we ploted several phase-portraits for different values of m in Figure 3b-e. Figure 3b shows the phase-portrait of the system (6) with m = −1.79 < m * 1 . Here, we have 0.5 = r < bȳ w 1+qȳ w = 0.50960 and 1.21 = m + n < β 2 ls (ȳ w + n) 2 = 9.94580. According to Theorem 3c, the prey extintion point W 2 = (0, 1.46) is asymptotically stable. Figure 3b confirms this behavior where all numerical solutions converge to the prey extinction pointW 2 = (0, 1.46). The phase portrait of the system (6) with m = −0.46 > m * 2 is depicted in Figure 3c. In this respect, we have χ 1w = 0.03488 > 0 and 0.87258 = α * < α = 0.98. Thus, the stability condition of the interior pointŴ = (1.21751, 2.31593) in Figure 3c, where all numerical solutions of the system (6) are convergent to a limit-cycle. The appearence of a stable limit-cycle indicates the presence of a Hopf bifurcation in the system. Next, we consider a case of m ∈ (m * 3 , m * 4 ) by choosing m = 0.4. The phase-portrait in Figure 3d shows that numerical solutions are convergent to the predator extinction point S 1 = (5, 0) or to the interior pointŜ 1 = (2.24039, 2.40121), depending on the initial conditions. Hence, the system (6) exhibits a bistability phenomenon. Finally, if we take m = 0.90 > m * 4 , then we have situation where the interior pointŜ disappears and the predator extinction point S 1 = (5, 0) becomes the only equilibrium point which is asymptotically stable. This behavior is plotted in Figure 3e.
Note that, when the Allee effect is weak, i.e., when −n < m < 0, the predator always exists as depicted in the interval [−2, m * 3 ] in Figure 3a. In this case, the predator growth rate is always positive and the prey population could suffer from extinction as the m value decreases. On the other hand, when the Allee effect is strong, i.e., when m > 0, there is a condition where the predator is always extinct as depicted by the blue solid line in the interval [m * 3 , 1] in Figure 3a. However, generally speaking, the system (6) could have a positive or negative growth rate on the predator population since there is a bistability phenomenon in that interval as explained before. These situations show that our system (6) may provide the condition for the existence of the predator population which is affected by the double Allee effect.

The Effects of the Order of Fractional System (α)
In the following simulations, we study the influence of the order of the fractional derivative (α) by considering the system (6) with the following hypothetical parameter values.
Based on parameter values (29) and Theorem 7, we can find a critical value α * ≈ 0.80631 such that the interior point is asymptotically stable if α < α * and it is unstable if α > α * . In order to observe such behavior, we plot the bifurcation diagram for α ∈ [0.7, 1.0], see Figure 4a. It is observed that the interior pointŜ = (0.38825, 1.80915) is asymptotically stable when α < α * . If α > α * , then the interior pointŜ loses its stability and all numerical solutions converge to a limit-cycle via a Hopf bifurcation. The Hopf bifurcation can also be observed from the phase-portraits shown in Figure 4b for α = 0.73 < α * and Figure 4c for α = 0.89 > α * .
In Figure 4a, we observe that the order of the fractional derivative does not affect the stability of the interior point as long as 0 < α < α * . To verify this property, we plot in Figure 5 the time series of both prey and predator populations for α = 0.5, 0.6, 0.7, 0.75. It is shown that all solutions are indeed convergent to the interior pointŜ = (0.38825, 1.80915), but the solution of system (6) with a higher-order fractional derivative has faster convergence.

Conclusions
In this paper, a fractional-order Leslie-Gower predator-prey model with Beddington-DeAngelis functional response and double Allee effect in the predator population is proposed and the dynamics of the model has been analyzed. First, the existence, uniqueness, non-negativity, and boundedness of the solution have been proven. We then determined all possible non-negative equilibrium points and their local and global stability properties. We found that the model has four types of biologically feasible equilibrium. The extinction of both prey and predator point is always unstable for both weak and strong Allee effect cases. The predator extinction point is always stable for the strong Allee effect case, but it is always unstable for the weak Allee effect case. The prey extinction point and the interior point are conditionally stable. Our numerical simulations showed that for the case of the weak Allee effect, there is a capturing rate threshold b * such that for b > b * the prey population is extinct while the predator population still survives. However, for the case of the strong Allee effect, the situation is also dependent on the initial value. Here, if the initial predator population is relatively low then the predator will become extinct but the prey will survive. Additionally, we also proved the existence of a Hopf bifurcation about the interior point driven by the order of the fractional derivative (α) and the critical α * of this bifurcation has been determined analytically. The occurrence of the Hopf bifurcation has been confirmed by our numerical simulations. The existence of Hopf bifurcation controlled by α has also been observed in [31,32]. This shows that the unstable interior point of the first order system (i.e., the system is convergent to a limit cycle) may become stable in the fractional order system (i.e., the system converges to the interior point).