A Discrete-Time Modified Leslie-Gower Model with Double Allee Effect and Its Stability and Bifurcation Analysis

Khan, M.; Abbas, Mujahid; Haq, Absar; Nazeer, Waqas

doi:10.61091/um119-10

Abstract

1. Introduction
2. The Fixed Points & Their Stability
3. Bifurcations
4. Numerical Simulations
5. Conclusion
Conflict of Interest

References

Utilitas Mathematica

Volume 119
Pages: 117-133

Research article

A Discrete-Time Modified Leslie-Gower Model with Double Allee Effect and Its Stability and Bifurcation Analysis

^¹,², ^¹,³, ^⁴, ^¹

¹Department of Mathematics, Government College University, Lahore 54000, Pakistan

²Department of Mathematics, Riphah International University-Lahore Campus, Islamabad, Pakistan

³Department of Medical Research, China Medical University, Taichung 40402, Taiwan

⁴Department of Basic Sciences and Humanities, University of Engineering and Technology, Lahore(NWL Campus), Pakistan

Received: 11/12/2023
Accepted: 27/02/2024
Published: 31/03/2024

Copyright Link
License

Abstract

Generally, all the models discussed so far are continuous time models. The continuous time models are quite apt at explaining the phenomena they are trying to predict and have known methods to get information from these type of models. But these models are not accurate for the physical systems which are observed over discreet time periods or which have non-continuous phenomena embedded in them, like production of new generation. Some species like salmon have non-overlapping generation characteristics since they have an annual spawning season and are born each year at a certain time. The discrete models are much more apt in describing the nature’s complex dynamics than the continuous models. A discrete-time modified Leslie-Gower system with double Allee effect is studied in this paper. The stability analysis of interior fixed points is performed. Using center manifold theorem it is shown that the system under consideration exhibits period-doubling and Neimark-Sacker bifurcations. The numerical simulations are provided to illustrate the consistency of the theoretical results.

Keywords: Double Allee effect, Bifurcation, Leslie-Gower Model, Stability

1. Introduction

The relationships between different species on Earth can be modeled as predator-prey interactions, which are fundamental to the survival of many organisms [1, 2]. Mathematical models can be used to approximate these dynamics. However, predator-prey models are not only useful for understanding the dynamics of predators and their prey but also for managing renewable resources [3]. While Lotka [4] and Volterra [5] were the first to model predator-prey populations, the Lotka-Volterra model was not realistic enough as it neglected many important aspects of predator-prey interactions. Thus, various modifications have been proposed by researchers to address these limitations [6].

Leslie [7] proposed an alternative to the Lotka-Volterra model that considered the carrying capacity of the predator’s environment and the number of prey as proportional. This model also accounted for the fact that predator and prey populations must have bounded increasing capacity. This model has a unique positive fixed point that has been shown to be globally asymptotically stable for any allowed parameters [8, 9, 10]. In [11, 12], the authors studied the uniqueness of limit cycles and the existence of Hopf bifurcation for this model.

However, this model still has limitations as it does not account for the fact that predators can switch between different prey depending on the conditions and needs of the environment [13]. To address this limitation, a modified Leslie-Gower predator-prey model was proposed [14]. This model has been used to model the dynamics of various systems such as prickly-pear cactus [15] and mite outbreaks in fruit trees [16, 17]. For further details on where this model has been applied and studied, refer to [18, 19] and the references therein.

Warder Clyde Allee is credited with introducing the concept of Allee effects, which are believed to be widespread in nature despite being difficult to detect [20]. Allee effects are characterized by a positive correlation between population fitness and size over a finite interval, and can lead to a critical population size below which the population cannot persist [21]. Strong Allee effects exhibit a threshold below which the population is driven to extinction, while weak Allee effects do not [22]. Allee effects have been observed in diverse natural phenomena [23], and various mathematical models have been proposed to describe them, some of which are topologically complex [24]. Researchers have explored the bifurcations of predator-prey systems subject to Allee effect [25, 26], and have found that many Allee effects act simultaneously on a single population, particularly in renewable resources [27].

Double Allee effects occur when two mechanisms act together on a single population, and have been observed in plants, vertebrates, and invertebrates [27]. The double Allee effect has also been observed in marine ecosystems [28]. Gonzlez-Olivares et al. investigated the growth of prey subject to double Allee effect using the Lotka-Volterra predator-prey model, and found that regardless of the state of Allee effect, two limit cycles exist [29]. The modified Rosenzweig-MacArthur model was used in [30] to study the two Allee effects, and the authors demonstrated that the positive fixed point is locally asymptotically stable for prey, and Hopf bifurcation may occur, generating a stable limit cycle. Researchers have also explored the effects of double Allee effect on prey population using some ratio-dependent predator-prey model, discussing the stability of equilibrium points and their bifurcations [31, 32, 33].

Continuous time models are commonly used to explain Allee effects, but they may not be accurate for physical systems observed over discrete time periods, or systems with non-continuous phenomena, such as the production of new generations. For instance, species like salmon have non-overlapping generation characteristics since they have an annual spawning season and are born each year at a certain time. Discrete models are therefore better suited to describing the complex dynamics of nature than continuous models [34, 35]. Discrete systems also lend themselves more readily to analytical solutions compared to continuous models [36, 37]. For further details, interested readers may refer to [38],[39] and the references therein.

In this paper, we will use modified Leslie-Gower predator-prey model with double Allee effect on prey, proposed by Singh et al. in [40]. The model is given in [40] as; $\begin{matrix} (1) & {\begin{cases} \frac{d x}{d t} = \frac{x (1 - x) (x - β)}{θ + x} - ξ \frac{x y}{γ + x}, \\ \frac{d y}{d t} = ρ y (1 - \frac{y}{γ + x}) . \end{cases} \end{matrix}$

The initial conditions are assumed to be positive, that is, $x (0) > 0$ and $y (0) > 0$ . Details of the other parameters are provided in [40]. For our purpose, all parameters are positive except for $β$ . A positive value of $β$ corresponds to a strong Allee effect, while a negative value of $β$ corresponds to a weak Allee effect. Singh $et. al$ demonstrated that the system (1) is both bounded and positive. They studied the stability dynamics of the equilibrium points and proved the existence of bifurcations, including fold, Hopf, and Bogdanov-Takens bifurcations. However, the model possesses much richer dynamics than previously established. This can be accomplished by examining the discrete version of system (1). Discretization is useful for ecosystems in which consecutive generations do not overlap. We utilize the forward Euler method to discretize system (1). The discrete version of system (1) with a step-size of $h$ is defined as follows: $\begin{matrix} (2) & {\begin{cases} x_{n + 1} = x_{n} + h x_{n} (\frac{(1 - x_{n}) (x_{n} - β)}{θ + x_{n}} - ξ \frac{y_{n}}{γ + x_{n}}), \\ y_{n + 1} = y_{n} + h ρ y_{n} (1 - \frac{y_{n}}{γ + x_{n}}) . \end{cases} \end{matrix}$

The paper focuses on studying a discrete-time modified Leslie-Gower system with double Allee effect, which is more suitable for describing complex dynamics in nature compared to continuous time models. The authors highlight that continuous time models are not accurate for physical systems observed over discreet time periods or those with non-continuous phenomena embedded in them, such as the annual spawning season of salmon. To address this issue, the paper performs stability analysis of interior fixed points and shows that the system exhibits period-doubling and Neimark-Sacker bifurcations using the center manifold theorem. Numerical simulations are also provided to demonstrate the consistency of the theoretical results. By studying this discrete-time model, the authors aim to better understand and predict the dynamics of natural systems.

In this paper, we investigate the dynamics of system (2) and aim to demonstrate the richness of its local dynamics, as well as the conditions under which the positive interior stationary points become non-hyperbolic. The non-hyperbolic point leads to period-doubling bifurcation as well as Neimark-Sacker bifurcation, at a fixed step-size. Additionally, we provide numerical examples to illustrate our findings. This article is structured as follows: in Section 2, we examine the stability of the stationary points. In Section 3, we discuss the existence and conditions for period-doubling and Neimark-Sacker bifurcations. In Section 4, we present numerical examples for both strong and weak Allee effects, accompanied by diagrams. Finally, we conclude our findings in Section 5.

2. The Fixed Points & Their Stability

The following system of simultaneous equations is solved to find the fixed points. $\begin{matrix} (3) & {\begin{cases} x (\frac{(1 - x) (x - β)}{θ + x} - ξ \frac{y}{γ + x}) = 0, \\ ρ y (1 - \frac{y}{γ + x}) = 0. \end{cases} \end{matrix}$

If we define $A = \frac{1 + β - ξ}{2}$ and $B = β + θ ξ$ , then for any scenario, the system has four stationary points on the boundary ${E_{1}, E_{2}, E_{3}, E_{4}} = {(0, 0), (0, γ), (1, 0), (β, 0)}$ . Let $E_{5^{*}} = (A, A + γ)$ , $E_{5} = ({\bar{x}}_{1}, {\bar{x}}_{1} + γ)$ and $E_{6} = ({\bar{x}}_{2}, {\bar{x}}_{2} + γ)$ , where, ${\bar{x}}_{1} = A + \sqrt{A^{2} - B}$ and ${\bar{x}}_{2} = A - \sqrt{A^{2} - B}$ . For positive fixed points, following scenarios are present.

(i) The system has no interior fixed point, if,

β \in (1 + ξ - 2 \sqrt{ξ (1 + θ)}, 1 + ξ + 2 \sqrt{ξ (1 + θ)}) .

(ii) If

β = 1 + ξ \pm 2 \sqrt{ξ (1 + θ)}

, then the system has unique interior equilibrium point,

E_{5^{*}}

. It will always exist for strong Allee effect and it will exist for weak Allee effect if

(1 - ξ)^{2} < 4 ξ θ

(iii) If

β < 1 + ξ - 2 \sqrt{ξ (1 + θ)}

β > 1 + ξ + 2 \sqrt{ξ (1 + θ)}

, i-e,

A^{2} - B > 0

, then

A < \sqrt{A^{2} - B}

θ ξ < - β

, which can only be true in the weak Allee effect. In this case, the only positive stationary point will be

E_{5}

On the other hand, $A > \sqrt{A^{2} - B} > 0$ only if $1 + β > ξ$ and $β + θ ξ > 0$ . In strong Allee effect it is always true. In the weak Allee effect, this is true if $1 > ξ - β > 0$ and $θ ξ > - β$ . In either case the system incorporates two positive fixed points $E_{5}$ and $E_{6}$ . All the scenarios can be observed in Figures 1a and 1b, where change in parametric values contributes to the number of positive fixed points the system may possess.

Figure 1:a) Nullclines for System (2) With Strong Allee Effect, Showing the Existence of Two, One and No Interior Equilibrium Points for Different Parametric Values, (b)) Nullclines for System (2) With Weak Allee Effect, Showing the Existence of Two, One and No Interior Equilibrium Points for Different Parametric Values

At fixed points, the respective jacobian matrices are given by $\begin{aligned} J_{1} & = [\begin{array}{cc} 1 - h \frac{β}{θ} & 0 \\ 0 & 1 + h ρ \end{array}], & J_{2} & = [\begin{array}{cc} 1 - h (ξ - \frac{β}{θ}) & 0 \\ h ρ & 1 - h ρ \end{array}], \\ J_{3} & = [\begin{array}{cc} 1 + h \frac{β - 1}{1 + θ} & - h \frac{ξ}{1 + γ} \\ 0 & 1 + h ρ \end{array}], & J_{4} & = [\begin{array}{cc} 1 - h \frac{β (β - 1)}{β + θ} & - h \frac{β ξ}{β + γ} \\ 0 & 1 + h ρ \end{array}], \\ J_{5} & = [\begin{array}{cc} 1 + h {\bar{x}}_{1} (D - C) & - h {\bar{x}}_{1} D \\ h ρ & 1 - h ρ \end{array}], & J_{6} & = [\begin{array}{cc} 1 + h {\bar{x}}_{2} (D - C) & - h {\bar{x}}_{2} D \\ h ρ & 1 - h ρ \end{array}], \end{aligned}$ where $C = - \frac{2 (A - {\bar{x}}_{i})}{{\bar{x}}_{i} + θ}$ and $D = \frac{ξ}{{\bar{x}}_{i} + γ}$ , $i = 1, 2$ . Note that, for $E_{5}$ , $C > 0$ and $D > 0$ whereas for $E_{6}$ , $- C > 0$ and $D > 0$ , as long as $A^{2} - B > 0$ . In order to find the stability of these fixed points, we will use the following two lemmas.

Lemma 1. Let $t_{i} = T r (J_{i})$ , $d e t_{i} = | J_{i} |$ and $p_{i} (z) = z^{2} - t_{i} z + d e t_{i}$ , $i = 1, \dots, 6$ . Suppose $z_{1}$ , $z_{2}$ are the roots of $p (z)$ . Then:

$| z_{1, 2} | < 1$ if and only if $p (- 1) > 0$ and $d e t < 1$ .
$| z_{1, 2} | > 1$ if and only if $p (- 1) > 0$ and $d e t > 1$ .
$| z_{1} | < 1$ and $| z_{2} | > 1$ or $| z_{1} | > 1$ and $| z_{2} | < 1$ if and only if $p (- 1) < 0$ .
$z_{1} = - 1$ with $| z_{2} | \neq 1$ iff $p (- 1) = 0$ and $t \neq - 2, 0$ .
$z_{1, 2}$ are complex and $| z_{1, 2} | = 1$ if and only if $4 d e t - t^{2} > 0$ and $d e t = 1$ .

Lemma 2. Let $z_{1}$ and $z_{2}$ to be the eigenvalues of $2 \times 2$ Jacobian matrix and $E$ be the positive fixed point. Then $E$ is called

(i)sink if $| z_{1} | < 1$ and $| z_{2} | < 1$ , so sink is locally asymptotically stable.

(ii)source if $| z_{1} | > 1$ and $| z_{2} | > 1$ , so source is locally unstable.

(iii)saddle if $| z_{1} | < 1$ and $| z_{2} | > 1$ (or $| z_{1} | > 1$ and $| z_{2} | < 1$ ).

(iv)non-hyperbolic if either $| z_{1} | = 1$ or $| z_{2} | = 1$ .

With the help of Lemmas 1 and 2, it is obvious that for strong Allee effect, $E_{1}$ is always saddle and the prey free equilibrium $E_{2}$ is always stable. For the equilibrium points $E_{3}$ and $E_{4}$ which represents the extinction of predator, $E_{3}$ is unstable if $β > 1$ , saddle if $β < 1$ and non-hyperbolic with one eigenvalue $1$ and the other not on the unit circle, if $β = 1$ , i-e, it may undergo transcritical or fold bifurcation. Whereas, $E_{4}$ is unstable if $β < 1$ , saddle if $β > 1$ and non-hyperbolic with one eigenvalue $1$ and the other not on the unit circle, if $β = 1$ , i-e, it may undergo transcritical or fold bifurcation.
For weak Allee effect, predator extinct equilibrium points $E_{1}$ , $E_{3}$ and $E_{4}$ are always unstable. The prey free equilibrium point $E_{2}$ is stable if $B > 0$ , saddle if $B < 0$ and there is a possibility of transcritical or fold bifurcation, since it is non-hyperbolic with one eigenvalue $1$ and the other not on the unit circle, if $B = 1$ .
For all of the above boundary equilibrium points, period-doubling or Neimark-Sacker bifurcations are not possible. The stability analysis of positive interior stationary points is richer and valid for our model.

Theorem 1. If $β = (1 + ξ \pm 2 \sqrt{(1 + θ) ξ})$ , then $E_{5^{*}}$ may undergo transcritical or fold bifurcation with eigenvalues, $z_{1} = 1$ and $| z_{2} | \neq 1$ if and only if $h \neq \frac{2}{ρ - A D}$ and $ρ \neq A D$ .

Proof. Suppose, $β = (1 + ξ \pm 2 \sqrt{(1 + θ) ξ})$ . Then, for strong Allee effect, $E_{5^{*}}$ will always exist, while for weak Allee effect, it will exist if additionally $(1 - ξ)^{2} < 4 ξ θ$ . The jacobian of the system (1), at $E_{5^{*}}$ is $J^{*} = (\begin{array}{cc} 1 + h A D & - h A D \\ h ρ & 1 - h ρ \end{array})$ and the eigenvalues are $z_{1} = 1$ and $z_{2} = 1 + h (a d - ρ)$ . Thus, $| z_{2} | \neq 1$ if and only if $h \neq \frac{2}{ρ - A D}$ and $ρ \neq A D$ .

Theorem 2. For the positive interior fixed points $E_{5}$ and $E_{6}$ , suppose $1 + β > ξ$ , $β + θ ξ > 0$ and $β < 1 + ξ - 2 \sqrt{ξ (1 + θ)}$ and let, $\begin{aligned} h_{1} & = \frac{(C - D) {\bar{x}}_{1} + ρ}{ρ C {\bar{x}}_{1}} - \frac{\sqrt{{((C - D) {\bar{x}}_{1} + ρ)}^{2} - 4 ρ C {\bar{x}}_{1}}}{ρ C {\bar{x}}_{1}}, \\ h_{2} & = \frac{(C - D) {\bar{x}}_{1} + ρ}{ρ C {\bar{x}}_{1}} + \frac{\sqrt{{((C - D) {\bar{x}}_{1} + ρ)}^{2} - 4 ρ C {\bar{x}}_{1}}}{ρ C {\bar{x}}_{1}} . \end{aligned}$ Then, for both strong and weak Allee effects, the following holds true.

(i) If $h \in (0, h_{1})$ , then the fixed point is a sink.

(ii) If $h \in (h_{2}, \infty)$ , then the fixed point is a source.

(iii) The fixed is saddle if $h \in (h_{1}, h_{2})$ .

(iv) The fixed point is non-hyperbolic with eigenvalues $z_{1} = - 1$ and $| z_{2} | \neq 1$ if $h = h_{1}$ or $h = h_{2}$ and $h \neq \frac{2}{ρ - {\bar{x}}_{1} (D - C)}$ and $h \neq \frac{4}{ρ - {\bar{x}}_{1} (D - C)}$ .

(v) If, $\begin{aligned} ρ & \in ((D - C) {\bar{x}}_{1}, \infty) \cap ({(\sqrt{C} - \sqrt{D})}^{2} {\bar{x}}_{1}, {(\sqrt{C} + \sqrt{D})}^{2} {\bar{x}}_{1}) \end{aligned}$ then, the fixed point is non-hyperbolic with complex conjugate eigenvalues, $| z_{1} | = 1 = | z_{2} |$ if and only if $h = \frac{ρ - {\bar{x}}_{1} (D - C)}{C ρ {\bar{x}}_{1}} < \frac{4}{ρ - {\bar{x}}_{1} (D - C)}$ .

Proof. The characteristic polynomial of the system at the stationary points is, $p (z) = z^{2} - t z + d e t,$ where, $t = 2 + h ({\bar{x}}_{1} (D - C) - ρ)$ and $d e t = 1 + h ({\bar{x}}_{1} (D - C) - ρ + h ρ {\bar{x}}_{1} C)$ . Then, we can use lemma 1, since for any $h > 0$ , $p (- 1) = ρ C {\bar{x}}_{1} h^{2} - 2 ((C - D) {\bar{x}}_{1} + ρ) h + 4,$ which shows that if we define $h_{1, 2}$ as given, then ${\begin{cases} p (- 1) > 0, & h \in (0, h_{1}) \cap (h_{2}, \infty); \\ p (- 1) < 0, & h \in (h_{1}, h_{2}); \\ p (- 1) = 0, & h = h_{1} or h = h_{2} . \end{cases}$ Also, ${\begin{cases} d e t > 1, & h \in (\frac{ρ + {\bar{x}}_{1} (C - D)}{ρ C {\bar{x}}_{1}}, \infty); \\ d e t < 1, & h \in (0, \frac{ρ + {\bar{x}}_{1} (C - D)}{ρ C {\bar{x}}_{1}}); \\ d e t = 1, & h = \frac{ρ + {\bar{x}}_{1} (C - D)}{ρ C {\bar{x}}_{1}} . \end{cases}$ Then, using Lemmas 1 and 2, we arrive at the desired results.

3. Bifurcations

Note that the mathematical method and results for proving the period-doubling and Neimark-Sacker bifurcations, for $E_{5}$ and $E_{6}$ , are exactly similar. So from now onward, we will only show the mathematical results for the equilibrium point $E_{5}$ . The mathematical result for $E_{6}$ are identical. Define, $\begin{aligned} Ω_{P D} & = {(β, θ, ξ, ρ, γ, h_{0}) \in R_{+}^{6} : h_{0} = h_{1} or h_{2}, 1 + β > ξ, h \neq \frac{2}{ρ - x_{1} (D - C)}, h \neq \frac{4}{ρ - x_{1} (D - C)}}, \\ Ω_{N S} & = {(β, θ, ξ, ρ, γ, h_{0}) \in R_{+}^{6} : 1 + β > ξ, h_{0} = \frac{ρ - {\bar{x}}_{1} (D - C)}{C ρ {\bar{x}}_{1}} < \frac{4}{ρ - {\bar{x}}_{1} (D - C)}} . \end{aligned}$

Let $H$ be the parameter for the mapping (2), such that $| H | ⋘ 1$ , to obtain bifurcations. We can write the perturbed mapping as $\begin{aligned} [\begin{array}{c} x \\ y \end{array}] \to [\begin{array}{c} x + (h + H) x (\frac{(1 - x) (x - β)}{x + θ} - ξ \frac{y}{x + γ}) \\ y + (h + H) ρ y (1 - \frac{y}{x + γ}) \end{array}] . \end{aligned}$

Define $X = x - {\bar{x}}_{1}$ and $Y = y - {\bar{y}}_{1}$ . Using the series expansion method, we can write $\begin{matrix} (4) & [\begin{matrix} X \\ Y \end{matrix}] \to [\begin{array}{cc} a_{11} & a_{12} \\ a_{21} & a_{22} \end{array}] [\begin{matrix} X \\ Y \end{matrix}] + [\begin{matrix} f (X, Y, H) \\ g (X, Y, H) \end{matrix}], \end{matrix}$ where $a_{11} = 1 + h b_{3}$ , $a_{12} = h b_{4}$ , $a_{21} = h ρ$ , $a_{22} = 1 - h ρ$ , $f (X, Y, H) = h b_{1} X^{2} + h b_{2} X Y + b_{3} X H + b_{4} Y H + b_{5} X^{3} + b_{6} X^{2} Y + b_{1} X^{2} H + b_{2} X Y H + O ({| X + Y + H |}^{4}),$ $g (X, Y, H) = - h b_{7} X^{2} + 2 h b_{7} X Y + ρ X H - h b_{7} Y^{2} + h a_{4} X^{3} - 2 h a_{4} X^{2} Y - a_{4} X^{2} H ρ Y H + 2 b_{7} X Y H + h a_{4} X Y^{2} - b_{7} Y^{2} H + O ({| X + Y + H |}^{4})$ and $b_{1} = θ ({\bar{x}}_{1} + θ) a_{1} + {\bar{y}}_{1} a_{3} - 1,$ $b_{2} = - {\bar{y}}_{1} a_{3},$ $b_{3} = a_{2} + {\bar{x}}_{1} (({\bar{x}}_{1} + θ)^{2} a_{1} - 1) - {\bar{y}}_{1}^{2} a_{3},$ $b_{4} = - {\bar{x}}_{1} {\bar{y}}_{1}^{2} a_{3},$ $b_{5} = - h (θ a_{1} + a_{3}),$ $b_{6} = h a_{3},$ $b_{7} = \frac{ρ}{{\bar{y}}_{1}},$ $a_{1} = \frac{(1 + θ) (β + θ)}{({\bar{x}}_{1} + θ)^{4}},$ $a_{2} = ξ \frac{{\bar{x}}_{1}}{{\bar{y}}_{1}},$ $a_{3} = \frac{γ ξ}{{\bar{y}}_{1}^{3}},$ $a_{4} = \frac{ρ}{{\bar{y}}_{1}^{2}} .$

3.1. Period-Doubling Bifurcation

Assume that the parameters $(β, θ, ξ, ρ, γ, h_{0}) \in Ω_{P D}$ . Then variation of $h$ in some small neighborhood of $h_{0}$ gives rise to period-doubling bifurcation. For conversion into the normal form of period-doubling bifurcation, let $[\begin{matrix} X \\ Y \end{matrix}] = T_{P D} [\begin{matrix} u \\ v \end{matrix}],$ where $T_{P D}$ is an invertible matrix defined as, $T_{P D} = [\begin{array}{cc} a_{12} & a_{12} \\ - 1 - a_{11} & z_{2} - a_{22} \end{array}] .$ Using $T_{P D}$ , we can transform our system to the one given below. $\begin{matrix} (5) & \begin{aligned} [\begin{array}{c} u \\ v \end{array}] \to [\begin{array}{cc} - 1 & 0 \\ 0 & z_{2} \end{array}] [\begin{array}{c} u \\ v \end{array}] + [\begin{array}{c} f_{P D} (H, u, v) \\ g_{P D} (H, u, v) \end{array}], \end{aligned} \end{matrix}$ where, $\begin{aligned} f_{P D} (H, u, v) = & d_{1} u^{2} + d_{2} u v + d_{3} H u + d_{4} H v + d_{5} v^{2} + d_{6} u^{3} + d_{7} u^{2} v + d_{8} H u^{2} \\ + d_{9} H u v + d_{10} u v^{2} + d_{11} H v^{2} + d_{12} v^{3} + O ({| H + u + v |}^{4}), \\ g_{P D} (H, u, v) = & d_{13} u^{2} + d_{14} u v + d_{15} H u + d_{16} H v + d_{17} v^{2} + d_{18} u^{3} + d_{19} u^{2} v + d_{20} H u^{2} \\ + d_{21} H u v + d_{22} u v^{2} + d_{23} H v^{2} + d_{24} v^{3} + O ({| H + u + v |}^{4}), \end{aligned}$ and $\begin{aligned} d_{1} = & - c_{1} h (a_{12}^{2} (b_{1} c_{2} - a_{4}) - (a_{11} + 1)^{2} b_{7} - a_{12} (a_{11} + 1) (b_{2} c_{2} + 2 b_{7})), \\ d_{2} = & - c_{1} h (2 a_{12}^{2} (a_{4} - b_{1} c_{2}) + a_{12} (a_{11} + a_{22} - z_{2} + 1) (b_{2} c_{2} + 2 b_{7}) + 2 (a_{11} + 1) b_{7} (a_{22} - z_{2})), \\ d_{3} = & c_{1} (a_{12} (b_{3} c_{2} + ρ) - (1 + a_{11}) (b_{4} c_{2} - ρ)), \\ d_{4} = & c_{1} (a_{12} (b_{3} c_{2} + ρ) + (z_{2} - a_{22}) (b_{4} c_{2} - ρ)), \\ d_{5} = & - c_{1} h (a_{12} (a_{22} - z_{2}) (b_{2} c_{2} + 2 b_{7}) + a_{12}^{2} (b_{7} - b_{1} c_{2}) + b_{7} (a_{22} - z_{2})^{2}), \\ d_{6} = & a_{12} c_{1} (- a_{12} (a_{11} + 1) (b_{6} c_{2} - 2 a_{4} h) + a_{12}^{2} (a_{4} h + b_{5} c_{2}) + (a_{11} + 1)^{2} a_{4} h), \\ d_{7} = & a_{12} c_{1} (+ 3 a_{12}^{2} (a_{4} h + b_{5} c_{2}) + a_{4} (a_{11} + 1) h (a_{11} + 2 a_{22} - 2 z_{2} + 1) \\ - a_{12} (2 a_{11} + a_{22} - z_{2} + 2) (b_{6} c_{2} - 2 a_{4} h)), \\ d_{8} = & c_{1} (a_{12}^{2} (b_{1} c_{2} - a_{4}) - (a_{11} + 1)^{2} b_{7} - a_{12} (a_{11} + 1) (b_{2} c_{2} + 2 b_{7})), \\ d_{9} = & - c_{1} (2 (a_{11} + 1) b_{7} (a_{22} - z_{2}) + 2 a_{12}^{2} (a_{4} - b_{1} c_{2}) + a_{12} (a_{11} + a_{22} - z_{2} + 1) (b_{2} c_{2} + 2 b_{7})), \\ d_{10} = & a_{12} c_{1} (3 a_{12}^{2} (a_{4} h + b_{5} c_{2}) - a_{12} (a_{11} + 2 a_{22} - 2 z_{2} + 1) (b_{6} c_{2} - 2 a_{4} h) \\ + a_{4} h (a_{22} - z_{2}) (2 a_{11} + a_{22} - z_{2} + 2)), \\ d_{11} = & - c_{1} (a_{12} (a_{22} - z_{2}) (b_{2} c_{2} + 2 b_{7}) + a_{12}^{2} (b_{7} - b_{1} c_{2}) + b_{7} (a_{22} - z_{2})^{2}), \\ d_{12} = & a_{12} c_{1} (- a_{12} (a_{22} - z_{2}) (b_{6} c_{2} - 2 a_{4} h) + a_{12}^{2} (a_{4} h + b_{5} c_{2}) + a_{4} h (a_{22} - z_{2})^{2}), \\ d_{13} = & - h c_{1} ((1 + a_{11} + a_{12})^{2} b_{7} + a_{12} (- a_{12} b_{1} + b_{2} + a_{11} b_{2}) c_{3}), \\ d_{14} = & - h c_{1} (2 a_{12}^{2} (b_{7} - b_{1} c_{3}) + 2 (1 + a_{11}) b_{7} (a_{22} - z_{2}) + a_{12} (2 b_{7} + b_{2} c_{3}) (1 + a_{11} + a_{22} - z_{2})), \\ d_{15} = & c_{1} (a_{12} (b_{3} c_{3} + ρ) + (1 + a_{11}) (- b_{4} c_{3} + ρ)), \\ d_{16} = & c_{1} ((- a_{22} + z_{2}) (b_{4} c_{3} - ρ) + a_{12} (b_{3} c_{3} + ρ)), \\ d_{17} = & - c_{1} h (a_{12}^{2} (b_{7} - b_{1} c_{3}) + b_{7} (a_{22} - z_{2})^{2} + a_{12} (2 b_{7} + b_{2} c_{3}) (a_{22} - z_{2})), \\ d_{18} = & a_{12} c_{1} (a_{12} (a_{12} b_{5} - (1 + a_{11}) b_{6}) c_{3} + (1 + a_{11} + a_{12})^{2} a_{4} h), \\ d_{19} = & a_{12} c_{1} ((1 + a_{11}) a_{4} h (1 + a_{11} + 2 a_{22} - 2 z_{2}) + 3 a_{12}^{2} (b_{5} c_{3} + a_{4} h) \\ - a_{12} (b_{6} c_{3} - 2 a_{4} h) (2 + 2 a_{11} + a_{22} - z_{2})), \\ d_{20} = & - c_{1} ((1 + a_{11})^{2} b_{7} + a_{12}^{2} (a_{4} - b_{1} c_{3}) + (1 + a_{11}) a_{12} (2 b_{7} + b_{2} c_{3})), \\ d_{21} = & - c_{1} (2 (1 + a_{11}) b_{7} (a_{22} - z_{2}) + 2 a_{12}^{2} (a_{4} - b_{1} c_{3}) + a_{12} (b_{6} c_{3} - 2 a_{4} h) (2 + 2 a_{11} + a_{22} - z_{2})), \\ d_{22} = & a_{12} c_{1} (3 a_{12}^{2} (b_{5} c_{3} + a_{4} h) - a_{12} (b_{6} c_{3} - 2 a_{4} h) (1 + a_{11} + 2 a_{22} - 2 z_{2}) \\ + a_{4} h (a_{22} - z_{2}) (2 + 2 a_{11} + a_{22} - z_{2})), \\ d_{23} = & - c_{1} (a_{12} (2 b_{7} + b_{2} c_{3}) (a_{22} - z_{2}) + b_{7} (a_{22} - z_{2})^{2} + a_{12}^{2} (a_{4} - b_{1} c_{3})), \\ d_{24} = & a_{12} c_{1} (- a_{12} (b_{6} c_{3} - 2 a_{4} h) (a_{22} - z_{2}) + a_{4} h (a_{22} - z_{2})^{2} + a_{12}^{2} (b_{5} c_{3} + a_{4} h)), \\ X = & a_{12} (u + v), \\ Y = & - (1 + a_{11}) u + (z_{2} - a_{22}) v . \end{aligned}$

In order to determine the center manifold of map (5) at fixed point $(0, 0)$ , define, $\begin{aligned} M_{P D} & = {(H, u, v) : v = G_{P D} (H, u), | u | < δ_{1}, | H | < δ_{2}, δ_{1}, δ_{2} \in R^{+}, G_{P D} (0, 0) = 0, D G_{P D} (0, 0) = 0}, \end{aligned}$ where $G_{P D} (H, u) = - \frac{d_{15}}{1 + z_{2}} H u + O ({| H + u |}^{3})$ . Then the following map illustrates the dynamics restricted to $M_{P D}$ . $\begin{matrix} (6) & \begin{aligned} F : u & \mapsto - u + d_{1} u^{2} + d_{3} H u - \frac{d_{4} d_{15}}{1 + z_{2}} H^{2} u + (d_{8} - \frac{d_{2} d_{15}}{1 + z_{2}}) H u^{2} + d_{6} u^{3} + O ({| H + u |}^{4}) . \end{aligned} \end{matrix}$

In order to undergo a period-doubling bifurcation for map (6), the following two quantities, $l_{1}$ and $l_{2}$ , must be non-zero, i-e, $\begin{aligned} l_{1} & = & {(\frac{\partial^{2} F}{\partial H \partial u} + \frac{1}{2} \frac{\partial F}{\partial H} \frac{\partial^{2} F}{\partial u^{2}}) |}_{(0, 0)} = d_{3} \neq 0, \\ l_{2} & = & {(\frac{1}{6} \frac{\partial^{3} F}{\partial u^{3}} + {(\frac{1}{2} \frac{\partial^{2} F}{\partial u^{2}})}^{2}) |}_{(0, 0)} = d_{1}^{2} + d_{6} \neq 0. \end{aligned}$

Thus, from the above analysis and the theorem given in [41], we obtain the following result.

Theorem 3. Suppose that $l_{1} \neq 0$ and $l_{2} \neq 0$ , then system exhibits period-doubling bifurcation at the interior equilibrium point, if $h$ varies around $h_{1}$ or $h_{2}$ . Moreover, if $l_{2} > 0$ , then stable otherwise unstable period-two orbits bifurcate from the fixed point.

3.2. Neimark-Sacker Bifurcation

Let $(β, θ, ξ, ρ, γ, h_{0}) \in Ω_{N S}$ . Then, variation of $h$ around $h_{0}$ , gives emergence to Neimark-Sacker bifurcation. Instead of (4), we can write the perturbed system as $\begin{matrix} (7) & [\begin{matrix} X \\ Y \end{matrix}] \to [\begin{array}{cc} c_{11} & c_{12} \\ c_{21} & c_{22} \end{array}] [\begin{matrix} X \\ Y \end{matrix}] + [\begin{matrix} \bar{f} (X, Y) \\ \bar{g} (X, Y) \end{matrix}], \end{matrix}$ where $c_{11} = 1 + (h + H) b_{3}$ , $c_{12} = (h + H) b_{4}$ , $c_{21} = (h + H) ρ$ , $c_{22} = 1 - (h + H) ρ$ , $\begin{aligned} \bar{f} (X, Y) & = (h + H) b_{1} X^{2} + (h + H) b_{2} X Y + b_{5} X^{3} + b_{6} X^{2} Y + O ({| X + Y + H |}^{4}), \\ \bar{g} (X, Y) & = - (h + H) b_{7} X^{2} + 2 (h + H) b_{7} X Y - (h + H) b_{7} Y^{2} + h a_{4} X^{3} - 2 h a_{4} X^{2} Y \\ + h a_{4} X Y^{2} + O ({| X + Y + H |}^{4}), \end{aligned}$ and the above mentioned coefficients are defined at the start of this section. Let $p_{N S} (z)$ be the characteristic polynomial of matrix in (9), given by $\begin{matrix} (8) & p_{N S} (z) = z^{2} - P (H) z + Q (H), \end{matrix}$ where $P (H) = c_{11} + c_{22} and Q (H) = c_{11} c_{22} - c_{21} c_{12} .$

Since, $(β, θ, ξ, ρ, γ, h) \in Ω_{N S}$ , the roots of (8) are complex conjugate $z_{1}$ , $z_{2}$ and $| z_{1} | = | z_{2} | = 1$ It follows that $z_{1, 2} = \frac{P (H)}{2} \pm \frac{ι}{2} \sqrt{4 Q (H) - P (H)^{2}}$ and $1 = | z_{1} | = | z_{2} | = \sqrt{Q (H)}$ . Also, ${(\frac{d | z_{1, 2} |}{d H})}_{H = 0} \neq 0$ if and only if $\frac{(θ + 1) (β + θ)}{({\bar{x}}_{1} + θ)^{2}} \neq 1$ .

We also need to check that $z_{1, 2}^{m} \neq 1$ , for $m = 1, 2, 3, 4$ at $H = 0$ , since we do not want the characteristic polynomial to lie in the intersection of unit circle of coordinate axis. This is the same as checking when $P (0) \neq - 2, 0, 1, 2$ . Since, $(β, θ, ξ, ρ, γ, h) \in Ω_{N S}$ , $P (H)^{2} - 4 Q (H) < 0$ , which implies that $P (H)^{2} < 4 Q (H)$ or $P (H) \in (- 2, 2)$ . Thus, $P (0) \neq - 2, 0, 2$ . Finally, $\begin{aligned} P (0) & = h [\bar{x} (\frac{(θ + 1) (β + θ)}{{(\bar{x} + θ)}^{2}} + \frac{ξ}{\bar{y}} - 1) - ρ] + 2, \end{aligned}$ which implies that $P (0) \neq 1$ if and only if $h \neq \frac{1}{ρ - \bar{x} (\frac{(θ + 1) (β + θ)}{{(\bar{x} + θ)}^{2}} + \frac{ξ}{\bar{y}} - 1)}$ .

Now, to convert the system (4) to the normal form of Neimark-Sacker bifurcation, let $R = \frac{P (0)}{2}$ and $S = \frac{1}{2} \sqrt{4 Q (0) - P (0)^{2}}$ and let $T_{N S} = (\begin{array}{cc} 0 & 1 \\ S & R \end{array}) .$ Here, $T_{N S}$ is an invertible matrix. Consider the transformation $(\begin{matrix} X \\ Y \end{matrix}) = T_{N S} (\begin{matrix} u \\ v \end{matrix}) .$ Using the above given transformation, we can write system (4) as, $\begin{matrix} (9) & \begin{aligned} [\begin{array}{c} u \\ v \end{array}] & \to [\begin{array}{cc} R & - S \\ S & R \end{array}] [\begin{array}{c} u \\ v \end{array}] + [\begin{array}{c} f_{N S} (u, v) \\ g_{N S} (u, v) \end{array}] . \end{aligned} \end{matrix}$

$\begin{aligned} f_{N S} (u, v) = & e_{1} u^{2} + e_{2} u v + e_{3} v^{2} + e_{4} u^{2} v + e_{5} u v^{2} + e_{6} v^{3} + O ({| u + v |}^{4}), \\ g_{N S} (u, v) = & e_{7} u v + e_{8} v^{2} + e_{9} u v^{2} + e_{10} v^{3} + O ({| u + v |}^{4}), \end{aligned}$ where, $\begin{aligned} e_{1} & = - (h + H) b_{7} S, \\ e_{2} & = - (h + H) (b_{2} R + 2 b_{7} (R - 1)), \\ e_{3} & = - \frac{(h + H) (R (b_{1} + b_{2} R) + b_{7} (R - 1)^{2})}{S}, \\ e_{4} & = h a_{4} S, \\ e_{5} & = 2 a_{4} h (R - 1) - b_{6} R, \\ e_{6} & = \frac{a_{4} h (R - 1)^{2} - R (b_{5} + b_{6} R)}{S}, \\ e_{7} & = (h + H) b_{2} S, \\ e_{8} & = (h + H) (b_{1} + b_{2} R), \\ e_{9} & = b_{6} S, \\ e_{10} & = b_{5} + b_{6} R, \\ X & = v, Y = S u + R v . \end{aligned}$ We need the following quantity to be non-zero for the map (9) to undergo Neimark-Sacker bifurcation: $\begin{matrix} (10) & \begin{aligned} ζ_{N S} & = {[- Re (ϖ_{20} ϖ_{11} \frac{(1 - 2 z_{1}) z_{2}^{2}}{1 - z_{1}}) - \frac{1}{2} {| ϖ_{11} |}^{2} - {| ϖ_{02} |}^{2} + Re (ϖ_{21} z_{2})]}_{H = 0}, \end{aligned} \end{matrix}$ where, $\begin{aligned} ϖ_{20} & = \frac{1}{8} [\frac{\partial^{2} f_{N S}}{\partial u^{2}} - \frac{\partial^{2} f_{N S}}{\partial v^{2}} + 2 \frac{\partial^{2} g_{N S}}{\partial u \partial v} + ι (\frac{\partial^{2} g_{N S}}{\partial u^{2}} - \frac{\partial^{2} g_{N S}}{\partial v^{2}} - 2 \frac{\partial^{2} f_{N S}}{\partial u \partial v})], \\ ϖ_{11} & = \frac{1}{4} [\frac{\partial^{2} f_{N S}}{\partial u^{2}} + \frac{\partial^{2} f_{N S}}{\partial v^{2}} + ι (\frac{\partial^{2} g_{N S}}{\partial u^{2}} + \frac{\partial^{2} g_{N S}}{\partial v^{2}})], \\ ϖ_{02} & = \frac{1}{8} [\frac{\partial^{2} f_{N S}}{\partial u^{2}} - \frac{\partial^{2} f_{N S}}{\partial v^{2}} - 2 \frac{\partial^{2} g_{N S}}{\partial u \partial v} + ι (\frac{\partial^{2} g_{N S}}{\partial u^{2}} - \frac{\partial^{2} g_{N S}}{\partial v^{2}} + 2 \frac{\partial^{2} f_{N S}}{\partial u \partial v})], \\ ϖ_{21} & = \frac{1}{16} [\frac{\partial^{3} f_{N S}}{\partial u^{3}} + \frac{\partial^{3} f_{N S}}{\partial u \partial v^{2}} + \frac{\partial^{3} g_{N S}}{\partial u^{2} \partial v} + \frac{\partial^{3} g_{N S}}{\partial v^{3}} + ι (\frac{\partial^{3} g_{N S}}{\partial u^{3}} + \frac{\partial^{3} g_{N S}}{\partial u \partial v^{2}} - \frac{\partial^{3} f_{N S}}{\partial u^{2} \partial v} - \frac{\partial^{3} f_{N S}}{\partial v^{3}})] . \end{aligned}$ Based on the prior analysis, we have the following result [41].

Theorem 4. If $\frac{(θ + 1) (β + θ)}{({\bar{x}}_{1} + θ)^{2}} \neq 1$ and $ζ_{N S}$ defined in (10) is non-zero and parameter $h$ changes in the small neighborhood of $h_{0}$ then the model (3) exhibits Neimark-Sacker bifurcation at the fixed point $(\bar{x_{1}}, \bar{y_{1}})$ . Moreover, an attracting (resp., repelling) invariant closed curve bifurcate from the positive fixed point if $ζ_{N S} < 0$ (resp., $ζ_{N S} > 0$ ).

4. Numerical Simulations

In this section, we will use specific numerical values of parameters $(β, θ, ξ, ρ, γ)$ to present some examples which will show the emergence of period-doubling and Neimark-Sacker bifurcations for the system (2). We will use $h$ as the bifurcation parameter. The illustration will be done using phase portraits and bifurcation diagrams. We will also ratify Theorems 3 and 4, by showing that the numerical examples are according to these theorems.

4.1. Strong Allee Effect

4.1.1. Period-Doubling Bifurcation

Example 1. Select $(β, θ, ξ, ρ, γ)$ $=$ $(0.1, 0.01, 0.05, 0.05, 0.3)$ , $h \in [2.4, 3.345]$ and the initial value $(x_{0}, y_{0}) = (1.2, 2)$ . The jacobian matrix is, $(\begin{array}{cc} 1 - 0.790244 h & - 0.0379371 h \\ 0.05 h & 1 - 0.05 h \end{array}),$ and $z_{1} = 1 - 0.787672 h$ and $z_{2} = 1 - 0.0525714 h$ are respective eigenvalues which are both less then $1$ , for any $h > 0$ . For $h = h_{0} = 2.53913$ , we have $z_{1} = - 1$ and $z_{2} \in (- 1, 1)$ . Thus the fixed point is always stable with given parametric conditions and for any $h \neq h_{0}$ , or, $h = 38.0435$ . If $h = h_{0}$ , the fixed point exhibits period-doubling bifurcation when $h$ moves around $h_{0}$ . We can also see in the phase portrait given in Figure 2 $(b)$ that the fixed point $(0.943479, 1.24348)$ is stable for $h < h_{0}$ and the period-doubling bifurcation rises at $h_{0}$ . At around $h = 3.1$ another stable period-two orbit bifurcates which can be seen in Figure 2 $(c)$ . With these parametric conditions, we calculate and get $l_{1} = - 0.0421168 \neq 0$ and $l_{2} = 0.000652147 \neq 0$ . These values also verify the results obtained in Theorem 3. Finally since $l_{2} > 0$ , the period-two orbits that bifurcate from the fixed point is stable.

Figure 2: Phase Portraits with Initial Conditions $(X_{0}, Y_{0}) = (1.2, 2)$ , at Parametric Values $(β, Θ, Ξ, ρ, Γ)$ $=$ $(0.1, 0.01, 0.05, 0.05, 0.3)$ , Around the Bifurcation Parameter $H_{0} = 2.53913$

Figure 3: Emergence of Period-Doubling Bifurcation for Positive Fixed Point for the Model (), in the Interval $H \in [2.4, 3.345]$ , with Initial Conditions $(X_{0}, Y_{0}) = (1.2, 2)$ and the Parametric Values $(β, θ, ξ, ρ, γ)$ $=$ $(0.1, 0.01, 0.05, 0.05, 0.3)$

4.1.2. Neimark-Sacker Bifurcation

Example 2. Let us assume that for $h \in [2.1, 2.34]$ and the initial value $(x_{0}, y_{0}) = (0.98, 0.6)$ , we have the following parameter values, $(β, θ, ξ, ρ, γ)$ $=$ $(0.15, 0.06, 0.05, 1.1, 0.1)$ . The jacobian matrix is, $(\begin{array}{cc} 1 - 0.681574 h & - 0.0451768 h \\ 1.1 h & 1 - 1.1 h \end{array}),$ and $z_{1} = - 0.985179 - 0.171532 ι$ and $z_{2} = - 0.985179 + 0.171532 ι$ are the corresponding eigenvalues, for $h = h_{0} = 2.2285674489338163$ . Both eigenvalues have modulus equal to $1$ . Also, $h_{0} \in (0, \frac{4}{ρ - x_{1} (D - C)}) = (0, 2.24521)$ , ${(\frac{d | z_{1} |}{d H})}_{H = 0} = {(\frac{d | z_{2} |}{d H})}_{H = 0} = 1.78157 \neq 0$ and $ζ_{N S} = - 7834.72 \neq 0$ . Thus, the fixed point exhibits Neimark-Sacker bifurcation, at the equilibrium point $(0.936652, 1.03665)$ , as can be seen in Figure 5. Moreover, for $h > h_{0}$ we have $ζ_{N S} < 0$ , which shows that the invariant closed curve that bifurcates from the equilibrium point is attracting.

We can see a stable spiral in the phase portrait in Figure 4 $(a)$ which increases in size as the $h$ increases (Figure 4 $(b)$ ). With further increase in $h$ , the spiral changes into a closed curve which is invariant and attracting. Due to Neimark-Sacker bifurcation, the attractor has rough edges. These rough edges and spiral are plotted in Figure 4 $(c)$ at $h = 2.245$ . At $h = 2.3$ , it can be observed that the edges begin to disappear and the invariant closed curves continue around the positive fixed point, as can be seen in 4 $(d)$ and 4 $(e)$ .

Figure 4: Phase Space Portraits with Initial Conditions $(X_{0}, Y_{0}) = (0.98, 0.6)$ , at Parametric Values $(β, Θ, Ξ, ρ, Γ)$ $=$ $(0.15, 0.06, 0.05, 1.1, 0.1)$ , around $h_{0}$

Figure 5: Emergence of Neimark-Sacker Bifurcation for Positive Fixed Points in the Interval $H \in [2.1, 2.34]$ , of Model () with Initial Conditions $(X_{0}, Y_{0}) = (0.98, 0.6)$ and the Parametric Values $(β, θ, ξ, ρ, γ)$ $=$ $(0.15, 0.06, 0.05, 1.1, 0.1)$

4.2. Weak Allee Effect

4.2.1. Period-Doubling Bifurcation

Example 3. Select $(β, θ, ξ, ρ, γ)$ $=$ $(- 0.02, - \frac{β}{ξ}, 0.09, 0.01, 1.5)$ , $h \in [2.9, 3.92]$ and the initial value $(x_{0}, y_{0}) = (0.6, 2)$ . The jacobian matrix is, $(\begin{array}{cc} 1 - 0.678663 h & - 0.0335146 h \\ 0.01 h & 1 - 0.01 h \end{array}),$ and $z_{1} = 1 - 0.678162 h$ and $z_{2} = 1 - 0.0105016 h$ are respective eigenvalues which are both less then $1$ , for any $h > 0$ . For $h = h_{0} = 2.94915$ , we have $z_{1} = - 1$ and $z_{2} \in (- 1, 1)$ . Thus the fixed point is always stable with given parametric conditions and for any $h \neq h_{0}$ , or, $h = 190.447$ . If $h = h_{0}$ , the fixed point exhibits period-doubling bifurcation when $h$ moves around $h_{0}$ . We can also see in the phase portrait given in Figure 6 $(b)$ that the fixed point $(0.943479, 1.24348)$ is stable for $h < h_{0}$ and the period-doubling bifurcation rises at $h_{0}$ . At around $h = 3.1$ another stable period-two orbit bifurcates which can be seen in Figure 6 $(c)$ . With these parametric conditions, we calculate and get $l_{1} = - 0.0151578 \neq 0$ and $l_{2} = 0.00090295 \neq 0$ . These values also verify the results obtained in Theorem 3. Finally since $l_{2} > 0$ , the period-two orbits that bifurcate from the fixed point is stable.

Figure 6: Phase Portraits with Initial Conditions $(X_{0}, Y_{0}) = (0.6, 2)$ , at Parametric Values $(β, Θ, Ξ, ρ, Γ)$ $=$ $(- 0.02, - \frac{β}{ξ}, 0.09, 0.01, 1.5)$ , Around the Bifurcation Parameter $H_{0} = 2.94915$

Figure 7: Emergence of Period-Doubling Bifurcation for Positive Fixed Point, for the Model (), in the Interval $H \in [2.9, 3.92]$ , with Initial Conditions $(X_{0}, Y_{0}) = (0.6, 2)$ and the Parametric Values $(β, θ, ξ, ρ, γ)$ $=$ $(- 0.02, - \frac{β}{ξ}, 0.09, 0.01, 1.5)$

4.2.2. Neimark-Sacker Bifurcation

Example 4. Let us assume that for $h \in [3.23, 3.2]$ and the initial value $(x_{0}, y_{0}) = (0.6, 0.4)$ , we have the following parameter values, $(β, θ, ξ, ρ, γ)$ $=$ $(- 0.08, 0.06, 0.1, 0.9, 0.01)$ . The jacobian matrix is, $(\begin{array}{cc} 1 - 0.316267 h & - 0.0987952 h \\ 0.9 h & 1 - 0.9 h \end{array}),$ and $z_{1} = - 0.355266 - 0.136097 ι$ and $z_{2} = - 0.355266 + 0.136097 ι$ are the corresponding eigenvalues, for $h = h_{0} = 3.25592$ . Both eigenvalues have modulus equal to $1$ . Also, $h_{0} \in (0, \frac{4}{ρ - x_{1} (D - C)}) = (0, 3.28875)$ , ${(\frac{d | z_{1} |}{d H})}_{H = 0} = {(\frac{d | z_{2} |}{d H})}_{H = 0} = 0.832494 \neq 0$ and $ζ_{N S} = - 2237.02 \neq 0$ . Thus, the fixed point exhibits Neimark-Sacker bifurcation, at the equilibrium point $(0.82, 0.83)$ , as can be seen in Figure 9. Moreover, for $h > h_{0}$ we have $ζ_{N S} < 0$ , which shows that the invariant closed curve that bifurcates from the fixed point is attracting.

We can see a stable spiral in the phase portrait in Figure 8 $(a)$ which increases in size as the $h$ increases (Figure 8 $(b)$ ). With further increase in $h$ , the spiral changes into a closed curve which is invariant and attracting. Due to Neimark-Sacker bifurcation, the attractor has rough edges. These rough edges and spiral are plotted in Figure 8 $(c)$ at $h = 3.27$ . At $h = 3.288$ , it can be observed that the edges begin to disappear and the invariant closed curves continue around the positive fixed point, as can be seen in Figure 8 $(d)$ .

Figure 8: Phase Space Portraits with Initial Conditions $(X_{0}, Y_{0}) = (0.6, 0.4)$ , at Parametric Values $(β, Θ, Ξ, ρ, Γ)$ $=$ $(- 0.08, 0.06, 0.1, 0.9, 0.01)$ , Around $H_{0}$

Figure 9: Emergence of Neimark-Sacker Bifurcation for Positive Fixed Point in the Interval $H \in [3.23, 3.28]$ , of Model () with Initial Conditions $(X_{0}, Y_{0}) = (0.6, 0.4)$ and Parametric Values $(β, θ, ξ, ρ, γ)$ $=$ $(- 0.08, 0.06, 0.1, 0.9, 0.01)$

5. Conclusion

In this paper, we have extended the analysis of the modified Leslie-Gower model with double Allee effect by studying its discrete version and exploring its rich dynamics. We investigated the existence and stability of interior equilibrium points for both strong and weak Allee effects, and showed that the system undergoes period-doubling and Neimark-Sacker bifurcations under certain parametric conditions. To this end, we employed the center manifold theorem and our main results are presented in Theorems 1, 2, 3, and 4, which were complemented by numerical simulations in Section 4.

Our findings revealed that the step size $h$ can serve as a bifurcation parameter, allowing for the emergence of a plethora of dynamical behaviors. Specifically, our numerical simulations in Figures 5 and 9 demonstrated that the system displays orbits of period- $1$ , $2$ , $4$ , and $8$ , while period-doubling and Neimark-Sacker bifurcations were shown to occur for both strong and weak Allee effects, as depicted in Figures 3, 7, 5, and 9. Furthermore, we confirmed the validity of Theorems 3 and 4 through our numerical simulations.

In summary, our study provides a comprehensive understanding of the dynamical behavior of the modified Leslie-Gower model with double Allee effect in its discrete form. Our findings contribute to the growing body of literature on Allee effect models and can inform future research in the field.

Conflict of Interest

The authors declare no conflict of interest.

References:

Khan, M. S., Abbas, M., Bonyah, E. and Qi, H., 2022. Michaelis-Menten-Type Prey Harvesting in Discrete Modified Leslie-Gower Predator-Prey Model. Journal of Function Spaces, 2022, p.9575638.
Allesina, S. and Tang, S., 2012. Stability criteria for complex ecosystems. Nature, 483(7388), pp.205-208.
Clark, C.W., 1976. Mathematical bioeconomics: The Optimal Management Resources. John Wiley & Sons.
Lotka, A.J., 1920. Analytical note on certain rhythmic relations in organic systems. Proceedings of the National Academy of Sciences, 6(7), pp.410-415.
Volterra, V., 1928. Variations and fluctuations of the number of individuals in animal species living together. ICES Journal of Marine Science, 3(1), pp.3-51.
Embree, D.G., 1965. The population dynamics of the winter moth in Nova Scotia, 1954–1962. The Memoirs of the Entomological Society of Canada, 97(S46), pp.5-57.
LESLIE, P.H., 1958. A stochastic model for studying the properties of certain biological systems by numerical methods. Biometrika, 45(1-2), pp.16-31.
Hsu, S.B. and Huang, T.W., 1995. Global stability for a class of predator-prey systems. SIAM Journal on Applied Mathematics, 55(3), pp.763-783.
Chen, L. and Chen, F., 2009. Global stability of a Leslie–Gower predator–prey model with feedback controls. Applied Mathematics Letters, 22(9), pp.1330-1334.
May, R. M., 1973. Stability and Complexity in Model Ecosystems. Princeton University Press.
Hsu, S.B. and Hwang, T.W., 1998. Uniqueness of limit cycles for a predator-prey system of Holling and Leslie type. Canad. Appl. Math. Quart, 6(2), pp.91-117.
Hsu, S.B. and Hwang, T.W., 1999. Hopf bifurcation analysis for a predator-prey system of Holling and Leslie type. Taiwanese Journal of Mathematics, 3(1), pp.35-53.
Huang, J., Ruan, S. and Song, J., 2014. Bifurcations in a predator–prey system of Leslie type with generalized Holling type III functional response. Journal of Differential Equations, 257(6), pp.1721-1752.
Aziz-Alaoui, M.A. and Okiye, M.D., 2003. Boundedness and global stability for a predator-prey model with modified Leslie-Gower and Holling-type II schemes. Applied Mathematics Letters, 16(7), pp.1069-1075.
Caughley, G., 1976. Plant-herbivore systems. In R. M. May (Ed.), Theoretical ecology: Principles and Applications, (pp. 94-113). Philadelphia, PA: W. B. Saunders Co.
Wollkind, D.J. and Logan, J.A., 1978. Temperature-dependent predator-prey mite ecosystem on apple tree foliage. Journal of Mathematical Biology, 6, pp.265-283.
Wollkind, D.J., Collings, J.B. and Logan, J.A., 1988. Metastability in a temperature-dependent model system for predator-prey mite outbreak interactions on fruit trees. Bulletin of Mathematical Biology, 50(4), pp.379-409.
Liu, Y. and Zeng, Z., 2019. Analysis of a predator-prey model with Crowley-Martin and modified Leslie-Gower schemes with stochastic perturbation. Journal of Applied Analysis & Computation, 9(6), pp.2409-2435.
Jiang, J. and Song, Y., 2013, January. Stability and bifurcation analysis of a delayed Leslie-Gower predator-prey system with nonmonotonic functional response. In Abstract and Applied Analysis (Vol. 2013), p.152459, 19 pages.
Drake, J.M. and Kramer, A.M., 2011. Allee effects. Nature Education Knowledge, 3(10), p.2.
Stephens, P.A. and Sutherland, W.J., 1999. Consequences of the Allee effect for behaviour, ecology and conservation. Trends in Ecology & Evolution, 14(10), pp.401-405.
Lewis, M.A. and Kareiva, P., 1993. Allee dynamics and the spread of invading organisms. Theoretical Population Biology, 43(2), pp.141-158.
Courchamp, F., Berec, L. and Gascoigne, J., 2008. Allee Effects in Ecology and Conservation. Oxford University Press.
González-Olivares, E.D.U.A.R.D.O., Gonzalez-Yanez, B., Mena-Lorca, J. and Ramos-Jiliberto, R., 2006. Modelling the Allee effect: are the different mathematical forms proposed equivalents. In Proceedings of the International Symposium on Mathematical and Computational Biology BIOMAT (Vol. 2007, pp. 53-71).
Kent, A., Doncaster, C.P. and Sluckin, T., 2003. Consequences for predators of rescue and Allee effects on prey. Ecological Modelling, 162(3), pp.233-245.
Zhou, S.R., Liu, Y.F. and Wang, G., 2005. The stability of predator–prey systems subject to the Allee effects. Theoretical Population Biology, 67(1), pp.23-31.
Berec, L., Angulo, E. and Courchamp, F., 2007. Multiple Allee effects and population management. Trends in Ecology & Evolution, 22(4), pp.185-191.
Gascoigne, J. and Lipcius, R.N., 2004. Allee effects in marine systems. Marine Ecology Progress Series, 269, pp.49-59.
González-Olivares, E., González-Yañez, B., Lorca, J.M., Rojas-Palma, A. and Flores, J.D., 2011. Consequences of double Allee effect on the number of limit cycles in a predator–prey model. Computers & Mathematics with Applications, 62(9), pp.3449-3463.
Huincahue-Arcos, J. and González-Olivares, E., 2013. The Rosenzweig-MacArthur predation model with double Allee effects on prey. In Proceedings of the International Conference on Applied Mathematics and Computational Methods in Engineering (pp. 206-211).
Flores, J.D. and González-Olivares, E., 2014. Dynamics of a predator–prey model with Allee effect on prey and ratio-dependent functional response. Ecological complexity, 18, pp.59-66.
Feng, P. and Kang, Y., 2015. Dynamics of a modified Leslie–Gower model with double Allee effects. Nonlinear Dynamics, 80, pp.1051-1062.
Pal, P.J. and Saha, T., 2015. Qualitative analysis of a predator–prey system with double Allee effect in prey. Chaos, Solitons & Fractals, 73, pp.36-63.
Kuznetsov, Y.A., Kuznetsov, I.A. and Kuznetsov, Y., 1998. Elements of Applied Bifurcation Theory (Vol. 112, pp. xx+-591). New York: Springer.
Gukenheimer, J. and Holmes, P., 1983. Nonlinear oscillations. Dynamical Systems, and Bifurcation of Vector Fields, Springer-Verlag, NY.
Turchin, P.: Complex Population Dynamics, p. 35. Princeton University Press, Princeton.
Brauer, F., Castillo-Chavez, C. and Castillo-Chavez, C., 2012. Mathematical Models in Population Biology and Epidemiology (Vol. 2, No. 40). New York: springer.
Rof, A. and Krishnamoorthy, A., 2022. On a queueing inventory with common life time and reduction sale consequent to increase in age. 3c Empresa: Investigación y Pensamiento Crítico, 11(2), pp.15-31.
Hu, S., Meng, Q., Xu, D. and Al-Juboori, U.A., 2021. The optimal solution of feature decomposition based on the mathematical model of nonlinear landscape garden features. Applied Mathematics and Nonlinear Sciences, 7(1), pp.751-760.
Singh, M.K., Bhadauria, B.S. and Singh, B.K., 2018. Bifurcation analysis of modified Leslie-Gower predator-prey model with double Allee effect. Ain Shams Engineering Journal, 9(4), pp.1263-1277.
Wiggins, S., 1990. Applied Nonlinear Dynamics and Chaos. New York: Spring-Verlag.

[1] Khan, M. S., Abbas, M., Bonyah, E. and Qi, H., 2022. Michaelis-Menten-Type Prey Harvesting in Discrete Modified Leslie-Gower Predator-Prey Model. Journal of Function Spaces, 2022, p.9575638.

[2] Allesina, S. and Tang, S., 2012. Stability criteria for complex ecosystems. Nature, 483(7388), pp.205-208.

[3] Clark, C.W., 1976. Mathematical bioeconomics: The Optimal Management Resources. John Wiley & Sons.

[4] Lotka, A.J., 1920. Analytical note on certain rhythmic relations in organic systems. Proceedings of the National Academy of Sciences, 6(7), pp.410-415.

[5] Volterra, V., 1928. Variations and fluctuations of the number of individuals in animal species living together. ICES Journal of Marine Science, 3(1), pp.3-51.

[6] Embree, D.G., 1965. The population dynamics of the winter moth in Nova Scotia, 1954–1962. The Memoirs of the Entomological Society of Canada, 97(S46), pp.5-57.

[7] LESLIE, P.H., 1958. A stochastic model for studying the properties of certain biological systems by numerical methods. Biometrika, 45(1-2), pp.16-31.

[8] Hsu, S.B. and Huang, T.W., 1995. Global stability for a class of predator-prey systems. SIAM Journal on Applied Mathematics, 55(3), pp.763-783.

[9] Chen, L. and Chen, F., 2009. Global stability of a Leslie–Gower predator–prey model with feedback controls. Applied Mathematics Letters, 22(9), pp.1330-1334.

[10] May, R. M., 1973. Stability and Complexity in Model Ecosystems. Princeton University Press.

[11] Hsu, S.B. and Hwang, T.W., 1998. Uniqueness of limit cycles for a predator-prey system of Holling and Leslie type. Canad. Appl. Math. Quart, 6(2), pp.91-117.

[12] Hsu, S.B. and Hwang, T.W., 1999. Hopf bifurcation analysis for a predator-prey system of Holling and Leslie type. Taiwanese Journal of Mathematics, 3(1), pp.35-53.

[13] Huang, J., Ruan, S. and Song, J., 2014. Bifurcations in a predator–prey system of Leslie type with generalized Holling type III functional response. Journal of Differential Equations, 257(6), pp.1721-1752.

[14] Aziz-Alaoui, M.A. and Okiye, M.D., 2003. Boundedness and global stability for a predator-prey model with modified Leslie-Gower and Holling-type II schemes. Applied Mathematics Letters, 16(7), pp.1069-1075.

[15] Caughley, G., 1976. Plant-herbivore systems. In R. M. May (Ed.), Theoretical ecology: Principles and Applications, (pp. 94-113). Philadelphia, PA: W. B. Saunders Co.

[16] Wollkind, D.J. and Logan, J.A., 1978. Temperature-dependent predator-prey mite ecosystem on apple tree foliage. Journal of Mathematical Biology, 6, pp.265-283.

[17] Wollkind, D.J., Collings, J.B. and Logan, J.A., 1988. Metastability in a temperature-dependent model system for predator-prey mite outbreak interactions on fruit trees. Bulletin of Mathematical Biology, 50(4), pp.379-409.

[18] Liu, Y. and Zeng, Z., 2019. Analysis of a predator-prey model with Crowley-Martin and modified Leslie-Gower schemes with stochastic perturbation. Journal of Applied Analysis & Computation, 9(6), pp.2409-2435.

[19] Jiang, J. and Song, Y., 2013, January. Stability and bifurcation analysis of a delayed Leslie-Gower predator-prey system with nonmonotonic functional response. In Abstract and Applied Analysis (Vol. 2013), p.152459, 19 pages.

[20] Drake, J.M. and Kramer, A.M., 2011. Allee effects. Nature Education Knowledge, 3(10), p.2.

[21] Stephens, P.A. and Sutherland, W.J., 1999. Consequences of the Allee effect for behaviour, ecology and conservation. Trends in Ecology & Evolution, 14(10), pp.401-405.

[22] Lewis, M.A. and Kareiva, P., 1993. Allee dynamics and the spread of invading organisms. Theoretical Population Biology, 43(2), pp.141-158.

[23] Courchamp, F., Berec, L. and Gascoigne, J., 2008. Allee Effects in Ecology and Conservation. Oxford University Press.

[24] González-Olivares, E.D.U.A.R.D.O., Gonzalez-Yanez, B., Mena-Lorca, J. and Ramos-Jiliberto, R., 2006. Modelling the Allee effect: are the different mathematical forms proposed equivalents. In Proceedings of the International Symposium on Mathematical and Computational Biology BIOMAT (Vol. 2007, pp. 53-71).

[25] Kent, A., Doncaster, C.P. and Sluckin, T., 2003. Consequences for predators of rescue and Allee effects on prey. Ecological Modelling, 162(3), pp.233-245.

[26] Zhou, S.R., Liu, Y.F. and Wang, G., 2005. The stability of predator–prey systems subject to the Allee effects. Theoretical Population Biology, 67(1), pp.23-31.

[27] Berec, L., Angulo, E. and Courchamp, F., 2007. Multiple Allee effects and population management. Trends in Ecology & Evolution, 22(4), pp.185-191.

[28] Gascoigne, J. and Lipcius, R.N., 2004. Allee effects in marine systems. Marine Ecology Progress Series, 269, pp.49-59.

[29] González-Olivares, E., González-Yañez, B., Lorca, J.M., Rojas-Palma, A. and Flores, J.D., 2011. Consequences of double Allee effect on the number of limit cycles in a predator–prey model. Computers & Mathematics with Applications, 62(9), pp.3449-3463.

[30] Huincahue-Arcos, J. and González-Olivares, E., 2013. The Rosenzweig-MacArthur predation model with double Allee effects on prey. In Proceedings of the International Conference on Applied Mathematics and Computational Methods in Engineering (pp. 206-211).

[31] Flores, J.D. and González-Olivares, E., 2014. Dynamics of a predator–prey model with Allee effect on prey and ratio-dependent functional response. Ecological complexity, 18, pp.59-66.

[32] Feng, P. and Kang, Y., 2015. Dynamics of a modified Leslie–Gower model with double Allee effects. Nonlinear Dynamics, 80, pp.1051-1062.

[33] Pal, P.J. and Saha, T., 2015. Qualitative analysis of a predator–prey system with double Allee effect in prey. Chaos, Solitons & Fractals, 73, pp.36-63.

[34] Kuznetsov, Y.A., Kuznetsov, I.A. and Kuznetsov, Y., 1998. Elements of Applied Bifurcation Theory (Vol. 112, pp. xx+-591). New York: Springer.

[35] Gukenheimer, J. and Holmes, P., 1983. Nonlinear oscillations. Dynamical Systems, and Bifurcation of Vector Fields, Springer-Verlag, NY.

[36] Turchin, P.: Complex Population Dynamics, p. 35. Princeton University Press, Princeton.

[37] Brauer, F., Castillo-Chavez, C. and Castillo-Chavez, C., 2012. Mathematical Models in Population Biology and Epidemiology (Vol. 2, No. 40). New York: springer.

[38] Rof, A. and Krishnamoorthy, A., 2022. On a queueing inventory with common life time and reduction sale consequent to increase in age. 3c Empresa: Investigación y Pensamiento Crítico, 11(2), pp.15-31.

[39] Hu, S., Meng, Q., Xu, D. and Al-Juboori, U.A., 2021. The optimal solution of feature decomposition based on the mathematical model of nonlinear landscape garden features. Applied Mathematics and Nonlinear Sciences, 7(1), pp.751-760.

[40] Singh, M.K., Bhadauria, B.S. and Singh, B.K., 2018. Bifurcation analysis of modified Leslie-Gower predator-prey model with double Allee effect. Ain Shams Engineering Journal, 9(4), pp.1263-1277.

[41] Wiggins, S., 1990. Applied Nonlinear Dynamics and Chaos. New York: Spring-Verlag.

Contents

Utilitas Mathematica

A Discrete-Time Modified Leslie-Gower Model with Double Allee Effect and Its Stability and Bifurcation Analysis

Abstract

1. Introduction

2. The Fixed Points & Their Stability

3. Bifurcations

3.1. Period-Doubling Bifurcation

3.2. Neimark-Sacker Bifurcation

4. Numerical Simulations

4.1. Strong Allee Effect

4.1.1. Period-Doubling Bifurcation

4.1.2. Neimark-Sacker Bifurcation

4.2. Weak Allee Effect

4.2.1. Period-Doubling Bifurcation

4.2.2. Neimark-Sacker Bifurcation

5. Conclusion

Conflict of Interest

References:

Information

Guidelines

CP Initiatives

Follow CP