1. Introduction
The ecosystem encompasses numerous interactions between species that are dynamic and complex in nature. One such interaction is the ever-changing relationship between predators and their prey. The study of dynamical behavior of predator–prey interactions is essential for understanding various aspects like ecosystem stability, ecological health, resource management, and the balance of natural systems. Lotka and Volterra [
1,
2] were the first to develop a model that represented the relationship shared by prey with its predator, which was later developed by various researchers to make the model more realistic. One such proposed model is the Leslie–Gower predator–prey (LGPP) model [
3], which incorporates the concept that a predator’s carrying capacity is a function of the prey population, restricting the increase in the predator’s growth rate. The coupled differential equations that govern the LGPP model are as follows:
with the initial conditions
and
. Here,
and
are the densities of the prey and predator at time
T, respectively, and
represents the functional response, depicting the rate at which prey are consumed by the predator per unit time for a given number of prey and predators. The parameters
r,
s,
K, and
b are positive and represent the intrinsic growth rate of the prey, intrinsic growth rate of the predator, environmental carrying capacity, and maximal per capita reduction rate, respectively.
The prey population
N is not abundant in the environment, so, in times of scarcity, the predator can switch to alternate food sources [
4]. So, the model itself is modified in various ways to account for this. This phenomenon can be represented in the model by adding a constant (positive) to the denominator of the Leslie–Gower term
[
5,
6]. So, the model (
1) is transformed into the following system:
with the initial conditions
and
. Here,
represents the extent of protection provided by the environment to the predators.
The functional response
can take various forms, among which Holling type-II
has been studied widely, as it best describes the feeding pattern of many predators by including the handling time they require to consume their prey [
6,
7,
8,
9,
10]. Introducing the Holling type-II functional response in the above model gives
with the initial conditions
and
. Here,
and
represent the maximal predator per capita consumption rate and the extent of protection provided by the environment to prey. This model has also been studied under the condition that environmental protections given to the prey and predators are the same, that is,
[
6,
11]. Ji et al. [
6] studied how this model will behave in the long term with stochastic perturbation. Gupta and Chandra [
11] undertook a bifurcation analysis of a modified LGPP model with non-linear harvesting.
The Allee effect [
12], coined by ecologist W.C. Allee, is associated with the occurrence of a decrease in the fitness of the population at low densities. It is often the result of a reduction in cooperative behavior, like mating, herding, and parental care, among individuals of a species. In the literature, this phenomenon is divided into two types: strong and weak. The key that differentiates between these two is the existence of a threshold or critical value. In the weak Allee effect, the absence of a threshold value leads to the growth rate staying positive all the time, while in the strong Allee effect, the growth rate becomes negative as the population goes below the critical value [
13,
14]. This phenomenon has been analyzed by many researchers, as it plays a crucial role in determining the survival or extinction of a species. The possibility of the existence of more than one Allee effect has also been explored by ecologists, and the joint influence that a few such phenomena have on the species is named the multiple Allee effect [
15,
16,
17]. González-Olivares et al. [
18] studied the LGPP model when the Allee effect is introduced in prey. They emphasized how the Allee effect changed the behavior of the original model in terms of the existence of equilibrium points, along with their stability. Cai et al. [
19] studied the LGPP model with the additive Allee effect and discovered how a new steady state comes into actuality when this phenomenon is introduced in the model. Feng and Kang [
20] focused their interest on studying stability and bifurcation when the Allee effect is introduced both in predator and prey species in the modified LGPP model. Pal and Saha [
21] studied a ratio-dependent predator–prey model with the double Allee effect. Singh et al. [
22] analyzed a modified LGPP model under the influence of strong and weak Allee effects. They established conditions for all possible equilibrium points and local bifurcations. Rahmi et al. [
23] analyzed the Beddington–DeAngelis functional response, memory effect, and double Allee effect in a modified LGPP model. Wang et al. [
24] explored a predator–prey model with the double Allee effect and self-diffusion terms. Various researchers studied the Allee effect phenomenon in predator–prey systems, along with other phenomena, like prey refuge [
25], delay [
26], non-local competition [
27], and hunting cooperation [
28].
For a very long time, the sole focus of researchers studying the impact that predators have on prey was direct killing. But over time, it was established that the prey population can be affected by just the mere presence of predators [
29,
30,
31]. Predators induce fear in their prey, which leads the prey to display anti-predator behavior in various forms, like foraging, herding, and abandoning high-risk areas [
30,
31]. So, along with direct predation, indirect interaction with the predator also alters the behavior of prey [
32], which makes the study of this phenomenon, the fear effect, crucial.
Wang et al. [
33] modeled the fear effect in the predator–prey model and established how the system becomes stable barring the existence of limit cycles with a high level of fear, while a low level of fear leads to multiple limit cycles, leading to bistability. Sarkar and Khajanchi [
34] studied the impact that the fear effect has on prey because of the predator in a predator–prey interaction. A ratio-dependent LGPP model with fear and the Allee effect was examined by Li et al. [
35]. Chen et al. [
36] studied a modified LGPP model with the fear effect. Wu et al. [
37] studied the fear effect, along with non-linear harvesting, in the LGPP model. They established the occurrence of two limit cycles, along with saddle-node and Bogdanov–Takens bifurcations of codimension three. A study of the LGPP model incorporating the fear effect was also undertaken, along with various phenomena like harvesting [
38], prey refuge [
39,
40,
41], disease transmission [
42], and hunting cooperation [
43,
44].
Singh et al. [
22] concluded in their study that the double Allee effect can accelerate the extinction of vulnerable species. Also, Preisser and Bolnick [
32] commented on how indirect interaction with the predator, that is, the presence of fear, can alter the behavior of prey. A study where both of these phenomena are explored simultaneously has not been conducted thus far. Therefore, exploring both effects together can uncover complex dynamics of predator and prey interactions. In this article, our aim is to unveil the dynamical complexity of a modified LGPP model when both of these phenomena occur together.
This entire work contains seven sections in total.
Section 2 includes the formulation of the proposed model by the incorporation of the double Allee effect and fear effect. It also presents the positivity and boundedness of the model. Following this,
Section 3 comprises the computation of the equilibrium points that emerge in the proposed system, along with their stability analysis.
Section 4 includes the emergence of various bifurcations, namely, saddle-node, Hopf, and Bogadanov–Takens bifurcations.
Section 5 comprises the numerical simulation and the phase portraits for the proposed system. Following this,
Section 6 delves into the impact of the Allee effect and fear effect parameters. This work concludes with
Section 7, where we discuss the dynamic behavior of the proposed model in detail.
3. Existence of Equilibrium Points and Local Stability
This section emphasizes establishing the presence of the equilibria of system (
5) and analyzes their local stability. Studying the system in the case of the strong Allee effect
, we find that system (
5) exhibits axial and interior equilibria.
All possible equilibria that the system can possess are as follows:
Axial: along with the origin, that is,
, system (
5) exhibits three other axial equilibrium points:
,
, and
.
Interior: The abscissa and ordinate of the interior equilibria can be acquired by solving the subsequent equations:
The ordinate of the interior equilibrium points is given by (
9). The discriminant of Equation (
8) is
. For the ease of calculation, consider
,
, and
, where
. For the feasibility of positive equilibrium points, we must have
,
, and
. Under these conditions, system (
5) exhibits the following:
Two equilibrium points and if , where , , , and .
A unique equilibrium point if , where and .
No equilibrium point if .
Now, we shift our focus toward studying the local stability of each equilibrium point. For system (
5), the Jacobian matrix at any given point
E comes out to be
where
,
,
, and
.
Theorem 1. The equilibrium point is consistently a saddle, is consistently an asymptotically stable point, is consistently a saddle, and is consistently unstable.
Proof. From the Jacobian matrix given in (
10), we have the following:
The eigenvalues at point are and , which clearly indicate that is a saddle.
The eigenvalues at point are − and , which clearly indicate that is an asymptotically stable point.
The eigenvalues at point are and , which clearly indicate that is a saddle, as .
The eigenvalues at point are and , which clearly indicate that is an unstable point, as .
□
Theorem 2. The equilibrium points , , and , if they exist, are stable if , consistently a saddle point, and a degenerate singularity, respectively.
Proof. The equilibrium points
,
, and
are obtained from the subsequent equations:
Using the above two equations and the Jacobian matrix in (
10), we have the determinant and trace as follows:
and
It is clear that at , . Now, the condition leads to . So, is asymptotically stable. At , . Therefore, constitutes a saddle point. And at , . Therefore, it is evident that constitutes a degenerate singularity. □
From the aforementioned theorem, we established that the point
is identified as a degenerate singularity. Hence, we proceed to explore system (
5)’s dynamics in the vicinity of this point.
Theorem 3. If it exists, the positive equilibrium point is characterized as follows:
- a.
A saddle node when holds.
- b.
A cusp of codimension two when , , and holds.
Proof. Consider the transformation
and
to shift
to the origin; then, using a Taylor series, system (
5) reduces to
where
,
,
, and
.
Now, the following two cases arise:
- a.
If , then and . This implies that one eigenvalue of the Jacobian matrix is zero, while the other is non-zero, which indicates that is a saddle node.
- b.
If
, then system (
11) reduces to
Introducing variables
, system (
12) becomes
where
,
, and
.
Taking
and
, system (
13) converts to
where
and
.
Consider the transformation
and
; system (
14) then reduces to
Consider the final transformation , and
; then, system (
15) reduces to
Now, if and , then in the plane, the origin presents as a cusp of codimension two, which corresponds to being a codimension-two cusp in the -plane. □
4. Bifurcation Analysis
This section emphasizes analyzing various bifurcations, such as saddle-node, Hopf, and Bogdanov–Takens bifurcations, that system (
5) may experience.
4.1. Saddle-Node Bifurcation
In
Section 3, we show that if
, system (
5) possesses two positive equilibria:
and
. And when
, the two equilibrium points coincide and give rise to the unique equilibrium point
. Also, system (
5) exhibits no equilibrium point when
. This kind of change, that is, the change in the number of equilibrium points in a system, occurs due to an inherent property called saddle-node bifurcation. Consider
. Treating
as a bifurcation parameter, to guarantee the emergence of a saddle-node bifurcation in system (
5), we use Sotomayor’s theorem.
Theorem 4. A saddle-node bifurcation emerges in system (5) around the unique interior equilibrium point with respect to parameter θ. Proof. Evidently . Thus, the Jacobian matrix has one eigenvalue that is zero. We suppose two eigenvectors, say, V and W, corresponding to matrix’s and , zero eigenvalue, respectively, then,
Now, we have , where , and , where , as
Thus, the transversality condition required for a saddle-node bifurcation is satisfied. □
4.2. Hopf Bifurcation
In
Section 3, we show that
, if it exists, is consistently a saddle point and
, if it exists, is stable whenever
. Also,
and
when
. Then, a purely imaginary nature is exhibited by the eigenvalues of
. Therefore,
will manifest as either a center or a weak focus.
Theorem 5. A Hopf bifurcation emerges at equilibrium point in system (
5)
if , and its direction is supercritical (or subcritical) if (or ). Proof. To establish the emergence of a Hopf bifurcation in system (
5), we show that the transversality conditions are fulfilled. Consider
as a bifurcation parameter; then, a critical value
exists for which
,
, and
.
Due to the occurrence of Hopf bifurcation, a limit cycle will arise in the system. We now explore the stability of the emergent limit cycles. We start by shifting
to the origin by considering
and
. Thus, system (
5) is reduced to the following (retaining
x and
y in place of
and
, respectively):
System (
17) can be represented as
, where
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
, and
.
The characteristic roots of the Jacobian matrix
J are
if
, where
and
.
and
if
, which implies that the eigenvalues
and
are purely imaginary. For further computation, due to symmetry, it is sufficient to take only one root of the characteristic equation of
J:
. We define a matrix
such that
is the eigenvector of
J corresponding to
, where
and
. The inverse of matrix
P is
. Consider the transformation
= T
, where we obtain
and
. Thus, system (
17) reduces to
where
,
, and
.
Rewriting system (
18) in polar coordinates, we obtain
Using Taylor’s expansion, at
, system (
19) reduces to
The stability of the Hopf bifurcation relies on the sign of
where
,
,
,
, , ,
, , , ,
, and
.
After simplification, we have
□
4.3. Bogdanov–Takens Bifurcation
From Theorem (3) in
Section 3, it is evident that point
is a cusp of codimension two whenever
,
, and
. There might be a chance of occurrence of a Bogdanov–Takens bifurcation of codimension two. The intersection of saddle-node and Hopf bifurcation curves is called the BT-point. Now, to explore the emergence of the Bogdanov–Takens bifurcation,
and
are designated as the bifurcation parameters.
Theorem 6. A Bogdanov–Takens bifurcation of codimension two emerges in system (
5)
around the equilibrium point in regard to parameters θ and ρ whenever , , and . Moreover, three bifurcation curves in the -plane exist through the BT-point and are given by the following:
- 1.
The saddle-node bifurcation curve:
;
- 2.
The Hopf bifurcation curve:
;
- 3.
The homoclinic bifurcation curve:
.
Proof. With the assumption that the bifurcation parameters
and
undergo small variations within the vicinity of the BT-point
, we move ahead with the analysis. Consider a point
in the vicinity of the BT-point
), where
and
are small. Then, system (
5) becomes
With respect to the variables
x and
y, system (
22) exhibits
smoothness in the vicinity of
.
Taking
and
, system (
22) turns out to be
where
,
,
,
,
, ,
, , , , and ; furthermore, and are power series in with powers and satisfying , and the coefficients smoothly depend upon and .
Now, we use the affine transformations
and
such that system (
23) becomes
where
,
,
,
,
, , , , , and ; furthermore, and are power series in with powers and satisfying , and the coefficients smoothly depend upon and .
Next, we perform the
transformations
and
such that system (
24) becomes
where
,
,
,
,
, and
; furthermore,
and
are power series in
with powers
and
satisfying
, and the coefficients smoothly depend upon
and
.
Now, consider a new time variable
where
. Rewriting
as
t, system (
25) can be rewritten as
Next, using the transformations
and
, system (
26) transforms to
where
,
,
,
, and
; furthermore,
and
are power series in
with powers
and
satisfying
, and the coefficients have a smooth dependence upon
and
.
When and are small, determination of the sign of is not possible. Hence, it is crucial to examine the following two cases:
Case 1:
Consider
,
, and
; then, system (
27) becomes
where
,
,
, and
; furthermore,
is a power series in
with powers
and
satisfying
, and the coefficients smoothly depend upon
and
.
Next, we use the affine transformations
and
; then, system (
28) becomes
where
,
, and
; furthermore,
is a power series in
with powers
and
satisfying
, and the coefficients smoothly depend upon
and
.
Making the change of variables
,
, and
, system (
29) becomes
where
furthermore,
is a power series in
with powers
and
satisfying
, and the coefficients smoothly depend upon
and
.
Case 2:
Consider
,
, and
; then, system (
27) becomes
where
,
,
, and
; furthermore,
is a power series in
with powers
and
satisfying
, and the coefficients smoothly depend upon
and
.
Next, we use the affine transformations
and
; then, system (
32) becomes
where
,
, and
; furthermore,
is a power series in
with powers
and
satisfying
, and the coefficients smoothly depend upon
and
.
Consider new variables
,
, and
; then, system (
33) becomes
where
furthermore,
is a power series in
with powers
and
satisfying
, and the coefficients smoothly depend upon
and
.
The number of cases can be minimized by using
and
to denote
and
in (
35). If
, that is, it is non-singular, then the parameter transformations in (
31) and (
35) are topologically equivalent in the small neighborhood of origin. So, it can be deduced that system (
22) experiences a Bogdanov–Takens bifurcation when
is in the vicinity of the origin by the result in Perko [
45]. The depictions of the bifurcation curves are as follows:
The saddle-node bifurcation curve ;
The Hopf bifurcation curve ;
The homoclinic bifurcation curve .
□
5. Numerical Simulation
To support and visualize all the above analysis, numerical simulation was undertaken using MATLAB R2013a and Mathematica version 12.0, and the phase portraits were drawn using MATLAB R2013a. The non-linear equations were solved using the “ode45” function present in MATLAB R2013a and the “NSolve” function present in Mathematica version 12.
We fixed the parameters as
,
,
,
, and
in system (
5). The acquired results are presented in
Table 1.
Table 1.
Number and characteristics of equilibria that resulted from the emergence of a saddle-node bifurcation.
Table 1.
Number and characteristics of equilibria that resulted from the emergence of a saddle-node bifurcation.
| Equilibrium Points (Stability) | Figure |
---|
Axial | Interior |
---|
| (saddle) | (stable), (saddle) | Figure 1a |
| (stable) | (degenerate singularity), emergence of saddle-node bifurcation | Figure 1b,d |
| (unstable), (saddle) | None | Figure 1c |
Figure 1.
Existence of equilibrium points. (
a) At
. (
b) At
. (
c) At
. The blue and red (separatrix) curves represent the solution trajectories of the system
5 and the arrows depict the direction of these trajectories in (
a–
c). (
d) Saddle-node bifurcation diagram (upper and lower curves represent stable equilibria and unstable equilibria, respectively).
Figure 1.
Existence of equilibrium points. (
a) At
. (
b) At
. (
c) At
. The blue and red (separatrix) curves represent the solution trajectories of the system
5 and the arrows depict the direction of these trajectories in (
a–
c). (
d) Saddle-node bifurcation diagram (upper and lower curves represent stable equilibria and unstable equilibria, respectively).
The occurrence of a limit cycle around the stable point, which led to a Hopf bifurcation, is shown in
Figure 2. We fixed the parameters as
,
,
,
, and
; the obtained results are in
Table 2.
Fixing the parameters as
,
,
, and
in (
22), we first found the threshold values, which come out to be
and
. Corresponding to these parameters, system (
5) exhibited a unique equilibrium point
. By incorporating all the specified parameters, we derived the subsequent system:
To transfer this equilibrium point to the origin, we made the transformations
and
. Then, we made the affine transformations
and
. So, the system reduced to
where
,
,
,
,
,
,
,
,
, and
; furthermore,
and
were power series in
with powers
and
satisfying
, and the coefficients had a smooth dependence upon
and
.
Next, we performed the
transformations
and
We made a change in the time variable
and rewrote
T as
t. Then, using the transformations
and
, system (
37) is reduced to
where
,
,
,
, and
; furthermore,
was a power series in
with powers
and
satisfying
, and the coefficients smoothly depended upon
and
.
When
and
were very small, then
. Now, we first considered
,
, and
; then, we used
and
. Then, making a change in the variables
,
, and
, the system reduced to
where
, , and ; furthermore, was a power series in with powers and that satisfied , and the coefficients smoothly depended upon and .
The determinant
; from this, it is evident that the rank of the matrix
E was 2. So, the parametric transformation (
31) was non-singular.
In the neighborhood of the BT-point
, the three bifurcation curves divided the
-plane into four distinct regions, as shown in
Figure 3a. The blue line represents the saddle-node bifurcation curve, the red line represents the Hopf bifurcation curve, and the green line represents the homoclinic bifurcation curve.
Behavior of the system in different regions:
When
, then system (
36) possessed a unique positive equilibrium point that was identified as a codimension-two cusp, as shown in
Figure 3b. When we varied the values of
and
in such a way that we entered region 1, system (
36) exhibited the absence of an interior equilibrium point, which implies that the system moved toward a prey-free equilibrium point, as shown in
Figure 3c. As
and
entered region 2 by crossing the saddle-node bifurcation curve, the system went from exhibiting no interior equilibrium point to exhibiting two interior equilibrium points. One of the equilibrium points was consistently a saddle, while the other was characterized as unstable, as shown in
Figure 3d. As the values of
and
entered region 3 by crossing the Hopf bifurcation curve, the positive equilibrium point was enclosed within an unstable limit cycle, as shown in
Figure 3e. As the values of
and
entered region 4 by crossing the homoclinic bifurcation curve, the system lost the limit cycle and attained an interior equilibrium point that was stable, along with a saddle interior equilibrium point (cf.
Figure 3f). The behavior of the system
5 in each region has been comprised in
Table 3.
7. Discussion
The conservation of the ecosystem relies on our capability to comprehend the delicate and dynamic intricacy of ecology. A better comprehension of the dynamics enables us to recognize potential threats to ecological stability. The ecosystem comprises various interactions between species that are complex and dynamic. In this work, we studied one such complex part of the ecosystem, which is the predator–prey interaction. For this purpose, we modified the system presented by Singh et al. [
22], which focuses on the study of the double Allee effect in the modified LGPP model. In their study, the authors considered direct killing as the only way in which predator populations alter the population of prey. In the current paper, we improved the realism of this system by incorporating the fear effect.
The relevant properties for the feasibility of the proposed model were studied. Analytically and numerically, the proposed system exhibits at most six equilibrium points, including both axial and interior equilibria. The existence of axial equilibria remains unaffected whatever the values of parameters, while the system admits zero, one, or two interior equilibria depending on the parametric conditions. After examining the stability, it is seen that the trivial equilibrium point is consistently a saddle; the prey-free equilibrium point is an asymptotically stable point; and one of the two predator-free equilibrium points is consistently a saddle, while the other is an unstable point. In terms of ecology, it can be said that the predator and prey population cannot go extinct simultaneously, but in some cases, the prey species will move toward extinction while predators will never go extinct. If there are two interior equilibrium points, one is consistently a saddle, while the other will switch between a stable point, unstable point, or surrounded by an unstable limit cycle depending on certain parametric restrictions. Ecologically speaking, this means that either both prey and predator populations will coexist, prey will cease to exist, or there will be oscillating behavior. If there is a unique interior equilibrium point, then it is either a saddle node or a cusp of codimension two, contingent on the parametric restrictions. Ecologically, we can say that either the species will coexist or the prey will go extinct, contingent upon the initial population that the system commenced with. The initial values of both species has a great influence in determining the future dynamics of the system.
The diverse bifurcations exhibited by the proposed system were also analyzed. The phenomenon under study, that is, the fear effect, plays a vital role in the occurrence of the bifurcation. The possibility of the existence of zero, one, or two positive equilibrium points depends on the critical value of the fear effect parameter in the proposed model. The emergence of saddle-node bifurcation was proved by using Sotomayer’s theorem. Ecologically, both populations will coexist if does not cross the critical value , and if the critical value is crossed, then the prey will go extinct; that is, after a certain degree of fear, the prey population will tend to extinction. The proposed model also exhibits a limit cycle because of the emergence of a Hopf bifurcation. The stability of this limit cycle and the presence of a homoclinic loop (formed by the convergence of a saddle point and limit cycle) were also studied numerically. The system also possesses a Bogdanov–Takens bifurcation of codimension two and is explored using both analytical methods and numerical simulation. Ecologically speaking, this bifurcation leads to the emergence of various regions in the vicinity of the BT-point, each possessing unique properties, like the extinction of prey species, the coexistence of both species, or the oscillating behavior of prey and predator.
The influence of fear and the double Allee effect on the species’ growth rate was also examined. The analysis of the system was undertaken using two scenarios: considering the predator as constant and then variable. Furthermore, to dig deeper into the study of the fear effect, the influence of fear with each Allee effect parameter, individually, was also explored. It was observed as the Allee effect parameters and fear effect were increased, the stable equilibrium point moved downward. And as these parameters crossed their critical value, the system moved toward the prey-free equilibrium point. From an ecological point of view, it can be said that the proposed system is highly dependent on the double Allee effect and fear effect parameters and the prey species can cease to exist if the degree of these phenomena is increased beyond certain values of these parameters.
In conclusion, the present study can help in better understanding the predator–prey system, as it incorporates the fear effect, which takes us one step closer to a more realistic ecological model.