Next Article in Journal
Generalized Galilean Rotations
Previous Article in Journal
Data-Sharing System with Attribute-Based Encryption in Blockchain and Privacy Computing
Previous Article in Special Issue
Bilateral Correlational Behavior of Pyroglutamate Aminopeptidase I Activity in Rat Photoneuroendocrine Locations During a Standard 12:12 h Light–Dark Cycle
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Dynamical Complexity of Modified Leslie–Gower Predator–Prey Model Incorporating Double Allee Effect and Fear Effect

by
Manoj Kumar Singh
1,
Arushi Sharma
1 and
Luis M. Sánchez-Ruiz
2,*
1
Department of Mathematics and Statistics, Banasthali Vidyapith, Tonk 304022, India
2
Departamento de Matemática Aplicada, Universitat Politècnica de València, 46022 Valencia, Spain
*
Author to whom correspondence should be addressed.
Symmetry 2024, 16(11), 1552; https://doi.org/10.3390/sym16111552
Submission received: 18 September 2024 / Revised: 1 November 2024 / Accepted: 12 November 2024 / Published: 19 November 2024
(This article belongs to the Special Issue Symmetry/Asymmetry in Life Sciences: Feature Papers 2024)

Abstract

:
This contribution concerns studying a realistic predator–prey interaction, which was achieved by virtue of formulating a modified Leslie–Gower predator–prey model under the influence of the double Allee effect and fear effect in the prey species. The initial theoretical work sheds light on the relevant properties of the solution, presence, and local stability of the equilibria. Both analytic and numerical approaches were used to address the emergence of diverse bifurcations, like saddle-node, Hopf, and Bogdanov–Takens bifurcations. It is noteworthy that while making the assumption that the characteristic equation of the Jacobian matrix J has a pair of imaginary roots C ( ρ ) ± ι D ( ρ ) , it is sufficient to consider only C ( ρ ) + ι D ( ρ ) due to symmetry. The impact of the fear effect on the proposed model is discussed. Numerical simulation results are provided to back up all the theoretical analysis. From the findings, it was established that the initial condition of the population, as well as the phenomena (fear effect) introduced, played a crucial role in determining the stability of the proposed model.

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:
d N d T = r N 1 N K f ( N ) P , d P d T = s P 1 b P N ,
with the initial conditions N ( 0 ) > 0 and P ( 0 ) > 0 . Here, N N ( T ) and P P ( T ) are the densities of the prey and predator at time T, respectively, and f ( N ) 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 b P N [5,6]. So, the model (1) is transformed into the following system:
d N d T = r N 1 N K f ( N ) P , d P d T = s P 1 b P a 1 + N ,
with the initial conditions N ( 0 ) > 0 and P ( 0 ) > 0 . Here, a 1 represents the extent of protection provided by the environment to the predators.
The functional response f ( N ) can take various forms, among which Holling type-II ( α N a + N ) 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
d N d T = r N 1 N K α N P a 2 + N , d P d T = s P 1 b P a 1 + N ,
with the initial conditions N ( 0 ) > 0 and P ( 0 ) > 0 . Here, α and a 2 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, a 1 = a 2 = a [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.

2. Mathematical Model

We studied model (3) with a double Allee effect and fear effect in the prey species, where we have
d N d T = r N 1 + f P 1 N K N m N + n α N P a + N , d P d T = s P 1 b P a + N ,
with the initial conditions N ( 0 ) > 0 and P ( 0 ) > 0 . Here, K < m < K and n > 0 represent the Allee threshold and auxiliary parameter with m > n . Parameter m has a significant influence in determining whether the prey population will survive or face extinction. When the population size goes below the threshold value of m, the species will move toward extinction. The parameter f indicates the degree of fear induced in the prey by the predator.
To investigate the model given in system (4) further, we introduce the non-dimensional variables x, y, and t defined as
N = K x , P = K b y , t = r T .
Using the new change of variables, system (4) becomes
d x d t = x ( 1 x ) ( x β ) ( 1 + θ y ) ( x + γ ) ξ x y x + δ , d y d t = ρ y 1 y x + δ ,
with the initial conditions x ( 0 ) > 0 and y ( 0 ) > 0 . Here, β = m K , γ = n K , θ = f K b , ξ = α b r , ρ = s r , and δ = a K . The parameter θ is directly proportional to the fear effect parameter f. The system is defined on the set Ω = { ( x , y ) R + 2 : x 0 , y 0 } .

Positivity and Boundedness

This part of the article demonstrates that the solution ( x ( t ) , y ( t ) ) of system (5) consistently maintains positivity and boundedness. Ecologically, this is significant to prove, as the population of a species can neither be negative nor infinite.
Lemma 1.
  • Under the initial conditions x ( 0 ) > 0 and y ( 0 ) > 0 , every solution ( x ( t ) , y ( t ) ) of system (5) is positive.
  • Under the initial conditions x ( 0 ) > 0 and y ( 0 ) > 0 , every solution ( x ( t ) , y ( t ) ) of system (5) is bounded.
Proof. 
  • To prove this, we start by integrating the equations in system (5) under the initial conditions x ( 0 ) > 0 and y ( 0 ) > 0 to obtain the solution ( x ( t ) , y ( t ) ) . So, system (5)’s solution can be expressed as
    x ( t ) = x ( 0 ) e x p 0 t ( 1 x ( s ) ) ( x ( s ) β ) ( 1 + θ y ( s ) ) ( x ( s ) + γ ) ξ y ( s ) x ( s ) + δ d s > 0 ,
    y ( t ) = y ( 0 ) e x p ρ 0 t y ( s ) 1 y ( s ) x ( s + δ ) d s > 0 .
    Therefore, we can conclude that every solution starting in the interior of set Ω , that is, I n t ( Ω ) , remains in I n t ( Ω ) for t 0 . Along with this, for all future time, any solution trajectory that starts from the positive x-axis (or y-axis) will remain on it. Hence, the set Ω is an invariant set.
  • Under the initial conditions x ( 0 ) > 0 and y ( 0 ) > 0 , if ( x ( t ) , y ( t ) ) represents any solution of the system, then we have two cases:
    Case 1: Let us suppose x ( 0 ) 1 and we make the claim that x ( t ) 1 t 1 . Otherwise, there will exist t 1 , t 2 R + such that t 2 > t 1 , x ( t 1 ) = 1 , and x ( t ) > 1 t ( t 1 , t 2 ) . Then, Equation (6) can be written as
    x ( t ) = x ( 0 ) e x p 0 t 1 ( 1 x ( s ) ) ( x ( s ) β ) ( 1 + θ y ( s ) ) ( x ( s ) + γ ) ξ y ( s ) x ( s ) + δ d s e x p t 1 t 2 ( 1 x ( s ) ) ( x ( s ) β ) ( 1 + θ y ( s ) ) ( x ( s ) + γ ) ξ y ( s ) x ( s ) + δ d s = x ( t 1 ) e x p t 1 t 2 ( 1 x ( s ) ) ( x ( s ) β ) ( 1 + θ y ( s ) ) ( x ( s ) + γ ) ξ y ( s ) x ( s ) + δ d s < x ( t 1 ) ,
    as ( 1 x ( s ) ) ( x ( s ) β ) ( 1 + θ y ( s ) ) ( x ( s ) + γ ) ξ y ( s ) x ( s ) + δ < 0 t ( t 1 , t 2 ) , which is a contradiction to our hypothesis. Thus, x ( t ) 1 for all t 0 .
    Case 2: Now, let us suppose x ( 0 ) > 1 ; then, as long as x ( t ) 1 ,
    x ( t ) = x ( 0 ) e x p 0 t ( 1 x ( s ) ) ( x ( s ) β ) ( 1 + θ y ( s ) ) ( x ( s ) + γ ) ξ y ( s ) x ( s ) + δ d s < x ( 0 ) ,
    as ( 1 x ( s ) ) ( x ( s ) β ) ( 1 + θ y ( s ) ) ( x ( s ) + γ ) ξ y ( s ) x ( s ) + δ < 0 for x ( t ) 1 .
    Hence, based on these two cases, we can deduce that every positive solution satisfies x ( t ) m a x { x ( 0 ) , 1 } M 1 t 0 .
    Next, we proceed to the predator equation of system (5):
    d y d t = ρ y 1 y x + δ ρ y 1 y M 1 + δ
    So, we have y ( t ) m a x { y ( 0 ) , M 1 + δ } t 0 , which completes our proof.

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 ( β > 0 ) , 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, E 0 = ( 0 , 0 ) , system (5) exhibits three other axial equilibrium points: E δ = ( 0 , δ ) , E 1 = ( 1 , 0 ) , and E β = ( β , 0 ) .
  • Interior: The abscissa and ordinate of the interior equilibria can be acquired by solving the subsequent equations:
    ( 1 + θ ξ ) x 2 ( 1 + β ξ ( 1 + θ δ + θ γ ) ) x + β + ξ γ ( 1 + θ δ ) = 0
    y = x + δ
    The ordinate of the interior equilibrium points is given by (9). The discriminant of Equation (8) is Δ ( θ ) = ξ 2 ( δ γ ) 2 θ 2 + 2 ξ ( 2 γ δ + 2 β + ( δ + γ ) ( 1 + β ) ξ ( δ γ ) ) θ + ( 1 + β ξ ) 2 4 β 4 γ ξ . For the ease of calculation, consider Δ 1 = ( 1 + β ξ ( 1 + θ δ + θ γ ) ,
    Δ 2 = 4 ( 1 + θ ξ ) ( β + ξ γ + θ ξ γ δ ) , and θ 1 = ( 2 γ δ + 2 β + ( δ + γ ) ( 1 + β ) ξ ( δ γ ) ) Δ * ξ ( δ γ ) 2 , where Δ * = ( 2 γ δ + 2 β + ( δ + γ ) ( 1 + β ) ξ ( δ γ ) ) 2 ( δ γ ) 2 ( ( 1 + β ξ ) 2 4 β 4 γ ξ ) . For the feasibility of positive equilibrium points, we must have Δ 1 > 0 , δ < γ , and ( 1 + β ξ ) 2 > 4 β + 4 γ ξ . Under these conditions, system (5) exhibits the following:
    • Two equilibrium points E 2 * = ( x 2 * , y 2 * ) and E 3 * = ( x 3 * , y 3 * ) if 0 < θ < θ 1 , where x 2 * = Δ 1 + Δ 1 2 Δ 2 2 ( 1 + θ ξ ) , x 3 * = Δ 1 Δ 1 2 Δ 2 2 ( 1 + θ ξ ) , y 2 * = x 2 * + δ , and y 3 * = x 3 * + δ .
    • A unique equilibrium point E * = ( x * , y * ) if θ = θ 1 , where x * = Δ 1 2 ( 1 + θ ξ ) and y * = x * + δ .
    • No equilibrium point if θ > θ 1 .
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
J E = a 11 a 12 a 21 a 22
where a 11 = x ( 1 + β 2 x ) ( 1 + θ y ) ( x + γ ) x ( 1 x ) ( x β ) ( 1 + θ y ) ( x + γ ) 2 + ξ x y ( x + δ ) 2 + ( 1 x ) ( x β ) ( 1 + θ y ) ( x + γ ) ξ y x + δ , a 12 = θ x ( 1 x ) ( x β ) ( x + γ ) ( 1 + θ y ) 2 ξ x x + δ , a 21 = ρ y 2 ( x + δ ) 2 , and a 22 = ρ 2 ρ y ( x + δ ) .
Theorem 1.
The equilibrium point E 0 is consistently a saddle, E δ is consistently an asymptotically stable point, E 1 is consistently a saddle, and E β is consistently unstable.
Proof. 
From the Jacobian matrix given in (10), we have the following:
  • The eigenvalues at point E 0 are β γ and ρ , which clearly indicate that E 0 is a saddle.
  • The eigenvalues at point E δ are − ( β γ ( 1 + δ θ ) + ξ ) and ρ , which clearly indicate that E δ is an asymptotically stable point.
  • The eigenvalues at point E 1 are 1 β 1 + γ and ρ , which clearly indicate that E 1 is a saddle, as 0 < β < 1 .
  • The eigenvalues at point E β are β 1 β β + γ and ρ , which clearly indicate that E β is an unstable point, as 0 < β < 1 .
Theorem 2.
The equilibrium points E 2 * , E 3 * , and E * , if they exist, are stable if x 2 * 1 + β 2 x 2 * ξ ( 1 + θ y 2 * ) ( 1 + θ y 2 * ) ( x 2 * + γ ) + ξ x 2 * x 2 * + δ ρ < 0 , consistently a saddle point, and a degenerate singularity, respectively.
Proof. 
The equilibrium points E 2 * , E 3 * , and E * are obtained from the subsequent equations:
( 1 x ) ( x β ) ( 1 + θ y ) ( x + γ ) ξ y x + δ = 0 , y = x + δ .
Using the above two equations and the Jacobian matrix in (10), we have the determinant and trace as follows:
d e t ( J E ) = ρ x ( 1 + θ y ) ( x + γ ) ( 1 + β ξ ( 1 + θ δ + θ γ ) 2 ( 1 + θ ξ ) x )
and
T r ( J E ) = x 1 + β 2 x ξ ( 1 + θ y ) ( 1 + θ y ) ( x + γ ) + ξ x x + δ ρ .
It is clear that at E 2 * , d e t ( J E 2 * ) > 0 . Now, the condition x 2 * 1 + β 2 x 2 * ξ ( 1 + θ y 2 * ) ( 1 + θ y 2 * ) ( x 2 * + γ ) + ξ x 2 * x 2 * + δ ρ < 0 leads to T r ( J E 2 * ) < 0 . So, E 2 * is asymptotically stable. At E 3 * , d e t ( J E 3 * ) < 0 . Therefore, E 3 * constitutes a saddle point. And at E * , d e t ( J E * ) = 0 . Therefore, it is evident that E * constitutes a degenerate singularity. □
From the aforementioned theorem, we established that the point E * 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 E * is characterized as follows:
a. 
A saddle node when x * ( 1 + β 2 x * ξ ( 1 + θ y * ) ) ( 1 + θ y * ) ( x * + γ ) + ξ x * x * + δ ρ holds.
b. 
A cusp of codimension two when x * ( 1 + β 2 x * ξ ( 1 + θ y * ) ) ( 1 + θ y * ) ( x * + γ ) + ξ x * x * + δ = ρ , α ¯ 20 0 , and ( 2 α 20 ¯ α ¯ 11 ) 0 holds.
Proof. 
Consider the transformation x ^ = x x * and y ^ = y y * to shift E * to the origin; then, using a Taylor series, system (5) reduces to
d x ^ d t = α 10 x ^ α 01 y ^ + α 20 x ^ 2 + α 02 y ^ 2 + α 11 x ^ y ^ + O | ( x , y ) 3 | d y ^ d t = ρ x ^ ρ y ^ ρ x * + δ x ^ 2 ρ x * + δ y ^ 2 + 2 ρ x * + δ x ^ y ^ + O | ( x , y ) 3 |
where α 10 = α 01 = ξ θ x * 1 + θ y * + ξ x * x * + δ , α 20 = x * ( 1 + θ y * ) ( x * + γ ) + γ ( 1 + β ξ ( 1 + θ y * ) 2 x * ) ( 1 + θ y * ) ( x * + γ ) 2 + 2 ξ δ ( x + δ ) 2 , α 02 = θ 2 ξ x * ( 1 + θ y * ) 2 , and α 11 = θ x * ( 1 2 x * + β ) ( 1 + θ y * ) 2 ( x * + γ ) θ γ ξ ( 1 + θ y * ) ( x * + γ ) ξ δ ( x * + δ ) 2 .
Now, the following two cases arise:
a.
If ξ θ x * 1 + θ y * + ξ x * x * + δ ρ , then T r ( J E * ) 0 and d e t ( J E * ) = 0 . This implies that one eigenvalue of the Jacobian matrix J E * is zero, while the other is non-zero, which indicates that E * is a saddle node.
b.
If ξ θ x * 1 + θ y * + ξ x * x * + δ = ρ , then system (11) reduces to
d x ^ d t = ρ x ^ ρ y ^ + α 20 x ^ 2 + α 02 y ^ 2 + α 11 x ^ y ^ + O | ( x , y ) 3 | , d y ^ d t = ρ x ^ ρ y ^ ρ x * + δ x ^ 2 ρ x * + δ y ^ 2 + 2 ρ x * + δ x ^ y ^ + O | ( x , y ) 3 | .
Introducing variables τ = ρ t , system (12) becomes
d x ^ d τ = x ^ y ^ + α 20 ^ x ^ 2 + α 02 ^ y ^ 2 + α 11 ^ x ^ y ^ + O | ( x , y ) 3 | , d y ^ d τ = x ^ y ^ 1 x * + δ x ^ 2 1 x * + δ y ^ 2 + 2 x * + δ x ^ y ^ + O | ( x , y ) 3 | ,
where α 20 ^ = 1 ρ α 20 , α 11 ^ = 1 ρ α 11 , and α 02 ^ = 1 ρ α 02 .
Taking x 1 = x ^ and x 2 = x ^ y ^ , system (13) converts to
d x 1 d τ = x 2 + α ¯ 20 x 1 2 + α 02 ^ x 2 2 α ¯ 11 x 1 x 2 + O | ( x 1 , x 2 ) 3 | , d x 2 d τ = α ¯ 20 x 1 2 + ( α 02 ^ + 1 x * + δ ) x 2 2 α ¯ 11 x 1 x 2 + O | ( x 1 , x 2 ) 3 | ,
where α ¯ 20 = α 20 ^ + α 11 ^ + α 02 ^ and α ¯ 11 = α 11 ^ + 2 α 02 ^ .
Consider the transformation y 1 = x 1 and y 2 = x 2 1 x * + δ x 1 x 2 ; system (14) then reduces to
d y 1 d τ = y 2 + α ¯ 20 y 1 2 + α 02 ^ y 2 2 + ( α ¯ 11 + 1 x * + δ ) y 1 y 2 + O | ( y 1 , y 2 ) 3 | , d y 2 d τ = α ¯ 20 y 1 2 + α 02 ^ y 2 2 α ¯ 11 y 1 y 2 + O | ( y 1 , y 2 ) 3 | .
Consider the final transformation z 1 = y 1 1 2 ( α ¯ 11 + 1 x * + δ ) y 1 2 + α 02 ^ y 2 2 , and
z 2 = y 2 + α ¯ 20 y 1 2 + O | ( z 1 , z 2 ) 3 | ; then, system (15) reduces to
d z 1 d τ = z 2 , d z 2 d τ = α ¯ 20 y 1 2 + α 02 ^ y 2 2 + ( 2 α 20 ¯ α ¯ 11 ) y 1 y 2 + O | ( z 1 , z 2 ) 3 | .
Now, if α ¯ 20 0 and ( 2 α 20 ¯ α ¯ 11 ) 0 , then in the z 1 z 2 plane, the origin presents as a cusp of codimension two, which corresponds to E * being a codimension-two cusp in the x y -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 0 < θ < θ 1 , system (5) possesses two positive equilibria: E 1 * and E 2 * . And when θ = θ 1 , the two equilibrium points coincide and give rise to the unique equilibrium point E * . Also, system (5) exhibits no equilibrium point when θ > θ 1 . 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 θ = θ S N = θ 1 . 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 E * = ( x * , y * ) with respect to parameter θ.
Proof. 
Evidently d e t ( J E * ) = 0 . Thus, the Jacobian matrix J E * has one eigenvalue that is zero. We suppose two eigenvectors, say, V and W, corresponding to matrix’s J E * and J E * T , zero eigenvalue, respectively, then,
V = 1 1 ; W = 1 1 ρ ( ξ θ x * 1 + θ y * + ξ x * x * + δ ) .
Now, we have W T F θ ( E * , θ S N ) = ξ ( x * + δ ) ( 1 + θ y * ) 0 , where F θ ( E * , θ S N ) = ξ ( x * + β ) ( 1 + θ y * ) 0 , and W T [ D 2 F ( E * , θ S N ) ( V , V ) ] = 2 + 2 θ δ + θ ( 1 + β + ξ ( 1 + θ δ θ γ ) ) ( 1 + θ y * ) 2 ( x * + γ ) 0 , where D 2 F ( E * , θ S N ) = 2 + 2 θ δ + θ ( 1 + β + ξ ( 1 + θ δ θ γ ) ) ( 1 + θ y * ) 2 ( x * + γ ) 0 0 , as 1 + β ξ ( ( 1 + θ δ + θ γ ) ) > 0
Thus, the transversality condition required for a saddle-node bifurcation is satisfied. □

4.2. Hopf Bifurcation

In Section 3, we show that E 3 * , if it exists, is consistently a saddle point and E 2 * , if it exists, is stable whenever x 2 * 1 + β 2 x 2 * ξ ( 1 + θ y 2 * ) ( 1 + θ y 2 * ) ( x 2 * + γ ) + ξ x 2 * x 2 * + δ < ρ . Also, d e t ( J E 2 * ) > 0 and T r ( J E 2 * ) = 0 when x 2 * 1 + β 2 x 2 * ξ ( 1 + θ y 2 * ) ( 1 + θ y 2 * ) ( x 2 * + γ ) + ξ x 2 * x 2 * + δ = ρ . Then, a purely imaginary nature is exhibited by the eigenvalues of J E 2 * . Therefore, E 2 * will manifest as either a center or a weak focus.
Theorem 5.
A Hopf bifurcation emerges at equilibrium point E 2 * in system (5) if ρ = ρ [ hf ] = x 2 * 1 + β 2 x 2 * ξ ( 1 + θ y 2 * ) ( 1 + θ y 2 * ) ( x 2 * + γ ) + ξ x 2 * x 2 * + δ , and its direction is supercritical (or subcritical) if Λ < 0 (or Λ > 0 ).
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 ρ = ρ [ hf ] = x 2 * 1 + β 2 x 2 * ξ ( 1 + θ y 2 * ) ( 1 + θ y 2 * ) ( x 2 * + γ ) + ξ x 2 * x 2 * + δ exists for which d e t ( J E 2 * ) > 0 , T r ( J E 2 * ) = 0 , and ρ ( T r ( J E 2 * ) ) = 1 0 .
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 E 2 * = ( x 2 * , y 2 * ) to the origin by considering x ^ = x x 2 * and y ^ = y y 2 * . Thus, system (5) is reduced to the following (retaining x and y in place of x ^ and y ^ , respectively):
d x d t = ( x + x 2 * ) ( 1 x x 2 * ) ( x + x 2 * β ) ( 1 + θ ( y + y 2 * ) ) ( x + x 2 * + γ ) ξ ( x + x 2 * ) ( y + y 2 * ) x + x 2 * + δ d y d t = ρ ( y + y 2 * ) ( 1 y + y 2 * x + x 2 * + δ )
System (17) can be represented as x t y t = J x y + f ( x , y ) g ( x , y ) , where J = J E 2 * = x 2 * 1 + β 2 x 2 * ξ ( 1 + θ y 2 * ) ( 1 + θ y 2 * ) ( x 2 * + γ ) + ξ x 2 * x 2 * + δ ξ θ x 2 * 1 + θ y 2 * ξ x 2 * x 2 * + δ ρ ρ , f ( x , y ) = α 20 x 2 + α 11 x y + α 02 y 2 + α 30 x 3 + α 21 x 2 y + α 12 x y 2 + α 03 y 3 + , g ( x , y ) = β 20 x 2 + β 11 x y + β 02 y 2 + β 30 x 3 + β 21 x 2 y + β 12 x y 2 + , α 20 = x ( 1 + θ y ) ( x + γ ) + γ ( 1 + β 2 x ξ ( 1 + θ y ) ) ( 1 + θ y ) ( x + γ ) 2 + ξ δ ( x + δ ) 2 , α 02 = θ 2 ξ x ( 1 + θ y ) 2 , α 11 = θ x ( 1 2 x + β ) ( 1 + θ y ) 2 ( x + γ ) θ γ ξ ( 1 + θ y ) ( x + γ ) ξ δ ( x + δ ) 2 , α 30 = γ ( 1 β + 2 x + ξ ( 1 + θ y ) ) ( 1 + θ y ) ( x + γ ) 3 γ ( 1 + θ y ) ( x + γ ) 2 δ ξ ( x + δ ) 3 , α 03 = θ 3 ξ x ( 1 + θ y ) 3 , α 21 = θ x ( 1 + θ y ) 2 ( x + γ ) θ γ ( 1 + β 2 x ξ ( 1 + θ y ) ) ( 1 + θ y ) 2 ( x + γ ) 2 + ξ δ ( x + δ ) 3 , α 12 = θ 2 x ( 1 2 x + β ) ( 1 + θ y ) 3 ( x + γ ) + θ 2 γ ξ ( 1 + θ y ) 2 ( x + γ ) , β 10 = ρ , β 01 = ρ , β 20 = b 02 = ρ x + δ , β 11 = 2 ρ x + δ , β 30 = ρ ( x + δ ) 2 , β 03 = 0 , β 21 = 2 ρ ( x + δ ) 2 , and β 12 = 2 ρ ( x + δ ) 2 .
The characteristic roots of the Jacobian matrix J are λ ( ρ ) = C ( ρ ) ± ι D ( ρ ) if d e t ( J ) > ( 1 2 T r ( J ) ) 2 , where C ( ρ ) = 1 2 T r ( J ) and D ( ρ ) = d e t ( J ) ( 1 2 T r ( J ) ) 2 . C ( ρ [ h f ] ) = 0 and D ( ρ [ h f ] ) = d e t ( J ) if ρ = ρ [ h f ] , which implies that the eigenvalues λ 1 and λ 2 are purely imaginary. For further computation, due to symmetry, it is sufficient to take only one root of the characteristic equation of J: λ = C ( ρ ) + ι D ( ρ ) . We define a matrix P = 1 0 M N such that 1 M ι N is the eigenvector of J corresponding to λ = A ( ρ ) + ι B ( ρ ) , where M = 1 ξ θ x 1 + θ y + ξ x x + δ ( x ( 1 + β ξ ( 1 + θ y ) 2 x ) ( 1 + θ y ) ( x + γ ) + ξ x x + δ A ( ρ ) ) and N = B ( ρ ) ξ θ x 1 + θ y + ξ x x + δ . The inverse of matrix P is P 1 = 1 0 M N 1 N . Consider the transformation x y = T u v , where we obtain x = u and y = M u + N v . Thus, system (17) reduces to
u t v t = J ( ρ ) u v + F ( u , v ) G ( u , v ) ,
where J ( ρ ) = C ( ρ ) D ( ρ ) D ( ρ ) C ( ρ ) , F ( u , v ) = f ( u , M u + N v ) = ( α 20 + α 11 M + α 02 M 2 ) u 2 + ( α 11 N + 2 α 02 M N ) u v + α 02 N v 2 + ( α 30 + α 21 M + α 12 M 2 + α 03 M 3 ) u 3 + ( α 21 N + 2 α 12 M N + 3 α 03 M 2 N ) u 2 v + ( α 12 N 2 + 3 α 03 M N 2 ) u v 2 + α 03 N 3 v 3 + , and G ( u , v ) = M N F ( u , v ) + 1 N g ( u , M u + N v ) = 1 N ( ( β 20 + β 11 M + β 02 M 2 α 20 M α 11 M 2 α 02 M 2 ) u 2 + ( β 11 N + 2 β 02 M N α 11 M N 2 α 02 M 2 N ) u v + ( β 02 N 2 α 02 M N ) v 2 + ( β 30 + β 21 M + β 12 M 2 + β 03 M 3 α 30 M α 21 M 2 α 12 M 3 α 03 M 4 ) u 3 + ( β 21 N + 2 β 12 M N + 3 β 03 M 2 N α 21 M N 2 α 12 M 2 N 3 α 03 M 3 N ) u 2 v + ( β 12 N + 3 β 03 M N 2 α 12 M N 2 3 α 03 M 2 N 2 ) u v 2 + ( β 03 N 2 α 03 M N 3 ) v 3 ) .
Rewriting system (18) in polar coordinates, we obtain
r ˙ = C ( ρ ) r + c ( ρ ) r 3 + , ϕ ˙ = D ( ρ ) + d ( ρ ) r 2 +
Using Taylor’s expansion, at ρ = ρ * , system (19) reduces to
r ˙ = C ( ρ * ) ( ρ ρ * ) r + c ( ρ ) r 3 + ϕ ˙ = D ( ρ * ) + D ( ρ * ) ( ρ ρ * ) + d ( ρ * ) r 2 +
The stability of the Hopf bifurcation relies on the sign of
Λ = c ( ρ * ) C ( ρ * ) ,
where c ( ρ * ) = 1 16 ( F x x x + F x y y + G x y y + G y y y ) ( 0 , 0 , ρ * ) + 1 16 B ( ρ * ) ( F x y ( F x x + F y y ) G x y ( G x x + G y y ) F x x G x x + F y y G y y ) ( 0 , 0 , ρ * ) , C ( ρ * ) = 1 , F x x = 2 ( α 20 + α 11 M + α 02 M 2 ) ,
F y y = 2 α 02 N , F x y = α 11 N + 2 α 02 M N , F x x x = 6 ( α 30 + α 21 M + α 12 M 2 + α 03 M 3 ) ,
F x y y = 2 ( α 12 N 2 + 3 α 03 M N 2 ) , G x x = 2 N ( β 20 + β 11 M + β 02 M 2 α 20 M α 11 M 2 α 02 M 2 ) , G y y = 2 ( β 02 N α 02 M ) , G x y = β 11 + 2 β 02 M α 11 M 2 α 02 M 2 ,
G x x x = 6 N ( β 30 + β 21 M + β 12 M 2 + β 03 M 3 α 30 M α 21 M 2 α 12 M 3 α 03 M 4 ) , and
G y y y = 6 ( β 03 N 2 α 03 M N 2 ) .
After simplification, we have
Λ = 1 8 ( α 30 + β 21 + B ( ρ * ) α 02 a 12 2 ( α 11 2 α 02 a 11 a 12 + 2 β 02 ) + B 2 ( ρ * ) a 12 2 ( α 21 + 3 β 03 ) 2 a 11 a 12 ( α 21 + β 21 ) + a 11 2 a 12 2 ( α 12 + 3 β 03 ) + β 11 β 02 a 12 + α 11 a 12 ( α 11 a 11 a 12 α 20 ) + α 02 a 11 a 12 2 ( 2 α 20 3 α 11 a 11 a 12 + 2 α 02 a 11 2 a 12 2 2 α 02 a 12 ) + β 02 a 11 a 12 2 ( α 11 2 β 02 2 α 02 a 11 a 12 ) ) + 1 8 B ( ρ * ) ( α 02 a 11 a 12 ( 2 β 02 a 11 a 12 α 11 a 11 a 12 + 2 α 02 a 11 2 a 12 2 β 11 ) ) + 1 8 B 2 ( ρ * ) ( β 11 ( β 02 a 12 β 11 a 11 ) + β 11 a 11 2 a 12 ( 3 β 02 α 02 ) + α 20 a 11 ( 2 α 20 β 11 ) + ( 2 α 20 a 11 2 β 02 a 11 2 a 12 α 11 a 11 2 a 12 ) ( β 02 a 11 a 12 + β 20 a 12 a 11 α 02 a 11 a 12 ) + α 11 a 11 2 a 12 ( α 11 a 11 a 12 3 α 20 + 2 β 02 a 11 a 12 ) α 20 a 11 ( 2 β 11 + 2 β 02 a 11 a 12 ) ) .

4.3. Bogdanov–Takens Bifurcation

From Theorem (3) in Section 3, it is evident that point E * is a cusp of codimension two whenever x * ( 1 + β 2 x * ξ ( 1 + θ y * ) ) ( 1 + θ y * ) ( x * + γ ) + ξ x * x * + δ ρ , α ¯ 20 0 , and ( 2 α 20 ¯ α ¯ 11 ) 0 . 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 E * = ( x * , y * ) in regard to parameters θ and ρ whenever x * ( 1 + β 2 x * ξ ( 1 + θ y * ) ) ( 1 + θ y * ) ( x * + γ )
+ ξ x * x * + δ ρ , α ¯ 20 0 , and ( 2 α 20 ¯ α ¯ 11 ) 0 . Moreover, three bifurcation curves in the λ 1 λ 2 -plane exist through the BT-point and are given by the following:
1. 
The saddle-node bifurcation curve:
S N = ( λ 1 , λ 2 ) : ζ 1 ( λ 1 , λ 2 ) = 0 , ζ 2 ( λ 1 , λ 2 ) 0 ;
2. 
The Hopf bifurcation curve:
H = ( λ 1 , λ 2 ) : ζ 2 ( λ 1 , λ 2 ) = ± ζ 1 ( λ 1 , λ 2 ) , ζ 1 ( λ 1 , λ 2 ) < 0 ;
3. 
The homoclinic bifurcation curve:
H L = ( λ 1 , λ 2 ) : ζ 2 ( λ 1 , λ 2 ) = ± 5 7 ζ 1 ( λ 1 , λ 2 ) , ζ 1 ( λ 1 , λ 2 ) < 0 .
Proof. 
With the assumption that the bifurcation parameters θ and ρ undergo small variations within the vicinity of the BT-point ( θ 0 , ρ 0 ) , we move ahead with the analysis. Consider a point ( θ 0 + λ 1 , ρ 0 + λ 2 ) in the vicinity of the BT-point ( θ 0 , ρ 0 ), where λ 1 and λ 2 are small. Then, system (5) becomes
d x d t = x ( 1 x ) ( x β ) ( 1 + ( θ + λ 1 ) y ) ( x + γ ) ξ x y x + δ , d y d t = ( ρ + λ 2 ) y ( 1 y x + δ ) .
With respect to the variables x and y, system (22) exhibits C smoothness in the vicinity of ( θ 0 , ρ 0 ) .
  • Taking u = x x * and v = y y * , system (22) turns out to be
    d u d t = α 00 + α 10 u + α 01 v + α 20 u 2 + α 11 u v + α 02 v 2 + R 1 ( u , v ) , d v d t = β 10 u + β 01 v + β 20 u 2 + β 11 u v + β 02 v 2 + R 2 ( u , v ) ,
    where α 00 = x * ( 1 x * ) ( x * β ) ( 1 + ( θ + λ 1 ) y * ) ( x * + γ ) ξ x * , α 01 = ( θ + λ 1 ) x * ( 1 x * ) ( x * β ) ( 1 + ( θ + λ 1 ) y * ) 2 ( x * + γ ) ξ x * x * + δ ,
α 10 = x * ( 1 2 x * + β ) ( 1 + ( θ + λ 1 ) y * ) ( x * + γ ) + γ ( 1 x * ) ( x * β ) ( 1 + ( θ + λ 1 ) y * ) ( x * + γ ) 2 ξ δ ( x * + δ ) 2 ,
α 20 = 2 x * ( 1 + ( θ + λ 1 ) y * ) ( x * + γ ) + 2 γ ( 1 2 x * + β ) ( 1 + ( θ + λ 1 ) y * ) ( x * + γ ) 2 2 γ ( 1 x * ) ( x * β ) ( 1 + ( θ + λ 1 ) y * ) ( x * + γ ) 3 + 2 ξ δ ( x * + δ ) 2 ,
α 11 = ( θ + λ 1 ) x * ( 1 2 x * + β ) ( 1 + ( θ + λ 1 ) y * ) 2 ( x * + γ ) ( θ + λ 1 ) γ ( 1 x * ) ( x * β ) ( 1 + ( θ + λ 1 ) y * ) 2 ( x * + γ ) 2 ξ δ ( x * + δ ) 2 , α 02 = 2 ( θ + λ 1 ) 2 x * ( 1 x ) ( 1 + ( θ + λ 1 ) y * ) 3 ( x * + γ ) ,
β 10 = ρ + λ 2 , β 01 = ( ρ + λ 2 ) , β 20 = ( ρ + λ 2 ) x + δ , β 02 = ( ρ + λ 2 ) x + δ , and β 11 = 2 ρ + λ 2 x + δ ; furthermore, R 1 and R 2 are power series in ( u , v ) with powers u i and v j satisfying i + j 3 , and the coefficients smoothly depend upon λ 1 and λ 2 .
Now, we use the affine transformations y 1 = u and y 2 = α 10 u + α 01 v such that system (23) becomes
d y 1 d t = y 2 + + ξ 00 + ξ 20 y 1 2 + ξ 11 y 1 y 2 + ξ 02 y 2 2 + R ˜ 1 ( y 1 , y 2 ) , d y 2 d t = η 00 + η 10 y 1 + η 01 y 2 + η 20 y 1 2 + η 11 y 1 y 2 + η 02 y 2 2 + R ˜ 2 ( y 1 , y 2 ) ,
where ξ 00 = α 00 , ξ 20 = α 20 α 11 α 10 α 01 + α 02 α 10 2 α 01 2 , ξ 02 = α 02 α 01 2 , ξ 11 = α 11 α 01 2 α 02 α 10 α 01 2 ,
η 00 = α 00 α 10 , η 10 = α 01 β 10 α 10 β 01 , η 01 = α 10 + β 01 , η 20 = α 20 α 10 + α 01 β 20 α 10 β 11 α 11 α 10 2 α 01 + α 10 2 β 02 α 01 + α 10 3 α 02 α 01 2 , η 02 = α 10 α 02 α 01 2 + β 02 α 01 , and η 11 = α 11 α 10 α 01 + β 11 2 β 02 α 10 α 01 2 α 02 α 10 2 α 01 2 ; furthermore, R ˜ 1 and R ˜ 2 are power series in ( y 1 , y 2 ) with powers y 1 i and y 2 j satisfying i + j 3 , and the coefficients smoothly depend upon λ 1 and λ 2 .
Next, we perform the C transformations x 1 = y 1 and x 2 = y 2 + + ξ 00 + ξ 20 y 1 2 + ξ 11 y 1 y 2 + ξ 02 y 2 2 such that system (24) becomes
d x 1 d t = x 2 + R ^ 1 ( x 1 , x 2 ) , d x 2 d t = γ 00 + γ 10 y 1 + γ 01 y 2 + γ 20 y 1 2 + γ 11 y 1 y 2 + γ 02 y 2 2 + R ^ 2 ( x 1 , x 2 ) ,
where γ 00 = η 00 ξ 00 η 01 + , γ 10 = η 10 ξ 00 η 11 + ξ 11 η 00 + , γ 01 = η 01 2 ξ 00 η 02 + ξ 02 η 00 ξ 00 ξ 11 + , γ 20 = η 20 + ξ 11 η 10 ξ 20 η 01 + , γ 01 = η 02 + ξ 11 ξ 02 η 01 + , and γ 11 = η 11 + 2 ξ 20 + 2 ξ 02 η 10 + ; furthermore, R ^ 1 and R ^ 2 are power series in ( x 1 , x 2 ) with powers x 1 i and x 2 j satisfying i + j 3 , and the coefficients smoothly depend upon λ 1 and λ 2 .
Now, consider a new time variable τ where d t = ( 1 γ 02 x 1 ) d τ . Rewriting τ as t, system (25) can be rewritten as
d x 1 d t = ( 1 γ 02 x 1 ) x 2 + R ^ 1 ( x 1 , x 2 ) , d x 2 d t = ( 1 γ 02 x 1 ) ( γ 00 + γ 10 y 1 + γ 01 y 2 + γ 20 y 1 2 + γ 11 y 1 y 2 + γ 02 y 2 2 + R ^ 2 ( x 1 , x 2 ) ) .
Next, using the transformations z 1 = x 1 and z 2 = ( 1 γ 02 x 1 ) x 2 + R ^ 1 ( x 1 , x 2 ) , system (26) transforms to
d z 1 d t = z 2 , d z 2 d t = δ 00 + δ 10 z 1 + δ 01 z 2 + δ 20 z 1 2 + δ 11 z 1 z 2 + R ¯ 2 ( z 1 , z 2 ) ,
where δ 00 = γ 00 , δ 10 = γ 10 2 γ 02 γ 00 , δ 01 = γ 01 , δ 20 = γ 20 + γ 00 γ 02 2 2 γ 02 γ 10 , and δ 11 = γ 11 γ 02 γ 01 ; furthermore, R ¯ 1 and R ¯ 2 are power series in ( z 1 , z 2 ) with powers z 1 i and z 2 j satisfying i + j 3 , and the coefficients have a smooth dependence upon λ 1 and λ 2 .
When λ 1 and λ 2 are small, determination of the sign of δ 20 is not possible. Hence, it is crucial to examine the following two cases:
Case 1: δ 20 > 0
  • Consider U = z 1 , V = z 2 δ 20 , and d τ = δ 20 d t ; then, system (27) becomes
    d U d τ = V , d V d τ = p 00 + p 10 U + p 01 V + U 2 + p 11 U V + R 3 ( U , V ) ,
    where p 00 = δ 00 δ 20 , p 10 = δ 10 δ 20 , p 01 = δ 01 δ 20 , and p 11 = δ 00 δ 20 ; furthermore, R 3 is a power series in ( U , V ) with powers U i and V j satisfying i + j 3 , and the coefficients smoothly depend upon λ 1 and λ 2 .
Next, we use the affine transformations v 1 = U + p 10 2 and v 2 = V ; then, system (28) becomes
d v 1 d τ = v 2 , d v 2 d τ = q 00 + q 01 v 2 + v 1 2 + q 11 v 1 v 2 + R 4 ( v 1 , v 2 ) ,
where q 00 = p 00 p 10 2 4 , q 01 = p 01 p 11 p 10 2 , and q 11 = p 11 ; furthermore, R 4 is a power series in ( v 1 , v 2 ) with powers v 1 i and v 2 j satisfying i + j 3 , and the coefficients smoothly depend upon λ 1 and λ 2 .
Making the change of variables w 1 = q 11 2 v 1 , w 2 = q 11 3 v 2 , and t = τ q 11 , system (29) becomes
d w 1 d t = w 2 , d w 2 d t = ζ 1 ( λ 1 , λ 2 ) + ζ 2 ( ( λ 1 , λ 2 ) ) w 2 + w 1 2 + w 1 w 2 + R 5 ( w 1 , w 2 ) ,
where
ζ 1 ( λ 1 , λ 2 ) = q 00 q 11 4 , ζ 2 ( λ 1 , λ 2 ) = q 01 q 11 ;
furthermore, R 5 is a power series in ( w 1 , w 2 ) with powers w 1 i and w 2 j satisfying i + j 3 , and the coefficients smoothly depend upon λ 1 and λ 2 .
Case 2: δ 20 < 0
Consider U ¯ = z 1 , V ¯ = z 2 δ 20 , and d τ = δ 20 d t ; then, system (27) becomes
d U ¯ d τ = V ¯ , d V ¯ d τ = p ¯ 00 + p ¯ 10 U ¯ + p ¯ 01 V ¯ + U ¯ 2 + p ¯ 11 U ¯ V ¯ + R ¯ 3 ( U ¯ , V ¯ ) ,
where p ¯ 00 = δ 00 δ 20 , p ¯ 10 = δ 10 δ 20 , p ¯ 01 = δ 01 δ 20 , and p ¯ 11 = δ 00 δ 20 ; furthermore, R ¯ 3 is a power series in ( U ¯ , V ¯ ) with powers U ¯ i and V ¯ j satisfying i + j 3 , and the coefficients smoothly depend upon λ 1 and λ 2 .
Next, we use the affine transformations v ¯ 1 = U ¯ p ¯ 10 2 and v ¯ 2 = V ¯ ; then, system (32) becomes
d v ¯ 1 d τ = v ¯ 2 , d v ¯ 2 d τ = q ¯ 00 + q ¯ 01 v ¯ 2 + v ¯ 1 2 + q ¯ 11 v ¯ 1 v ¯ 2 + R ¯ 4 ( v ¯ 1 , v ¯ 2 ) ,
where q ¯ 00 = p ¯ 00 + p ¯ 10 2 4 , q ¯ 01 = p ¯ 01 + p ¯ 11 p ¯ 10 2 , and q ¯ 11 = p ¯ 11 ; furthermore, R ¯ 4 is a power series in ( v ¯ 1 , v ¯ 2 ) with powers v ¯ 1 i and v ¯ 2 j satisfying i + j 3 , and the coefficients smoothly depend upon λ 1 and λ 2 .
Consider new variables w ¯ 1 = q ¯ 11 2 v ¯ 1 , w ¯ 2 = q ¯ 11 3 v ¯ 2 , and t = τ q ¯ 11 ; then, system (33) becomes
d w ¯ 1 d t = w ¯ 2 , d w ¯ 2 d t = ζ ¯ 1 ( λ 1 , λ 2 ) + ζ ¯ 2 ( λ 1 , λ 2 ) w ¯ 2 + w ¯ 1 2 + w ¯ 1 w ¯ 2 + R ¯ 5 ( w ¯ 1 , w ¯ 2 ) ,
where
ζ ¯ 1 ( λ 1 , λ 2 ) = q ¯ 00 q ¯ 11 4 , ζ ¯ 2 ( λ 1 , λ 2 ) = q ¯ 01 q ¯ 11 ;
furthermore, R ¯ 5 is a power series in ( w ¯ 1 , w ¯ 2 ) with powers w ¯ 1 i and w ¯ 2 j satisfying i + j 3 , and the coefficients smoothly depend upon λ 1 and λ 2 .
The number of cases can be minimized by using ζ 1 ( λ 1 , λ 2 ) and ζ 2 ( λ 1 , λ 2 ) to denote ζ ¯ 1 ( λ 1 , λ 2 ) and ζ ¯ 2 ( λ 1 , λ 2 ) in (35). If | ( ζ 1 , ζ 2 ) ( λ 1 , λ 2 ) | λ 1 = λ 2 = 0 0 , 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 ( λ 1 , λ 2 ) 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 S N = ( λ 1 , λ 2 ) : ζ 1 ( λ 1 , λ 2 ) = 0 , ζ 2 ( λ 1 , λ 2 ) 0 ;
  • The Hopf bifurcation curve H = ( λ 1 , λ 2 ) : ζ 2 ( λ 1 , λ 2 ) = ± ζ 1 ( λ 1 , λ 2 ) , ζ 1 ( λ 1 , λ 2 ) < 0 ;
  • The homoclinic bifurcation curve H L = ( λ 1 , λ 2 ) : ζ 2 ( λ 1 , λ 2 ) = ± 5 7 ζ 1 ( λ 1 , λ 2 ) , ζ 1 ( λ 1 , λ 2 ) < 0 .

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 β = 0.1 , γ = 0.3 , δ = 0.1 , ρ = 2 , and ξ = 0.2 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
    AxialInterior
    0 < θ < θ S N = 0.576118 E 0 = ( 0 , 0 ) (saddle) E 2 * = ( 0.620942 , 0.720942 ) (stable), E 3 * = ( 0.253568 , 0.353568 ) (saddle)Figure 1a
    θ = θ S N E δ = ( 0 , δ ) (stable) E * = ( 0.382843 , 0.482843 ) (degenerate singularity), emergence of saddle-node bifurcationFigure 1b,d
    θ > θ S N E β = ( β , 0 ) (unstable), E 1 = ( 1 , 0 ) (saddle)NoneFigure 1c
    Figure 1. Existence of equilibrium points. (a) At θ = 0.1 . (b) At θ = θ S N = 0.576118 . (c) At θ = 0.6 . The blue and red (separatrix) curves represent the solution trajectories of the system 5 and the arrows depict the direction of these trajectories in (ac). (d) Saddle-node bifurcation diagram (upper and lower curves represent stable equilibria and unstable equilibria, respectively).
    Figure 1. Existence of equilibrium points. (a) At θ = 0.1 . (b) At θ = θ S N = 0.576118 . (c) At θ = 0.6 . The blue and red (separatrix) curves represent the solution trajectories of the system 5 and the arrows depict the direction of these trajectories in (ac). (d) Saddle-node bifurcation diagram (upper and lower curves represent stable equilibria and unstable equilibria, respectively).
    Symmetry 16 01552 g001
  • 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 θ = 0.4 , β = 0.1 , γ = 0.3 , δ = 0.1 , and ξ = 0.2 ; the obtained results are in Table 2.
  • Fixing the parameters as β = 0.1 , γ = 0.3 , δ = 0.1 , and ξ = 0.2 in (22), we first found the threshold values, which come out to be θ 0 = 0.576118 and ρ 0 = 0.193091 . Corresponding to these parameters, system (5) exhibited a unique equilibrium point E * = ( 0.382843 , 0.482843 ) . By incorporating all the specified parameters, we derived the subsequent system:
    d x d t = x ( 1 x ) ( x 0.1 ) ( 1 + ( 0.576118 + λ 1 ) y ) ( x + 0.3 ) 0.2 x y x + 0.1 , d y d t = ( 0.193091 + λ 2 ) y ( 1 y x + 0.1 ) .
    To transfer this equilibrium point to the origin, we made the transformations u = x 0.382843 and v = y 0.482843 . Then, we made the affine transformations y 1 = u and y 2 = α 10 u + α 01 v . So, the system reduced to
    d y 1 d t = y 2 + + ξ 00 + ξ 20 y 1 2 + ξ 11 y 1 y 2 + ξ 02 y 2 2 + R ˜ 1 ( y 1 , y 2 ) d y 2 d t = η 00 + η 10 y 1 + η 01 y 2 + η 20 y 1 2 + η 11 y 1 y 2 + η 02 y 2 2 + R ˜ 2 ( y 1 , y 2 )
    where ξ 00 = 0.0765685 + 0.097868 1.27817 + 0.482843 λ 1 , ξ 20 = ( 0.100088 0.0369706 λ 1 2 ) 2 , ξ 02 = ( 0.195736 0.0369706 λ 1 2 ) 2 , ξ 11 = ( 0.0666203 0.0369706 λ 1 2 ) 2 , η 00 = ( 3.42781 10 15 + 0.482843 λ 1 ) 2 , η 10 = ( 2.60784 10 13 + 0.482843 λ 1 ) 2 , η 01 = ( 3.9968 10 15 + 0.482843 λ 1 , η 20 = ( 0.024702 0.0369706 λ 1 2 ) 2 , η 02 = ( 0.222134 + 0.0369706 λ 1 2 ) 2 , and η 11 = ( 0.0164421 + 0.0369706 λ 1 2 ) 2 ; furthermore, R ˜ 1 and R ˜ 2 were power series in ( y 1 , y 2 ) with powers y 1 i and y 2 j satisfying i + j 3 , and the coefficients had a smooth dependence upon λ 1 and λ 2 .
    Next, we performed the C transformations x 1 = y 1 and x 2 = y 2 + ξ 00 + ξ 20 y 1 2 + ξ 11 y 1 y 2 + ξ 02 y 2 2 . We made a change in the time variable d t = ( 1 γ 02 x 1 ) and rewrote T as t. Then, using the transformations z 1 = x 1 and z 2 = ( 1 γ 02 x 1 ) x 2 + R ^ 1 ( x 1 , x 2 ) , system (37) is reduced to
    d z 1 d t = z 2 d z 2 d t = δ 00 + δ 10 z 1 + δ 01 z 2 + δ 20 z 1 2 + δ 11 z 1 z 2 + R ¯ 2 ( z 1 , z 2 )
    where δ 00 = ( 8.45997 10 16 + 0.482843 λ 1 ) 2 , δ 10 = ( 4.31962 10 14 + 0.0369706 λ 1 2 ) 2 ) , δ 01 = ( 1.24357 10 15 + 0.0369706 λ 1 2 ) 2 ) , δ 20 = ( 0.00656111 0.0369706 λ 1 2 ) 4 ) , and δ 11 = ( 0.0298712 0.0369706 λ 1 2 ) 4 ) ; furthermore, R ¯ 2 was a power series in ( z 1 , z 2 ) with powers z 1 i and z 2 j satisfying i + j 3 , and the coefficients smoothly depended upon λ 1 and λ 2 .
    When λ 1 and λ 2 were very small, then δ 20 = 0.151939 . Now, we first considered U = z 1 , V = z 2 δ 20 , and d τ = δ 20 d t ; then, we used v 1 = U p 10 2 and v 2 = V . Then, making a change in the variables w 1 = q 11 2 v 1 , w 2 = q 11 3 v 2 , and t = τ q 11 , the system reduced to
    d w 1 d t = w 2 d w 2 d t = ζ 1 ( λ 1 , λ 2 ) + ζ 2 ( λ 1 , λ 2 ) w 2 + w 1 2 + w 1 w 2 + R ¯ 5 ( w 1 , w 2 )
    where
    ζ 1 ( λ 1 , λ 2 ) = q 00 q 11 4 , ζ 2 ( λ 1 , λ 2 ) = q 01 q 11 ,
    q 00 = ( 1.46715 10 19 3.16471 10 7 λ 1 14 λ 2 3 ) 2 , q 01 = ( 6.37002 10 16 0.0369706 λ 1 2 ) 4 ) , and q 11 = ( 0.0298712 0.0369706 λ 1 2 ) 4 ) ; furthermore, R ¯ 5 was a power series in ( w 1 , w 2 ) with powers w 1 i and w 2 j that satisfied i + j 3 , and the coefficients smoothly depended upon λ 1 and λ 2 .
The determinant E = | ( ζ 1 , ζ 2 ) ( λ 1 , λ 2 ) | λ 1 = λ 2 = 0 = 65.9058 0 ; 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 ( 0 , 0 ) , the three bifurcation curves divided the λ 1 λ 2 -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 λ 1 = λ 2 = 0 , 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 λ 1 and λ 2 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 λ 1 and λ 2 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 λ 1 and λ 2 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 λ 1 and λ 2 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.

6. Impact of Fear Effect

The emphasis of this section is on exploring the influence that fear and the double Allee effect have on prey species. The impact was studied in two parts: the first part dealt with a scenario where the predator’s population was constant, and the second part dealt with a scenario when the predator’s population was governed by the second differential equation of model (4).
  • Consider the prey equation of model (4) while the predator was taken as a constant. In this scenario, model (4) is reduced to
    d N d T = r N 1 + f P 1 N K N m N + n α N P a + N .
    In reality, this scenario where the predators remain constant in an environment is highly unlikely, but there do exist some cases in which keeping a species’ population constant is possible, and this allows us to control the dynamics of an ecological system to some extent. To fulfill the aim of this study, the predator species was considered constant so that the impact of the fear effect could be explored at great length. Since fear is a feeling that is induced in prey due to the presence of predators only, we cannot disregard the existence of predators even while doing a study on just the growth function of the prey.
    In addition to the fear effect, we also considered the influence of the double Allee effect on the growth rate of prey species. This purpose could be served by plotting the per capita growth rate of the prey, i.e., 1 N d N d T against the population of the prey N (see Figure 4). From Figure 4a, we can see that the incorporation of phenomena like the double Allee effect and the fear effect led to a decline in the growth rate of the prey species. Figure 4b–d reflects that as the values of the parameters m, n, and f increased, respectively, the per capita growth rate of prey species decreased. The prey population possessed a threshold value for these parameters, which was crucial in determining the extinction or survival of the species.
    Consider N = K x , P = K y , and T = t r ; then, system (41) reduces to the following non-dimensional form:
    d x d t = x ( 1 x ) ( x β ) ( 1 + θ y ) ( x + γ ) ξ x y x + δ = f ( x ) ,
    where β = m K , γ = n K , δ = a K , θ = f K , and ξ = α r .
    Except x = 0 , system (42) has non-trivial equilibrium points given by
    x 3 J x 2 K x + L = 0 ,
    where J = 1 + β δ , K = β δ + δ β ξ y ξ θ y 2 , and L = β δ + γ ( ξ y + ξ θ y 2 ) . If any one of the conditions J > 0 , K 0 ; J > 0 , K 0 ; or J 0 , K 0 holds, then by Descartes’ sign rule, it is clear that (43) has a negative root. Let this root be ν , ν > 0 . Dividing (43) by ( x + ν ) , we obtain the following quadratic equation:
    x 2 ( J + ν ) x + ν ( J + ν ) K = 0 ,
    Taking Δ 1 = ( J + ν ) 2 4 ( ν ( J + ν ) K ) , the ensuing outcome is described in the following theorem.
    Theorem 7.
    If any of the conditions ( i ) J > 0 , K 0 , ( i i ) J > 0 , K 0 , or ( i i i ) J 0 , K 0 hold, then system (42) has the following:
    (a) 
    Two equilibrium points E ¯ 1 = x ¯ 1 and E ¯ 2 = x ¯ 2 when Δ 1 > 0 , where x ¯ 1 = J + ν Δ 1 2 and x ¯ 2 = J + ν + Δ 1 2 ;
    (b) 
    One equilibrium point E ¯ 3 = x ¯ 3 when Δ 1 = 0 , where x ¯ 3 = J + ν 2 ;
    (c) 
    No equilibrium point when Δ 1 < 0 .
    The stability of system (42) can be determined by the sign of f ( x ) at various equilibrium points [46]. Differentiating f ( x ) given in Equation (42), we obtain
    f ( x ) = x ( 1 2 x + β ) ( 1 + θ y ) ( x + γ ) x ( 1 x ) ( x β ) ( 1 + θ y ) ( x + γ ) 2 + ξ x y ( x + δ ) 2 + ( 1 x ) ( x β ) ( 1 + θ y ) ( x + γ ) ξ y x + δ
    The theorem stated below and Figure 5 show the stability of system (42).
    Theorem 8.
    (a)  The trivial equilibrium point is always stable, as f ( x ) < 0 .
    (b) 
    The equilibrium points E ¯ 1 , E ¯ 2 , and E ¯ 3 , if they exist, are unstable when x ¯ 1 ( 1 2 x ¯ 1 + β ) ( 1 + θ y ) ( x ¯ 1 + γ ) + ξ x ¯ 1 y ( γ δ ) ( x ¯ 1 + γ ) ( x ¯ 1 + δ ) 2 > 0 , stable when x ¯ 2 ( 1 2 x ¯ 2 + β ) ( 1 + θ y ) ( x ¯ 2 + γ ) + ξ x ¯ 2 y ( γ δ ) ( x ¯ 2 + γ ) ( x ¯ 2 + δ ) 2 < 0 , and semi-stable (saddle-node bifurcation), respectively.
  • Now, coming back to the original system (5), we present the influence of fear and the double Allee effect. Figure 6a–c illustrate that as the Allee effect parameters ( β and γ ) and fear effect parameter ( θ ), respectively, were increased, the equilibrium point E 2 * moved downward. This implies that an increase in the double Allee effect parameters and the fear effect led to the existence of an equilibrium point at a lower population. As soon as the critical value of these parameters was crossed, the system moved toward the prey-free equilibrium, which is stable in nature.
    Now, to observe the collective influence of the fear effect and the double Allee effect, 3D phase portraits were used. Figure 7a shows the behavior of the prey population at the stable equilibrium point E 2 * against the parameters θ and β , and it can be seen that the prey population decreased as the values of the parameters θ and β increased. Similarly, Figure 7b reflects the same result with respect to the parameters θ and γ . Now, we know that y 2 * = x 2 * + δ , which implies the predator population will also behave in the same way as the prey population regarding the parameters θ and β , as well as θ and γ . From this, we could conclude that the higher the fear and value of the Allee effect parameters, the lower the population where system (5) attained stability. Ecologically, one can infer that escalation in the fear effect leads to a decline in prey and predator populations. As this fear intensifies, the prey population indulges less in foraging, and their interaction with each other decreases, which ultimately leads to a lower reproduction rate and, consequently, a reduced overall population. The diminishing predator–prey interaction is attributed to a declining prey population, which results in a decline in the population of predators. The populations of both species drop to the point where prey go extinct and predators survive (due to the presence of additional food).

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 θ = θ S N , 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 E 2 * 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.

Author Contributions

Conceptualization: M.K.S. and A.S.; Methodology: M.K.S., A.S. and L.M.S.-R.; Formal Analysis: M.K.S. and A.S.; Investigation: M.K.S., A.S. and L.M.S.-R.; Writing—Review and Editing: M.K.S., A.S. and L.M.S.-R.; Visualization: M.K.S. and A.S.; Supervision: M.K.S. and L.M.S.-R. All authors read and agreed to the published version of this manuscript.

Funding

This research received no external funding.

Data Availability Statement

Data are contained within the article.

Acknowledgments

The author Arushi Sharma expresses gratitude to Banasthali Vidyapith for providing the opportunity and inspiration to compose this research paper.

Conflicts of Interest

The authors declare there is no conflicts of interest.

References

  1. Lotka, A.J. Elements of Physical Biology; Williams & Wilkins: Baltimore, MD, USA, 1925. [Google Scholar]
  2. Volterra, V. Fluctuations in the abundance of a species considered mathematically. Nature 1926, 118, 558–560. [Google Scholar] [CrossRef]
  3. Leslie, P.; Gower, J. The properties of a stochastic model for the predator-prey type of interaction between two species. Biometrika 1960, 47, 219–234. [Google Scholar] [CrossRef]
  4. Aziz-Alaoui, M.; Okiye, M.D. Boundedness and global stability for a predator-prey model with modified Leslie-Gower and Holling-type II schemes. Appl. Math. Lett. 2003, 16, 1069–1075. [Google Scholar] [CrossRef]
  5. Zhu, Y.; Wang, K. Existence and global attractivity of positive periodic solutions for a predator–prey model with modified Leslie–Gower Holling-type II schemes. J. Math. Anal. Appl. 2011, 384, 400–408. [Google Scholar] [CrossRef]
  6. Ji, C.; Jiang, D.; Shi, N. Analysis of a predator–prey model with modified Leslie–Gower and Holling-type II schemes with stochastic perturbation. J. Math. Anal. Appl. 2009, 359, 482–498. [Google Scholar] [CrossRef]
  7. Jiang, G.; Lu, Q.; Qian, L. Complex dynamics of a Holling type II prey–predator system with state feedback control. Chaos Solitons Fractals 2007, 31, 448–461. [Google Scholar] [CrossRef]
  8. 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]
  9. Nindjin, A.; Aziz-Alaoui, M.; Cadivel, M. Analysis of a predator–prey model with modified Leslie–Gower and Holling-type II schemes with time delay. Nonlinear Anal. Real World Appl. 2006, 7, 1104–1118. [Google Scholar] [CrossRef]
  10. Song, X.; Li, Y. Dynamic behaviors of the periodic predator–prey model with modified Leslie-Gower Holling-type II schemes and impulsive effect. Nonlinear Anal. Real World Appl. 2008, 9, 64–79. [Google Scholar] [CrossRef]
  11. Gupta, R.; Chandra, P. Bifurcation analysis of modified Leslie–Gower predator–prey model with Michaelis–Menten type prey harvesting. J. Math. Anal. Appl. 2013, 398, 278–295. [Google Scholar] [CrossRef]
  12. Allee, W.C. Co-operation among animals. Am. J. Sociol. 1931, 37, 386–398. [Google Scholar] [CrossRef]
  13. Dennis, B. Allee effects: Population growth, critical density, and the chance of extinction. Nat. Resour. Model. 1989, 3, 481–538. [Google Scholar] [CrossRef]
  14. Lewis, M.A.; Kareiva, P. Allee dynamics and the spread of invading organisms. Theor. Popul. Biol. 1993, 43, 141–158. [Google Scholar] [CrossRef]
  15. Courchamp, F.; Berec, L.; Gascoigne, J. Allee Effects in Ecology and Conservation; OUP Oxford: New York, NY, USA, 2008. [Google Scholar]
  16. Angulo, E.; Roemer, G.W.; Berec, L.; Gascoigne, J.; Courchamp, F. Double Allee effects and extinction in the island fox. Conserv. Biol. 2007, 21, 1082–1091. [Google Scholar] [CrossRef] [PubMed]
  17. Berec, L.; Angulo, E.; Courchamp, F. Multiple Allee effects and population management. Trends Ecol. Evol. 2007, 22, 185–191. [Google Scholar] [CrossRef]
  18. González-Olivares, E.; Mena-Lorca, J.; Rojas-Palma, A.; Flores, J.D. Dynamical complexities in the Leslie–Gower predator–prey model as consequences of the Allee effect on prey. Appl. Math. Model. 2011, 35, 366–381. [Google Scholar] [CrossRef]
  19. Cai, Y.; Zhao, C.; Wang, W.; Wang, J. Dynamics of a Leslie–Gower predator–prey model with additive Allee effect. Appl. Math. Model. 2015, 39, 2092–2106. [Google Scholar] [CrossRef]
  20. Feng, P.; Kang, Y. Dynamics of a modified Leslie–Gower model with double Allee effects. Nonlinear Dyn. 2015, 80, 1051–1062. [Google Scholar] [CrossRef]
  21. 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]
  22. Singh, M.K.; Bhadauria, B.; 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]
  23. 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. [Google Scholar] [CrossRef]
  24. Wang, F.; Yang, R.; Zhang, X. Turing patterns in a predator–prey model with double Allee effect. Math. Comput. Simul. 2024, 220, 170–191. [Google Scholar] [CrossRef]
  25. Yin, W.; Li, Z.; Chen, F.; He, M. Modeling Allee Effect in the Leslie-Gower Predator–Prey System Incorporating a Prey Refuge. Int. J. Bifurc. Chaos 2022, 32, 2250086. [Google Scholar] [CrossRef]
  26. Cruz, R.L. Stability of a Leslie-Gower type predator-prey model with a strong Allee effect with delay. Sel. Mat. 2022, 9, 24–33. [Google Scholar] [CrossRef]
  27. Wang, F.; Yang, R. Dynamics of a delayed reaction–diffusion predator–prey model with nonlocal competition and double Allee effect in prey. Int. J. Biomath. 2023, 2350097. [Google Scholar] [CrossRef]
  28. Liu, Y.; Zhang, Z.; Li, Z. The Impact of Allee Effect on a Leslie–Gower Predator–Prey Model with Hunting Cooperation. Qual. Theory Dyn. Syst. 2024, 23, 88. [Google Scholar] [CrossRef]
  29. Lima, S.L. Nonlethal effects in the ecology of predator-prey interactions. Bioscience 1998, 48, 25–34. [Google Scholar] [CrossRef]
  30. Creel, S.; Christianson, D. Relationships between direct predation and risk effects. Trends Ecol. Evol. 2008, 23, 194–201. [Google Scholar] [CrossRef]
  31. Cresswell, W. Predation in bird populations. J. Ornithol. 2011, 152, 251–263. [Google Scholar] [CrossRef]
  32. Preisser, E.L.; Bolnick, D.I. The many faces of fear: Comparing the pathways and impacts of nonconsumptive predator effects on prey populations. PLoS ONE 2008, 3, e2465. [Google Scholar] [CrossRef]
  33. Wang, X.; Zanette, L.; Zou, X. Modelling the fear effect in predator–prey interactions. J. Math. Biol. 2016, 73, 1179–1204. [Google Scholar] [CrossRef] [PubMed]
  34. Sarkar, K.; Khajanchi, S. Impact of fear effect on the growth of prey in a predator-prey interaction model. Ecol. Complex. 2020, 42, 100826. [Google Scholar] [CrossRef]
  35. Li, Y.; He, M.; Li, Z. Dynamics of a ratio-dependent Leslie–Gower predator–prey model with Allee effect and fear effect. Math. Comput. Simul. 2022, 201, 417–439. [Google Scholar] [CrossRef]
  36. Chen, M.; Takeuchi, Y.; Zhang, J.F. Dynamic complexity of a modified Leslie–Gower predator–prey system with fear effect. Commun. Nonlinear Sci. Numer. Simul. 2023, 119, 107109. [Google Scholar] [CrossRef]
  37. Wu, H.; Li, Z.; He, M. Dynamic analysis of a Leslie-Gower predator-prey model with the fear effect and nonlinear harvesting. Math. Biosci. Eng. 2023, 20, 18592–18629. [Google Scholar] [CrossRef]
  38. Al-Momen, S.; Naji, R.K. The dynamics of modified Leslie-Gower predator-prey model under the influence of nonlinear harvesting and fear effect. Iraqi J. Sci 2022, 63, 259–282. [Google Scholar] [CrossRef]
  39. Vinoth, S.; Vadivel, R.; Hu, N.T.; Chen, C.S.; Gunasekaran, N. Bifurcation Analysis in a Harvested Modified Leslie–Gower Model Incorporated with the Fear Factor and Prey Refuge. Mathematics 2023, 11, 3118. [Google Scholar] [CrossRef]
  40. Jamil, A.R.M.; Naji, R.K. Modeling and Analysis of the Influence of Fear on the Harvested Modified Leslie–Gower Model Involving Nonlinear Prey Refuge. Mathematics 2022, 10, 2857. [Google Scholar] [CrossRef]
  41. Halder, S.; Bhattacharyya, J.; Pal, S. Comparative studies on a predator–prey model subjected to fear and Allee effect with type I and type II foraging. J. Appl. Math. Comput. 2020, 62, 93–118. [Google Scholar] [CrossRef]
  42. Purnomo, A.S.; Darti, I.; Suryanto, A.; Kusumawinahyu, W.M. Fear Effect on a Modified Leslie-Gower Predator-Prey Model with Disease Transmission in Prey Population. Eng. Lett. 2023, 31. [Google Scholar]
  43. Qi, H.; Meng, X. Dynamics of a stochastic predator-prey model with fear effect and hunting cooperation. J. Appl. Math. Comput. 2023, 69, 2077–2103. [Google Scholar] [CrossRef]
  44. Pal, D.; Kesh, D.; Mukherjee, D. Cross-diffusion mediated Spatiotemporal patterns in a predator–prey system with hunting cooperation and fear effect. Math. Comput. Simul. 2024, 220, 128–147. [Google Scholar] [CrossRef]
  45. Perko, L. Differential Equations and Dynamical Systems; Springer Science & Business Media: New York, NY, USA, 2013; Volume 7. [Google Scholar]
  46. Layek, G. An Introduction to Dynamical Systems and Chaos; Springer: New Delhi, India, 2015; Volume 449. [Google Scholar]
Figure 2. Phase potraits of system 5. The blue, red (separatrix) and green (limit cycle) curves represent the solution trajectories of the system 5 and the arrows depict the direction of these trajectories. (a) At ρ = 0.15 , E 3 * was consistently a saddle and E 2 * was stable. (b) At ρ = 0.1031 , E 2 * was enclosed within a homoclinic loop. (c) At ρ = ρ h f = 0.0845887 , an unstable limit cycle emerged at point E 2 * . (d) At ρ = 0.01 , E 2 * was unstable.
Figure 2. Phase potraits of system 5. The blue, red (separatrix) and green (limit cycle) curves represent the solution trajectories of the system 5 and the arrows depict the direction of these trajectories. (a) At ρ = 0.15 , E 3 * was consistently a saddle and E 2 * was stable. (b) At ρ = 0.1031 , E 2 * was enclosed within a homoclinic loop. (c) At ρ = ρ h f = 0.0845887 , an unstable limit cycle emerged at point E 2 * . (d) At ρ = 0.01 , E 2 * was unstable.
Symmetry 16 01552 g002
Figure 3. (a) Bifurcation diagram. (b) Cusp at λ 1 = λ 2 = 0 . (c) Region 1: no interior equilibrium point. (d) Region 2: two positive equilibrium points (one exhibiting instability and the other characterize as a saddle point). (e) Region 3: two positive equilibrium points (unstable limit cycle that enclosed one interior point and a saddle point). (f) Region 4: two positive equilibrium points (one that exhibited stability and the other characterized as a saddle point). The blue curves represent the solution trajectories of the system 5 and the arrows depict the direction of these trajectories in (bf).
Figure 3. (a) Bifurcation diagram. (b) Cusp at λ 1 = λ 2 = 0 . (c) Region 1: no interior equilibrium point. (d) Region 2: two positive equilibrium points (one exhibiting instability and the other characterize as a saddle point). (e) Region 3: two positive equilibrium points (unstable limit cycle that enclosed one interior point and a saddle point). (f) Region 4: two positive equilibrium points (one that exhibited stability and the other characterized as a saddle point). The blue curves represent the solution trajectories of the system 5 and the arrows depict the direction of these trajectories in (bf).
Symmetry 16 01552 g003
Figure 4. (a) Population growth rate of prey species with constant predators when governed by logistic growth only, by double Allee effect, by fear effect, and by both double Allee and fear effect. (b) Effect of m on growth curve of system (41). (c) Effect of n on growth curve of system (41). (d) Effect of f on growth curve of system (41).
Figure 4. (a) Population growth rate of prey species with constant predators when governed by logistic growth only, by double Allee effect, by fear effect, and by both double Allee and fear effect. (b) Effect of m on growth curve of system (41). (c) Effect of n on growth curve of system (41). (d) Effect of f on growth curve of system (41).
Symmetry 16 01552 g004
Figure 5. Behavior of the solution for system (42), where blue curves represent the constant solution of system (42) and green curves represent the non-constant solutions of system (42).
Figure 5. Behavior of the solution for system (42), where blue curves represent the constant solution of system (42) and green curves represent the non-constant solutions of system (42).
Symmetry 16 01552 g005
Figure 6. (a) Impact of β on E 2 * . (b) Impact of γ on E 2 * . (c) Impact of θ on E 2 * .
Figure 6. (a) Impact of β on E 2 * . (b) Impact of γ on E 2 * . (c) Impact of θ on E 2 * .
Symmetry 16 01552 g006
Figure 7. (a) Effects of θ and β on x 2 * . (b) Effects of θ and γ on x 2 * .
Figure 7. (a) Effects of θ and β on x 2 * . (b) Effects of θ and γ on x 2 * .
Symmetry 16 01552 g007
Table 2. Number and characteristics of equilibrium points that resulted from the emergence of a Hopf bifurcation.
Table 2. Number and characteristics of equilibrium points that resulted from the emergence of a Hopf bifurcation.
ρ Equilibrium Points (Stability)Figure
AxialInterior
ρ = 0.15 E 0 = ( 0 , 0 ) (saddle) E 2 *  =  ( 0.507277 , 0.607277 ) (stable), E 3 *  =  ( 0.296426 , 0.396426 ) (saddle)Figure 2a
ρ = 0.1031 E δ = ( 0 , δ ) (stable)Homoclinic loop enclosing E 2 *  =  ( 0.507277 , 0.607277 ) (stable), E 3 *  =  ( 0.296426 , 0.396426 ) (saddle)Figure 2b
ρ = ρ h f = 0.0845887 E β = ( β , 0 ) (unstable)Unstable limit cycle (as Λ = 0.66667 > 0 ) around E 2 *  =  ( 0.507277 , 0.607277 ) (stable), E 3 *  =  ( 0.296426 , 0.396426 ) (saddle)Figure 2c
ρ = 0.01 E 1 = ( 1 , 0 ) (saddle) E 2 *  =  ( 0.507277 , 0.607277 ) (unstable), E 3 *  =  ( 0.296426 , 0.396426 ) (saddle)Figure 2d
Table 3. Behavior of system in various regions.
Table 3. Behavior of system in various regions.
RegionNumber of Interior Equilibrium PointsDynamic Behavior Exhibited by RegionsFigure
OriginOneCuspFigure 3b
Region 1Zero Figure 3c
Region 2TwoOne was consistently a saddle while the other was unstableFigure 3d
Region 3TwoOne was consistently a saddle while the other was enclosed by an unstable limit cycleFigure 3e
Region 4TwoOne was consistently a saddle while the other was stableFigure 3f
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Singh, M.K.; Sharma, A.; Sánchez-Ruiz, L.M. Dynamical Complexity of Modified Leslie–Gower Predator–Prey Model Incorporating Double Allee Effect and Fear Effect. Symmetry 2024, 16, 1552. https://doi.org/10.3390/sym16111552

AMA Style

Singh MK, Sharma A, Sánchez-Ruiz LM. Dynamical Complexity of Modified Leslie–Gower Predator–Prey Model Incorporating Double Allee Effect and Fear Effect. Symmetry. 2024; 16(11):1552. https://doi.org/10.3390/sym16111552

Chicago/Turabian Style

Singh, Manoj Kumar, Arushi Sharma, and Luis M. Sánchez-Ruiz. 2024. "Dynamical Complexity of Modified Leslie–Gower Predator–Prey Model Incorporating Double Allee Effect and Fear Effect" Symmetry 16, no. 11: 1552. https://doi.org/10.3390/sym16111552

APA Style

Singh, M. K., Sharma, A., & Sánchez-Ruiz, L. M. (2024). Dynamical Complexity of Modified Leslie–Gower Predator–Prey Model Incorporating Double Allee Effect and Fear Effect. Symmetry, 16(11), 1552. https://doi.org/10.3390/sym16111552

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop