Next Article in Journal
A Numerical Study of Nonlinear Fractional Order Partial Integro-Differential Equation with a Weakly Singular Kernel
Next Article in Special Issue
State of Charge Estimation of Lithium-Ion Batteries Based on Fuzzy Fractional-Order Unscented Kalman Filter
Previous Article in Journal
Fractional Differential Equation-Based Instantaneous Frequency Estimation for Signal Reconstruction
Previous Article in Special Issue
Fractional Integral Inequalities for Exponentially Nonconvex Functions and Their Applications
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Department of Mathematics, Faculty of Mathematics and Natural Sciences, University of Brawijaya, Malang 65145, Indonesia
2
Department of Mathematics, Faculty of Mathematics and Natural Sciences, State University of Gorontalo, Bone Bolango 96119, Indonesia
*
Author to whom correspondence should be addressed.
Fractal Fract. 2021, 5(3), 84; https://doi.org/10.3390/fractalfract5030084
Submission received: 29 June 2021 / Revised: 24 July 2021 / Accepted: 28 July 2021 / Published: 1 August 2021
(This article belongs to the Special Issue Fractional Order Systems and Their Applications)

Abstract

:
In this paper, a modified 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 sufficient 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 confirm the theoretical results and can explore more dynamical behaviors of the system, such as the bi-stability and forward bifurcation.

1. 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.

2. 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]:
d x d t = f ( x ) x g ( x , y ) y , d y d t = k g ( x , y ) y μ y ,
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.
d x d t = r ^ x β ^ 1 x 2 b ^ x y 1 + c ^ x + q ^ y , d y d t = s ^ y σ ^ y 2 .
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 parameters r ^ , β ^ 1 , b ^ , c ^ , q ^ , s ^ , σ ^ 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; s ^ and σ ^ are the intrinsic growth rate of the predator and the predator intraspecific competition coefficient, respectively; b ^ and c ^ 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 s ^ . 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:
d x d t = r ^ x β ^ 1 x 2 b ^ x y 1 + c ^ x + q ^ y , d y d t = s ^ y m ^ y + n ^ y β ^ 2 l ^ + x y 2 ,
where m ^ is the Allee threshold, n ^ is the auxiliary Allee effect constant with n ^ > 0 and n ^ > m ^ . In the second equation of model (3), the intrinsic growth rate of the predator s ^ is affected by double Allee effects. Without the intraspecific competition for the predator, the per capita growth rate of the predator is reduced from s ^ to s ^ y m ^ y + n ^ due to the Allee effect [45]. Therefore, the Allee effect is strong if m ^ > 0 and weak if n ^ < m ^ 0 [9,18].
In order to seize the entire time population growth condition, we consider a fractional-order derivative to the left-hand side of the classical derivative system (3) as follows:
D * α x = r ^ x β ^ 1 x 2 b ^ x y 1 + c ^ x + q ^ y , D * α y = s ^ y m ^ y + n ^ y β ^ 2 l ^ + x y 2 ,
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:
D * α f ( t ) = 1 Γ ( n α ) 0 t f n ( τ ) ( t τ ) n α 1 d τ ,
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 ( t i m e 1 ). Thus, we have a new system.
D * α x = r ^ α x β ^ 1 α x 2 b ^ α x y 1 + c ^ x + q ^ y , D * α y = s ^ α y m ^ y + n ^ y β ^ 2 α l ^ + x y 2 .
For simplification, we replace system (5) with the redefined parameters as follows:
D * α x = r x β 1 x 2 b x y 1 + c x + q y , D * α y = s y m y + n y β 2 l + x y 2 ,
where
r ^ α = r , β ^ 1 α = β 1 , b ^ α = b , c ^ = c , q ^ = q , s ^ α = s , m ^ = m , n ^ = n , β ^ 2 α = β 2 , l ^ = l .
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.

3. 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.

3.1. Existence and Uniqueness

The existence and uniqueness of the solution of system (6) are examined in the region Ω × [ 0 , T ] , where Ω = { ( x , y ) R 2 : max { | x | , | y | } M } and T < + . Let a mapping F ( Z ) = ( F 1 ( Z ) , F 2 ( Z ) ) with the following.
F 1 ( Z ) = r x β 1 x 2 b x y 1 + c x + q y , F 2 ( Z ) = s y 2 y + n m s y y + n β 2 y 2 l + x .
For any Z = ( x , y ) , Z ¯ = ( x ¯ , y ¯ ) , Z , Z ¯ Ω , it follows from (7) that the following is the case:
| | F ( Z ) F ( Z ¯ ) | | = | F 1 ( Z ) F 1 ( Z ¯ ) | + | F 2 ( Z ) F 2 ( Z ¯ ) | = r ( x x ¯ ) β 1 ( x 2 x ¯ 2 ) b x y 1 + c x + q y x ¯ y ¯ 1 + c x ¯ + q y ¯ + s y 2 y + n y ¯ 2 y ¯ + n m s y y + n y ¯ y ¯ + n β 2 y 2 l + x y ¯ 2 l + x ¯ r | x x ¯ | + β 1 | x 2 x ¯ 2 | + b x y ( 1 + c x ¯ + q y ¯ ) x ¯ y ¯ ( 1 + c x + q y ) ( 1 + c x + q y ) ( 1 + c x ¯ + q y ¯ ) + s y 2 ( y ¯ + n ) y ¯ 2 ( y + n ) ( y + n ) ( y ¯ + n ) + | m | s y ( y ¯ + n ) y ¯ ( y + n ) ( y + n ) ( y ¯ + n ) + β 2 y 2 ( l + x ¯ ) y ¯ 2 ( l + x ) ( l + x ) ( l + x ¯ ) = r | x x ¯ | + β 1 | x + x ¯ | | x x ¯ | + b ( x + c x x ¯ ) ( y y ¯ ) + ( y ¯ + q y y ¯ ) ( x x ¯ ) ( 1 + c x + q y ) ( 1 + c x ¯ + q y ¯ ) + s ( y y ¯ + n ( y + y ¯ ) ) ( y y ¯ ) ( y + n ) ( y ¯ + n ) + | m | s n ( y y ¯ ) ( y + n ) ( y ¯ + n ) + β 2 ( l + x ) ( y + y ¯ ) ( y y ¯ ) y 2 ( x x ¯ ) ( l + x ) ( l + x ¯ ) r | x x ¯ | + β 1 | x + x ¯ | | x x ¯ | + b | ( x + c x x ¯ ) ( y y ¯ ) | + b | ( y ¯ + q y y ¯ ) ( x x ¯ ) | + s n 2 | ( y y ¯ + n ( y + y ¯ ) ) ( y y ¯ ) | + | m | s n | y y ¯ | + β 2 l 2 | ( l + x ) ( y + y ¯ ) ( y y ¯ ) y 2 ( x x ¯ ) | r | x x ¯ | + 2 β 1 M | x x ¯ | + ( b M + b c M 2 ) | y y ¯ | + ( b M + b q M 2 ) | x x ¯ | + s M 2 n 2 + 2 s M n | y y ¯ | + | m | s n | y y ¯ | + 2 β 2 M l + 2 β 2 M 2 l 2 | y y ¯ | + β 2 M 2 l 2 | x x ¯ | = L 1 | x x ¯ | + L 2 | y y ¯ | L | | Z Z ¯ | |
where
L 1 = r + ( 2 β 1 + b ) M + b q + β 2 l 2 M 2 , L 2 = | m | s n + b + 2 s n + 2 β 2 l M + b c + s n 2 + 2 β 2 l 2 M 2 , L = max L 1 , L 2 .
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 .

3.2. Non-Negativity and Boundedness

In order to prove that all solutions of system (6) are non-negative and bounded, let R + 2 = { W = ( x , y ) T R 2 | x ( t ) 0 , y ( t ) 0 } be the non-negative quadrant on the x y 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.
Proof. 
Let W ( t 0 ) = x ( t 0 ) y ( t 0 ) R + 2 and assume that W ( t ) = x ( t ) y ( t ) for t t 0 be the solutions of system (6).
Suppose that assumption is false, then there exists t * > t 0 such that W ( t ) > 0 for t 0 t < t * , W ( t * ) = 0 , and W ( t + * ) < 0 for t + * > t * . From system (6), one has the following.
D * α W ( t ) | t = t * = 0 .
Based on Lemma 1 in [51], we have W ( t + * ) = 0 , which contradicts with W ( t + * ) < 0 for t + * > t * . Therefore, we conclude W ( t ) 0 for all t 0 .
Next, we prove the boundedness of all solutions of system (6). From the first equation of system (6), we have the following.
D * α x ( t ) + x ( t ) = r x β 1 x 2 b x y 1 + c x + q y + x = β 1 x 2 + ( 1 + r ) x b x y 1 + c x + q y = β 1 x ( 1 + r ) 2 β 1 2 + ( 1 + r ) 2 4 β 1 b x y 1 + c x + q y ( 1 + r ) 2 4 β 1 .
By Lemma 3 in [51], we have the following:
x ( t ) x ( t 0 ) ( 1 + r ) 2 4 β 1 E α ( t t 0 ) α + ( 1 + r ) 2 4 β 1 ( 1 + r ) 2 4 β 1 , t ,
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.
Ω 1 = x ( t ) ( 1 + r ) 2 4 β 1 + ϵ 1 = γ 1 , ϵ 1 > 0 .
From the second equation of system (6), we have the following.
D * α y ( t ) + s y ( t ) = s y 2 y + n m s y y + n β 2 y 2 l + x + s y .
We also have x ( t ) γ 1 from (10), then the following obtains.
D * α y ( t ) + s y ( t ) s y 2 y + n m s y y + n β 2 y 2 l + γ 1 + s y s y 2 y m s y y + n β 2 y 2 l + γ 1 + s y = β 2 y 2 l + γ 1 + 2 s y m s y y + n = β 2 l + γ 1 y ( l + γ 1 ) s β 2 2 + ( l + γ 1 ) s 2 β 2 m s y y + n ( l + γ 1 ) s 2 β 2
Again by using Lemma 3 in [51], for the strong Allee effect ( m > 0 ) , we have the following.
y ( t ) y ( t 0 ) ( l + γ 1 ) s β 2 E α s ( t t 0 ) α + ( l + γ 1 ) s β 2 ( l + γ 1 ) s β 2 , t .
However, for the weak Allee effect ( m < 0 ) , we have the following.
y ( t ) y ( t 0 ) ( l + γ 1 ) s β 2 m E α s ( t t 0 ) α + ( l + γ 1 ) s β 2 m ( l + γ 1 ) s β 2 m , t .
Therefore, the solution of y ( t ) with initial condition y ( t 0 ) are confined to region Ω 2 where
Ω 2 = y ( t ) γ 2 ,
and where the following is the case.
γ 2 = ( l + γ 1 ) s β 2 + ϵ 2 , ϵ 2 > 0 , m > 0 ( l + γ 1 ) s β 2 m + ϵ 2 , ϵ 2 > 0 , m < 0

4. 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.
(1)
The equilibrium points of system (6) for the weak Allee effect ( m < 0 ) are as follows:
(a)
Both prey and predator extinction point W 0 = ( 0 , 0 ) , which always exists;
(b)
The predator extinction point W 1 = ( r β 1 , 0 ) , which always exists;
(c)
The prey extinction point W 2 = 0 , y ¯ w where we have the following.
y ¯ w = l s β 2 n 2 + ( l s β 2 n ) 2 4 m l s β 2 2
Denote W 2 as always existing.
(d)
The interior point W ^ = ( x ^ w , y ^ w ) where x ^ w = β 2 y ^ w ( y ^ w + n ) s ( y ^ w m ) l and y ^ w are all positive roots of the quartic equation (14):
a 1 y 4 + a 2 y 3 + a 3 y 2 + a 4 y + a 5 = 0 ,
where
a 1 = β 1 β 2 q s + c β 1 β 2 2 , a 2 = ( b q r l q β 1 ) s 2 + ( β 1 β 2 ( 1 + n q m q 2 c l ) c r β 2 ) s + 2 c β 1 β 2 2 n , a 3 = ( c β 1 l 2 + 2 l m q β 1 + l c r + 2 m q r 2 b m l β 1 r ) s 2 + ( β 1 β 2 ( 2 c l m 2 c l n n m q m + n ) + c r β 2 ( m n ) ) s + c β 1 β 2 2 n 2 , a 4 = ( 2 c l 2 m β 1 l m 2 q β 1 2 c l m r m 2 q r + b m 2 + 2 l m β 1 + 2 m r ) s 2 + ( 2 c l m n β 1 β 2 + c m n r β 2 m n β 1 β 2 ) s , a 5 = m 2 ( c l 1 ) ( l β 1 + r ) s 2 .
(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 , y ¯ s where y ¯ s is the positive solution of the quadratic equation y 2 l s β 2 n y + m l s β 2 = 0 . The existence of S 2 , 3 is described as follows:
(i)
If l s β 2 n 2 < 4 m l s β 2 , then the prey extinction point does not exist.
(ii)
If l s β 2 n 2 = 4 m l s β 2 and n < l s β 2 , then there exists a unique prey extinction point, S 2 = S 3 = 0 , 1 2 l s β 2 n ) .
(iii)
If l s β 2 n 2 > 4 m l s β 2 , then there exist two prey extinction points, i.e., the following is the case.
S 2 , 3 = 0 , l s β 2 n 2 ± ( l s β 2 n ) 2 4 m l s β 2 2
(d)
The interior point S ^ = ( x ^ s , y ^ s ) exists if y ^ s > m where x ^ s = β 2 y ^ s ( y ^ s + n ) s ( y ^ s m ) l and y ^ 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.
J ( x * , y * ) = J 1 , 1 J 1 , 2 J 2 , 1 J 2 , 2
J 1 , 1 = r 2 β 1 x * b y * 1 + c x * + q y * + b c x * y * ( 1 + c x * + q y * ) 2 J 1 , 2 = b x * 1 + c x * + q y * + b q x * y * ( 1 + c x * + q y * ) 2 J 2 , 1 = β 2 y * 2 ( l + x * ) 2 J 2 , 2 = s ( y * m ) y * + n + s y * 1 y * + n y * m ( y * + n ) 2 2 β 2 y * l + x * .
Theorem 3.
The stability properties of trivial and axial equilibrium points of system (6) for the weak Allee effect ( m < 0 ) are as follows:
(a) 
W 0 = ( 0 , 0 ) is always unstable.
(b) 
W 1 = ( r β 1 , 0 ) is always a saddle point.
(c) 
W 2 = ( 0 , y ¯ w ) is locally asymptotically stable if r < b y ¯ w 1 + q y ¯ w and m + n < β 2 l s ( y ¯ w + n ) 2 .
Proof. 
(a)
By substituting W 0 to (15), we obtain the Jacobian matrix.
J ( W 0 ) = r 0 0 m s n .
Therefore, the eigenvalues of (16) are λ 1 = r > 0 and λ 2 = m s n > 0 , since | arg ( λ 1 ) | = | arg ( λ 2 ) | = 0 < α π 2 , E 0 is always unstable by the Matignon condition [25].
(b)
By substituting W 1 to (15), we obtain the Jacobian matrix.
J ( W 1 ) = r b r c r + β 1 0 m s n .
The Jacobian matrix (17) has eigenvalues λ 1 = r < 0 and λ 2 = m s n > 0 , showing that | arg ( λ 1 ) | = π > α π 2 and | arg ( λ 2 ) | = 0 < α π 2 . Hence, W 1 is always a saddle point.
(c)
By evaluating (15) at W 2 , we obtain the following.
J ( W 2 ) = r b y ¯ w 1 + q y ¯ w 0 β 2 y ¯ w 2 l 2 s ( m + n ) ( y ¯ w + n ) 2 β 2 l y ¯ w .
The eigenvalues of (18) are as follows.
λ 1 = r b y ¯ w 1 + q y ¯ w and λ 2 = s ( m + n ) ( y ¯ w + n ) 2 β 2 l y ¯ w .
Thus, | arg ( λ 1 ) | = π > α π 2 and | arg ( λ 2 ) | = π > α π 2 , whenever r < b y ¯ w 1 + q y ¯ w and m + n < β 2 l s ( y ¯ w + n ) 2 .
Theorem 4.
Suppose that m < 0 (weak Allee effect) and the following is the case.
χ 1 w = β 1 x ^ w + β 2 y ^ w l + x ^ w + b c x ^ w y ^ w ( 1 + c x ^ w + q y ^ w ) 2 + s ( m + n ) y ^ w ( y ^ w + n ) 2 χ 2 w = β 1 x ^ w y ^ w β 2 l + x ^ w s ( m + n ) ( y ^ w + n ) 2 + b s x ^ w ( 1 + c x ^ w + q y ^ w ) 2 β 2 c ( m + n ) y ^ w 2 + s ( 1 + c x ^ w ) ( y ^ w m ) 2 β 2 ( y ^ w + n ) 2 β 2 c y ^ w 2 s ( l + x ^ w ) α * = 2 π tan 1 4 χ 2 w ( χ 1 w ) 2 χ 1 w .
The interior point W ^ = ( x ^ w , y ^ w ) is locally asymptotically stable if the following is the case:
(i) 
χ 1 w 2 4 χ 2 w , χ 1 w < 0 , and χ 2 w > 0 .
(ii) 
χ 1 w 2 < 4 χ 2 w , and if χ 1 w < 0 , or χ 1 w > 0 and α < α * .
Proof. 
By evaluating (15) at the interior equilibrium point W ^ = ( x ^ w , y ^ w ) , we obtain the following.
J ( W ^ ) = β 1 x ^ w + b c x ^ w y ^ w ( 1 + c x ^ w + q y ^ w ) 2 b x ^ w ( 1 + c x ^ w ) ( 1 + c x ^ w + q y ^ w ) 2 s 2 ( y ^ w m ) 2 β 2 ( y ^ w + n ) 2 s ( m + n ) y ^ w ( y ^ w + n ) 2 β 2 y ^ w l + x ^ w
The Jacobian matrix (19) has polynomial characteristic λ 2 χ 1 w λ + χ 2 w = 0 . By utilizing the Routh–Hurwitz criterion for Caputo fractional-order [52], it follows that W ^ is locally asymptotically stable if condition ( i ) or ( i i ) is satisfied. □
Theorem 5.
The stability properties of trivial and axial equilibrium points of system (6) for strong Allee effect ( m > 0 ) are as follows:
(a) 
S 0 = ( 0 , 0 ) is a saddle point.
(b) 
S 1 = ( r β 1 , 0 ) is always locally asymptotically stable.
(c) 
S 2 , 3 = ( 0 , y ¯ s ) is locally asymptotically stable if r < b y ¯ s 1 + q y ¯ s and m + n < β 2 l s ( y ¯ s + n ) 2 .
Proof. 
(a)
By substituting S 0 to (15), we obtain the following.
J ( S 0 ) = r 0 0 m s n .
It is clear that the eigenvalues of (20) are λ 1 = r > 0 and λ 2 = m s n < 0 , and | arg ( λ 1 ) | = 0 < α π 2 and | arg ( λ 2 ) | = π > α π 2 . Thus, S 0 is a saddle point.
(b)
The Jacobian matrix (15) evaluated at S 1 is the following:
J ( S 1 ) = r b r c r + β 1 0 m s n ,
where its eigenvalues are λ 1 = r < 0 and λ 2 = m s n < 0 , since | arg ( λ 1 , 2 ) | = π > α π 2 , S 1 is always locally asymptotically stable.
(c)
By evaluating (15) at S 2 , 3 , we acquire the following.
J ( S 2 , 3 ) = r b y ¯ s 1 + q y ¯ s 0 β 2 y ¯ s 2 l 2 s ( m + n ) ( y ¯ s + n ) 2 β 2 l y ¯ s .
The eigenvalues of (22) are as follows.
λ 1 = r b y ¯ s 1 + q y ¯ s and λ 2 = s ( m + n ) ( y ¯ s + n ) 2 β 2 l y ¯ s .
Therefore S 2 , 3 is locally asymptotically stable if r < b y ¯ s 1 + q y ¯ s and m + n < β 2 l s ( y ¯ w + n ) 2 because in this case | arg ( λ 1 ) | = | arg ( λ 2 ) | = π > α π 2 .
Theorem 6.
For the case of strong Allee effect ( m > 0 ) , the interior point S ^ = ( x ^ s , y ^ s ) is locally asymptotically stable if the following is the case:
(i) 
χ 1 s 2 4 χ 2 s , χ 1 s < 0 , and χ 2 s > 0 .
(ii) 
χ 1 s 2 < 4 χ 2 s , and if χ 1 s < 0 , or χ 1 s > 0 and α < α * ;
where the following is the case.
χ 1 s = β 1 x ^ s + β 2 y ^ s l + x ^ s + b c x ^ s y ^ s ( 1 + c x ^ s + q y ^ s ) 2 + s ( m + n ) y ^ s ( y ^ s + n ) 2 χ 2 s = β 1 x ^ s y ^ s β 2 l + x ^ s s ( m + n ) ( y ^ s + n ) 2 + b s x ^ s ( 1 + c x ^ s + q y ^ s ) 2 β 2 c ( m + n ) y ^ s 2 + s ( 1 + c x ^ s ) ( y ^ s m ) 2 β 2 ( y ^ s + n ) 2 β 2 c y ^ s 2 s ( l + x ^ s ) α * = 2 π tan 1 4 χ 2 s ( χ 1 s ) 2 χ 1 s .
Proof. 
The Jacobian matrix (15) at interior equilibrium point S ^ = ( x ^ s , y ^ s ) is given by the following.
J ( S ^ ) = β 1 x ^ s + b c x ^ s y ^ s ( 1 + c x ^ s + q y ^ s ) 2 b x ^ s ( 1 + c x ^ s ) ( 1 + c x ^ s + q y ^ s ) 2 s 2 ( y ^ s m ) 2 β 2 ( y ^ s + n ) 2 s ( m + n ) y ^ s ( y ^ s + n ) 2 β 2 y ^ s l + x ^ s .
The Jacobian matrix (24) has polynomial characteristic λ 2 χ 1 s λ + χ 2 s = 0 . By using the Routh–Hurwitz criterion for Caputo fractional-order [52], the stability condition is completely proven. □
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.
Theorem 7
(Existence of Hopf bifurcation [53]). Let χ 1 w 2 < 4 χ 2 w ( o r χ 1 s 2 < 4 χ 2 s ) and χ 1 w > 0 ( o r χ 1 s > 0 ) . System (6) undergoes a Hopf bifurcation around the interior point W ^ ( o r S ^ ) when α crosses α * .
Proof. 
Based on Theorem 6, when χ 1 s 2 < 4 χ 2 s and χ 1 s > 0 , the eigenvalues of system (6) at S ^ are a pair of complex conjugate numbers with the real parts are positive. We also confirm that ϕ 1 , 2 ( α * ) = 0 and d ϕ ( α ) d α | α = α * 0 where ϕ i ( α ) = α π 2 min 1 i 2 arg ( λ i ( α ) ) . Based on Theorem 3 in [53], the equilibrium point S ^ undergoes a Hopf bifurcation when α crosses α * . The similar proof works for the weak Allee effect case. □

5. Global Stability

5.1. 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 , y ¯ w ) and the interior point W ^ = ( x ^ w , y ^ w ) are conditionally locally asymptotically stable. In the following, we study the global asymptotic stability of those equilibrium points.
Theorem 8.
If n y ¯ w γ 2 + n 1 m γ 2 + n s y ¯ w r 2 4 β 1 + β 2 y ¯ w 2 l , then W 2 = ( 0 , y ¯ w ) is globally asymptotically stable.
Proof. 
We consider the following positive definite Lyapunov function.
V 1 ( x , y ) = x + y y ¯ w y ¯ w ln y y ¯ w
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.
D * α V 1 ( x , y ) D * α x + y y ¯ w y D * α y = r x β 1 x 2 b x y 1 + c x + q y + ( y y ¯ w ) s ( y m ) y + n β 2 y l + x = β 1 x r 2 β 1 2 + r 2 4 β 1 b x y 1 + c x + q y + s y 2 y + n m s y y + n s y ¯ w y y + n + m s y ¯ w y + n β 2 y ¯ w y l + x + β 2 y ¯ w 2 l + x β 2 l + x ( y y ¯ w ) 2 r 2 4 β 1 + s y 2 y + n m s y y + n s y ¯ w y y + n + m s y ¯ w y + n + β 2 y ¯ w 2 l + x β 2 l + x ( y y ¯ w ) 2 r 2 4 β 1 + s y m s y n s y ¯ w y γ 2 + n + m s y ¯ w γ 2 + n + β 2 y ¯ w 2 l β 2 l + x ( y y ¯ w ) 2 = r 2 4 β 1 + m s y ¯ w γ 2 + n + β 2 y ¯ w 2 l + s 1 m n y ¯ w γ 2 + n y β 2 l + x ( y y ¯ w ) 2
Since n y ¯ w γ 2 + n 1 m γ 2 + n s y ¯ w r 2 4 β 1 + β 2 y ¯ w 2 l , we obtain the following.
D * α V 1 β 2 l + x ( y y ¯ w ) 2 .
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 = y ¯ w . Substituting y = y ¯ w to the second equation of system (6) we obtain the following.
0 = D * α y ¯ w = s y ¯ w y ¯ w m y ¯ w + n β 2 y ¯ w 2 l + x .
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. □
Theorem 9.
The interior point W ^ of model (6) is globally asymptotically stable if the following is the case:
(i) 
n y ^ w γ 2 + n + β 2 y ^ w s ( l + γ 1 ) 1 m β 2 ( γ 2 + n ) y ^ w l s and;
(ii) 
b < min 2 β 1 ( 1 + c x ^ w + q y ^ w ) ( 1 + c γ 1 + q γ 2 ) 1 + c x ^ w + 2 c y ^ w ( 1 + c γ 1 + q γ 2 ) , 2 β 1 β 2 r ( 1 + c γ 1 + q γ 2 ) ( 1 + c x ^ w ) ( l + γ 1 ) .
Proof. 
Consider the following positive definite Lyapunov function.
V 2 ( x , y ) = r β 1 ( 1 + c x ^ w + q y ^ w ) x x ^ w x ^ w ln x x ^ w + y y ^ w y ^ w ln y y ^ w .
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.
D * α V 2 r β 1 ( 1 + c x ^ w + q y ^ w ) x x ^ w x D * α x + y y ^ w y D * α y = r β 1 ( 1 + c x ^ w + q y ^ w ) ( x x ^ w ) r β 1 x b y 1 + c x + q y + ( y y ^ w ) s ( y m ) y + n β 2 y l + x = r ( 1 + c x ^ w + q y ^ w ) ( x x ^ w ) 2 b r ( 1 + c x ^ w ) ( x x ^ w ) ( y y ^ w ) β 1 ( 1 + c x + q y ) + b c r y ^ w ( x x ^ w ) 2 β 1 ( 1 + c x + q y ) + s y 2 y + n m s y y + n s y ^ w y y + n + m s y ^ w y + n β 2 y ^ w y l + x + β 2 y ^ w 2 l + x β 2 ( y y ^ w ) 2 l + x b c r β 1 y ^ w r ( 1 + c x ^ w + q y ^ w ) ( x x ^ w ) 2 + b r β 1 ( 1 + c x ^ w ) ( 1 + c γ 1 + q γ 2 ) ( x x ^ w ) 2 + ( y y ^ w ) 2 2 + s y m s y n s y ^ w y γ 2 + n + m s y ^ w γ 2 + n β 2 y ^ w y l + γ 1 + β 2 y ^ w 2 l β 2 ( y y ^ w ) 2 l + γ 1 = b c r β 1 y ^ w r ( 1 + c x ^ w + q y ^ w ) + b r ( 1 + c x ^ w ) 2 β 1 ( 1 + c γ 1 + q γ 2 ) ( x x ^ w ) 2 β 2 l + γ 1 b r ( 1 + c x ^ w ) 2 β 1 ( 1 + c γ 1 + q γ 2 ) ( y y ^ w ) 2 + s m s n s y ^ w γ 2 + n β 2 y ^ w l + γ 1 y + m s y ^ w γ 2 + n + β 2 y ^ w 2 l .
It is clear that if
n y ^ w γ 2 + n + β 2 y ^ w s ( l + γ 1 ) 1 m β 2 ( γ 2 + n ) y ^ w l s
and
b < min 2 β 1 ( 1 + c x ^ w + q y ^ w ) ( 1 + c γ 1 + q γ 2 ) ( 1 + c x ^ w + 2 c y ^ w ( 1 + c γ 1 + q γ 2 ) ) , 2 β 1 β 2 r ( 1 + c γ 1 + q γ 2 ) ( 1 + c x ^ w ) ( l + γ 1 ) ,
then D * α V 2 ( x , y ) 0 for all ( x , y ) R + 2 . Moreover, D * α V 2 ( x , y ) = 0 implies that ( x , y ) = ( x ^ , y ^ ) . Based on that, the only invariant set on which D * α V 2 ( x , y ) = 0 is the singleton { W ^ } . By Lemma 4.6 in [55], it follows that W ^ is globally asymptotically stable. □

5.2. 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 b r β 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 .
V 3 ( x , y ) = x r β 1 r β 1 ln β 1 x r + y .
According to Lemma 3.1. in [54], the α -order derivative of V 3 ( x , y ) along the solution of system (6) satisfies the following.
D * α V 3 ( x , y ) x r β 1 x D * α x + D * α y = x r β 1 x r β 1 x b y 1 + c x + q y x + y s ( y m ) y + n β 2 y l + x = β 1 x r β 1 2 b x y 1 + c x + q y + b r y β 1 ( 1 + c x + q y ) + s y 2 y + n m s y y + n β 2 y 2 l + x β 1 x r β 1 2 + b r y β 1 + s y m s y y + n .
Based on Theorem 2, we have y ( t ) γ 2 and, thus, we have the following.
D * α V 3 ( x , y ) β 1 x r β 1 2 + b r y β 1 + s y m s y γ 2 + n = β 1 x r β 1 2 + b r β 1 + s m s γ 2 + n y .
We note that if m b r β 1 s + 1 ( γ 2 + n ) , then D * α V 3 ( x , y ) 0 for all ( x , y ) R + 2 . Furthermore, D * α V 3 ( x , y ) = 0 implies that ( x , y ) = ( r β 1 , 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. □

6. 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 fractional-order 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.

6.1. 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 .
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 w 1 * 0.24103 , the interior point W ^ 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 χ 1 w 2 4 χ 2 w = 0.06299 < 0 and χ 1 w = 0.08727 < 0 . According to Theorem 4, the interior point W ^ = ( 4.30171 , 8.80734 ) is asymptotically stable. The stability of the interior point W ^ is also achieved in the interval b ( b w 2 * = 0.53846 , b w 3 * = 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 W ^ for b ( b w 1 * , b w 2 * ) 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 W ^ = ( 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 w 3 * , the interior point W ^ 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.
r = 0.5 , β 1 = 0.1 , c = 1 , q = 0.1 , s = 0.1 , l = 1 , β 2 = 0.05 , m = 0.4 , and n = 0.2 .
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 Figure 1 and Figure 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.
Remark 1.
Based on Figure 1 and Figure 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 w 3 * (for the weak Allee effect case) or b > b s 3 * (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.

6.2. 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 .
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 W ^ appears. This shows that the system (6) exhibits a forward bifurcation. The interior point W ^ 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 W ^ ; 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 S ^ 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 S ^ 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 y ¯ w 1 + q y ¯ w = 0.50960 and 1.21 = m + n < β 2 l s ( y ¯ 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 point W 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 χ 1 w = 0.03488 > 0 and 0.87258 = α * < α = 0.98 . Thus, the stability condition of the interior point W ^ = ( 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 S ^ 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 S ^ 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.

6.3. 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.
r = 0.5 , β 1 = 0.1 , b = 0.4 , c = 1 , q = 0.1 , s = 0.1 , l = 1 , β 2 = 0.05 , m = 0.5 , n = 0.2 .
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 S ^ = ( 0.38825 , 1.80915 ) is asymptotically stable when α < α * . If α > α * , then the interior point S ^ 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 S ^ = ( 0.38825 , 1.80915 ) , but the solution of system (6) with a higher–order fractional derivative has faster convergence.

7. 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).

Author Contributions

Conceptualization, I.D. and E.R.; methodology, I.D. and A.S.; software, E.R.; validation, I.D. and T.; formal analysis, E.R., I.D. and A.S.; investigation, E.R. and A.S.; resources, I.D.; data curation, E.R. and T.; writing—original draft preparation, E.R.; writing—review and editing, I.D., A.S. and T.; visualization, E.R. and A.S.; supervision, I.D., A.S. and T.; project administration, A.S.; funding acquisition, I.D. All authors have read and agreed to the published version of the manuscript.

Funding

This research is funded by FMIPA-UB via PNBP-University of Brawijaya according to DIPA-UB number DIPA-023.17.2.677512/2021 and the Professor and Doctoral Research Grant 2021 under contract number 1604/UN10.F09/PN/2021.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The authors thank the editors and the reviewers for their critical suggestions, contributions, and comments that greatly improved this manuscript.

Conflicts of Interest

The authors declare no conflict of interest in this paper.

References

  1. Cushing, J.M. Backward bifurcations and strong Allee effects in matrix models for the dynamics of structured populations. J. Biol. Dyn. 2014, 8, 57–73. [Google Scholar] [CrossRef] [Green Version]
  2. Buffoni, G.; Groppi, M.; Soresina, C. Dynamics of predator-prey models with a strong Allee effect on the prey and predator-dependent trophic functions. Nonlinear Anal. Real World Appl. 2016, 30, 143–169. [Google Scholar] [CrossRef]
  3. Sasmal, S.K.; Chattopadhyay, J. An eco-epidemiological system with infected prey and predator subject to the weak Allee effect. Math. Biosci. 2013, 246, 260–271. [Google Scholar] [CrossRef]
  4. Saifuddin, M.; Biswas, S.; Samanta, S.; Sarkar, S.; Chattopadhyay, J. Complex dynamics of an eco-epidemiological model with different competition coefficients and weak Allee in the predator. Chaos Solitons Fractals 2016, 91, 270–285. [Google Scholar] [CrossRef]
  5. Pal, S.; Sasmal, S.K.; Pal, N. Chaos control in a discrete-time predator-prey model with weak Allee effect. Int. J. Biomath. 2018, 11, 1–26. [Google Scholar] [CrossRef]
  6. Arancibia-Ibarra, C.; Flores, J.; Bode, M.; Pettet, G.; Van Heijster, P. A modified May–Holling–Tanner predator-prey model with multiple Allee effects on the prey and an alternative food source for the predator. Discret. Contin. Dyn. Syst. 2021, 26, 943–962. [Google Scholar] [CrossRef]
  7. Berec, L.; Angulo, E.; Courchamp, F. Multiple Allee effects and population management. Trends Ecol. Evol. 2007, 22, 185–191. [Google Scholar] [CrossRef]
  8. Pavlová, V.; Berec, L.; Boukal, D.S. Caught between two Allee effects: Trade-off between reproduction and predation risk. J. Theor. Biol. 2010, 264, 787–798. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  9. Feng, P.; Kang, Y. Dynamics of a modified Leslie–Gower model with double Allee effects. Nonlinear Dyn. 2015, 80, 1051–1062. [Google Scholar] [CrossRef]
  10. Lai, L.; Zhu, Z.; Chen, F. Stability and bifurcation in a predator–prey model with the additive Allee effect and the fear effect. Mathematics 2020, 8, 1280. [Google Scholar] [CrossRef]
  11. Martínez-Jeraldo, N.; Aguirre, P. Allee effect acting on the prey species in a Leslie–Gower predation model. Nonlinear Anal. Real World Appl. 2019, 45, 895–917. [Google Scholar] [CrossRef]
  12. Zhang, C.; Yang, W. Dynamic behaviors of a predator–prey model with weak additive Allee effect on prey. Nonlinear Anal. Real World Appl. 2020, 55, 103137. [Google Scholar] [CrossRef]
  13. Xiao, Z.; Xie, X.; Xue, Y. Stability and bifurcation in a Holling type II predator–prey model with Allee effect and time delay. Adv. Differ. Equ. 2018, 2018, 288. [Google Scholar] [CrossRef]
  14. Teixeira Alves, M.; Hilker, F.M. Hunting cooperation and Allee effects in predators. J. Theor. Biol. 2017, 419, 13–22. [Google Scholar] [CrossRef]
  15. Guan, X.; Liu, Y.U.; Xie, X. Stability analysis of a Lotka-Volterra type predator-prey system with Allee effect on the predator species. Commun. Math. Biol. Neurosci. 2018, 2018, 9. [Google Scholar] [CrossRef]
  16. Bodine, E.N.; Yust, A.E. Predator–prey dynamics with intraspecific competition and an Allee effect in the predator population. Lett. Biomath. 2017, 4, 23–38. [Google Scholar] [CrossRef]
  17. Pal, P.J.; Saha, T. Qualitative analysis of a predator-prey system with double Allee effect in prey. Chaos Solitons Fractals 2015, 73, 36–63. [Google Scholar] [CrossRef]
  18. Singh, M.K.; Bhadauria, B.S.; Singh, B.K. Bifurcation analysis of modified Leslie-Gower predator-prey model with double Allee effect. Ain Shams Eng. J. 2018, 9, 1263–1277. [Google Scholar] [CrossRef] [Green Version]
  19. Courchamp, F.; Clutton-Brock, T.; Grenfell, B. Multipack dynamics and the Allee effect in the African wild dog, Lycaon pictus. Anim. Conserv. 2000, 3, 277–285. [Google Scholar] [CrossRef]
  20. El-Shahed, M.; Ahmed, A.M.; Abdelstar, I.M.E. Dynamics of a plant-herbivore model with fractional order. Progr. Fract. Differ. Appl. 2017, 3, 59–67. [Google Scholar] [CrossRef]
  21. Suryanto, A.; Darti, I.; Panigoro, H.S.; Kilicman, A. A fractional-order predator–prey model with ratio-dependent functional response and linear harvesting. Mathematics 2019, 7, 1100. [Google Scholar] [CrossRef] [Green Version]
  22. Khoshsiar Ghaziani, R.; Alidousti, J.; Eshkaftaki, A.B. Stability and dynamics of a fractional order Leslie–Gower prey–predator model. Appl. Math. Model. 2016, 40, 2075–2086. [Google Scholar] [CrossRef]
  23. Panigoro, H.S.; Suryanto, A.; Kusumawinahyu, W.M.; Darti, I. Continuous threshold harvesting in a Gause-type predator-prey model with fractional-order. AIP Conf. Proc. 2020, 2264, 040001. [Google Scholar] [CrossRef]
  24. Panigoro, H.S.; Suryanto, A.; Kusumawinahyu, W.M.; Darti, I. Dynamics of an eco-epidemic predator–prey model involving fractional derivatives with power-law and Mittag–Leffler kernel. Symmetry 2021, 13, 785. [Google Scholar] [CrossRef]
  25. Petras, I. Fractional-Order Nonlinear Systems; Springer: London, UK, 2011. [Google Scholar]
  26. Kumar, S.; Kumar, R.; Cattani, C.; Samet, B. Chaotic behaviour of fractional predator-prey dynamical system. Chaos Solitons Fractals 2020, 135, 109811. [Google Scholar] [CrossRef]
  27. Singh, J.; Ganbari, B.; Kumar, D.; Baleanu, D. Analysis of fractional model of guava for biological pest control with memory effect. J. Adv. Res. 2021, in press. [Google Scholar] [CrossRef]
  28. El-Ajou, A.; Al-Smadi, M.; Oqielat, M.N.; Momani, S.; Hadid, S. Smooth expansion to solve high-order linear conformable fractional PDEs via residual power series method: Applications to physical and engineering equations. Ain Shams Eng. J. 2020, 11, 1243–1254. [Google Scholar] [CrossRef]
  29. Ali, H.M.; Habib, M.A.; Miah, M.M.; Akbar, M.A. Solitary wave solutions to some nonlinear fractional evolution equations in mathematical physics. Heliyon 2020, 6, e03727. [Google Scholar] [CrossRef] [PubMed]
  30. Shi, G.; Gao, J. European option pricing problems with fractional uncertain processes. Chaos Solitons Fractals 2021, 143, 110606. [Google Scholar] [CrossRef]
  31. Suryanto, A.; Darti, I.; Anam, S. Stability analysis of a fractional order modified Leslie-Gower model with additive Allee effect. Int. J. Math. Math. Sci. 2017, 2017, 8273430. [Google Scholar] [CrossRef] [Green Version]
  32. Baisad, K.; Moonchai, S. Analysis of stability and Hopf bifurcation in a fractional Gauss-type predator–prey model with Allee effect and Holling type-III functional response. Adv. Differ. Equ. 2018, 2018, 82. [Google Scholar] [CrossRef] [Green Version]
  33. Indrajaya, D.; Suryanto, A.; Alghofari, A.R. Dynamics of modified Leslie-Gower predator-prey model with Beddington-DeAngelis functional response and additive Allee effect. Int. J. Ecol. Dev. 2016, 31, 60–71. [Google Scholar]
  34. Ginzburg, L.R. Assuming reproduction to be a function of consumption raises doubts about some popular predator-prey models. J. Anim. Ecol. 1998, 67, 325–327. [Google Scholar] [CrossRef]
  35. Leslie, P.H.; Gower, J.C. The properties of a stochastic model for the predator-prey type of interaction between two species. Biometrika 1960, 47, 219–234. [Google Scholar] [CrossRef]
  36. Mallet, J. The struggle for existence: How the notion of carrying capacity, K, obscures the links between demography, Darwinian evolution, and speciation. Evol. Ecol. Res 2012, 14, 627–665. [Google Scholar]
  37. Gabriel, J.P.; Saucy, F.; Bersier, L.-F. Paradoxes in the logistic equation? Ecol. Model. 2005, 185, 147–151. [Google Scholar] [CrossRef] [Green Version]
  38. Arditi, R.; Bersier, L.-F.; Rohr, R.P. The perfect mixing paradox and the logistic equation: Verhulst vs. Lotka. Ecosphere 2016, 7, e01599. [Google Scholar] [CrossRef] [Green Version]
  39. Aziz-Alaoui, M.A.; Okiye, M.D. Boundedness and Global Stability for a Predator-Prey Model with Modified and Holling-Type II Schemes. Appl. Math. Lett. 2002, 16, 1069–1075. [Google Scholar] [CrossRef] [Green Version]
  40. Yu, S. Global stability of a modified Leslie-Gower model with Beddington-DeAngelis functional response. Adv. Differ. Equ. 2014, 2014, 84. [Google Scholar] [CrossRef] [Green Version]
  41. Vera-Damián, Y.; Vidal, C.; González-Olivares, E. Dynamics and bifurcations of a modified Leslie-Gower–type model considering a Beddington-DeAngelis functional response. Math. Methods Appl. Sci. 2019, 42, 3179–3210. [Google Scholar] [CrossRef]
  42. Melese, D.; Feyissa, S. Stability and bifurcation analysis of a diffusive modified Leslie-Gower prey-predator model with prey infection and Beddington DeAngelis functional response. Heliyon 2021, 7, e06193. [Google Scholar] [CrossRef]
  43. Arancibia-Ibarra, C.; Flores, J. Dynamics of a Leslie–Gower predator–prey model with Holling type II functional response, Allee effect and a generalist predator. Math. Comput. Simul. 2021, 188, 1–22. [Google Scholar] [CrossRef]
  44. Tiwari, B.; Raw, S.N. Dynamics of Leslie–Gower model with double Allee effect on prey and mutual interference among predators. Nonlinear Dyn. 2021, 1229–1257, 1229–1257. [Google Scholar] [CrossRef]
  45. Zhou, S.R.; Liu, Y.F.; Wang, G. The stability of predator–prey systems subject to the Allee effects. Theor. Popul. Biol. 2005, 67, 23–31. [Google Scholar] [CrossRef]
  46. Diethelm, K. A fractional calculus based model for the simulation of an outbreak of dengue fever. Nonlinear Dyn. 2013, 71, 613–619. [Google Scholar] [CrossRef]
  47. Moustafa, M.; Mohd, M.H.; Ismail, A.I.; Abdullah, F.A. Dynamical analysis of a fractional-order eco-epidemiological model with disease in prey population. Adv. Differ. Equ. 2020, 2020, 48. [Google Scholar] [CrossRef]
  48. Panigoro, H.S.; Suryanto, A.; Kusumawinahyu, W.M.; Darti, I. A Rosenzweig–MacArthur model with continuous threshold harvesting in predator involving fractional derivatives with power law and mittag–leffler kernel. Axioms 2020, 9, 122. [Google Scholar] [CrossRef]
  49. Rahmi, E.; Darti, I.; Suryanto, A.; Trisilowati; Panigoro, H.S. Stability analysis of a fractional-order Leslie-Gower model with Allee effect in predator. J. Phys. Conf. Ser. 2021, 1821, 012051. [Google Scholar] [CrossRef]
  50. Li, Y.; Chen, Y.Q.; Podlubny, I. Stability of fractional-order nonlinear dynamic systems: Lyapunov direct method and generalized Mittag-Leffler stability. Comput. Math. Appl. 2010, 59, 1810–1821. [Google Scholar] [CrossRef] [Green Version]
  51. Li, H.L.; Zhang, L.; Hu, C.; Jiang, Y.L.; Teng, Z. Dynamical analysis of a fractional-order predator-prey model incorporating a prey refuge. J. Appl. Math. Comput. 2017, 54, 435–449. [Google Scholar] [CrossRef]
  52. Ahmed, E.; El-Sayed, A.; El-Saka, H.A. On some Routh–Hurwitz conditions for fractional order differential equations and their applications in Lorenz, Rössler, Chua and Chen systems. Phys. Lett. A 2006, 358, 1–4. [Google Scholar] [CrossRef]
  53. Li, X.; Wu, R. Hopf bifurcation analysis of a new commensurate fractional-order hyperchaotic system. Nonlinear Dyn. 2014, 78, 279–288. [Google Scholar] [CrossRef]
  54. Vargas-De-León, C. Volterra-type Lyapunov functions for fractional-order epidemic systems. Commun. Nonlinear. Sci. Numer. Simul. 2015, 24, 75–85. [Google Scholar] [CrossRef]
  55. Huo, J.; Zhao, H.; Zhu, L. The effect of vaccines on backward bifurcation in a fractional order HIV model. Nonlinear Anal. Real World Appl. 2015, 26, 289–305. [Google Scholar] [CrossRef]
  56. Arenas, A.J.; González-Parra, G.; Chen-Charpentier, B.M. Construction of nonstandard finite difference schemes for the SI and SIR epidemic models of fractional order. Math. Comput. Simul 2016, 121, 48–63. [Google Scholar] [CrossRef]
  57. Barman, D.; Roy, J.; Alrabaiah, H.; Panja, P.; Mondal, S.P.; Alam, S. Impact of predator incited fear and prey refuge in a fractional order prey predator model. Chaos Solitons Fractals 2020, 142, 110420. [Google Scholar] [CrossRef]
  58. Das, M.; Samanta, G.P. A prey-predator fractional order model with fear effect and group defense. Int. J. Dyn. Control 2020, 9, 334–349. [Google Scholar] [CrossRef]
  59. Alidousti, J.; Ghafari, E. Dynamic behavior of a fractional order prey-predator model with group defense. Chaos Solitons Fractals 2020, 134, 109688. [Google Scholar] [CrossRef]
  60. Jena, R.M.; Chakraverty, S.; Baleanu, D.; Jena, S.K. Analysis of time-fractional dynamical model of romantic and interpersonal relationships with non-singular kernels: A comparative study. Math. Meth. Appl. Sci. 2020, 44, 2183–2199. [Google Scholar] [CrossRef]
  61. Ullah, I.; Ahmad, S.; Rahman, M.; Arfan, M. Investigation of fractional order tuberculosis (TB) model via Caputo derivative. Chaos Solitons Fractals 2020, 142, 110479. [Google Scholar] [CrossRef]
  62. Mahmood, S.; Shah, R.; Khan, H.; Arif, M. Laplace adomian decomposition method for multi dimensional time fractional model of Navier-Stokes equation. Symmetry 2019, 11, 149. [Google Scholar] [CrossRef] [Green Version]
  63. Diethelm, K.; Ford, N.J.; Freed, A.D. A predictor-corrector approach for the numerical solution of fractional differential equations. Nonlinear Dyn. 2002, 29, 3–22. [Google Scholar] [CrossRef]
Figure 1. Dynamics of system (6) with parameter values (26) and α = 0.98 . (a) Bifurcation diagram of system (6) driven by b for the case of the weak Allee effect ( m = 1 < 0 ) . (be) Phase-portrait of system (6) with (b) b = 0.20 , (c) b = 0.27 , (d) b = 0.54 , and (e) b = 0.57 . The predator extinction point is always unstable.
Figure 1. Dynamics of system (6) with parameter values (26) and α = 0.98 . (a) Bifurcation diagram of system (6) driven by b for the case of the weak Allee effect ( m = 1 < 0 ) . (be) Phase-portrait of system (6) with (b) b = 0.20 , (c) b = 0.27 , (d) b = 0.54 , and (e) b = 0.57 . The predator extinction point is always unstable.
Fractalfract 05 00084 g001
Figure 2. Dynamics of system (6) with parameter values (27) and α = 0.98 . (a) Bifurcation diagram of system (6) driven by b for the case of strong Allee effect ( m = 0.4 > 0 ) . (be) Phase-portrait of system (6) with (b) b = 0.20 , (c) b = 0.40 , (d) b = 0.50 , and (e) b = 0.60 . The predator extinction point is always locally stable.
Figure 2. Dynamics of system (6) with parameter values (27) and α = 0.98 . (a) Bifurcation diagram of system (6) driven by b for the case of strong Allee effect ( m = 0.4 > 0 ) . (be) Phase-portrait of system (6) with (b) b = 0.20 , (c) b = 0.40 , (d) b = 0.50 , and (e) b = 0.60 . The predator extinction point is always locally stable.
Fractalfract 05 00084 g002
Figure 3. Dynamics of system (6) with parameter values (28) and α = 0.98 . (a) Bifurcation diagram of system (6) driven by m. (be) Phase-portrait of system (6) with (b) m = 1.79 , (c) m = 0.46 , (d) m = 0.40 , and (e) m = 0.90 . The predator extinction point is unstable for the weak Allee effect case ( m < 0 ) and is locally asymptotically stable for the strong Allee effect case ( m < 0 ) .
Figure 3. Dynamics of system (6) with parameter values (28) and α = 0.98 . (a) Bifurcation diagram of system (6) driven by m. (be) Phase-portrait of system (6) with (b) m = 1.79 , (c) m = 0.46 , (d) m = 0.40 , and (e) m = 0.90 . The predator extinction point is unstable for the weak Allee effect case ( m < 0 ) and is locally asymptotically stable for the strong Allee effect case ( m < 0 ) .
Fractalfract 05 00084 g003
Figure 4. Dynamics of system (6) with parameter values (29). (a) Bifurcation diagram of system (6) driven by α . (b,c) Phase-portrait of system (6) with (b) α = 0.73 and (c) α = 0.89 .
Figure 4. Dynamics of system (6) with parameter values (29). (a) Bifurcation diagram of system (6) driven by α . (b,c) Phase-portrait of system (6) with (b) α = 0.73 and (c) α = 0.89 .
Fractalfract 05 00084 g004
Figure 5. Time series of system (6) with parameter values (29) and the initial conditions taken around the interior point S ^ = ( 0.38825 , 1.80915 ) .
Figure 5. Time series of system (6) with parameter values (29) and the initial conditions taken around the interior point S ^ = ( 0.38825 , 1.80915 ) .
Fractalfract 05 00084 g005
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Rahmi, E.; Darti, I.; Suryanto, A.; Trisilowati. A Modified Leslie–Gower Model Incorporating Beddington–DeAngelis Functional Response, Double Allee Effect and Memory Effect. Fractal Fract. 2021, 5, 84. https://doi.org/10.3390/fractalfract5030084

AMA Style

Rahmi E, Darti I, Suryanto A, Trisilowati. A Modified Leslie–Gower Model Incorporating Beddington–DeAngelis Functional Response, Double Allee Effect and Memory Effect. Fractal and Fractional. 2021; 5(3):84. https://doi.org/10.3390/fractalfract5030084

Chicago/Turabian Style

Rahmi, Emli, Isnani Darti, Agus Suryanto, and Trisilowati. 2021. "A Modified Leslie–Gower Model Incorporating Beddington–DeAngelis Functional Response, Double Allee Effect and Memory Effect" Fractal and Fractional 5, no. 3: 84. https://doi.org/10.3390/fractalfract5030084

APA Style

Rahmi, E., Darti, I., Suryanto, A., & Trisilowati. (2021). A Modified Leslie–Gower Model Incorporating Beddington–DeAngelis Functional Response, Double Allee Effect and Memory Effect. Fractal and Fractional, 5(3), 84. https://doi.org/10.3390/fractalfract5030084

Article Metrics

Back to TopTop