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.
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; \[\label{mainsys}\tag{1} \begin{cases} \dfrac{dx}{dt}= \dfrac{x\left(1-x\right)\left(x-\beta\right)}{\theta + x} – \xi \dfrac{x y}{\gamma + x}, \\ \dfrac{dy}{dt}= \rho y \left(1 – \dfrac{y}{\gamma + x}\right). \end{cases}\]
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 \(\beta\). A positive value of \(\beta\) corresponds to a strong Allee effect, while a negative value of \(\beta\) corresponds to a weak Allee effect. Singh \(\textit{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: \[\label{disc-Sys}\tag{2} \begin{cases} x_{n+1} = x_n + h x_n\left(\dfrac{\left(1 – x_n\right)\left(x_n – \beta\right)}{\theta + x_n} – \xi \dfrac{y_n}{\gamma + x_n}\right), \\ y_{n+1} = y_n + h\rho y_n\left(1 – \dfrac{y_n}{\gamma + x_n}\right). \end{cases}\]
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.
The following system of simultaneous equations is solved to find the fixed points. \[\label{main-model}\tag{3} \begin{cases} x\left(\dfrac{\left(1 – x\right)\left(x – \beta\right)}{\theta + x} – \xi \dfrac{y}{\gamma + x} \right) = 0,\\ \rho y\left(1 – \dfrac{y}{\gamma + x}\right) = 0. \end{cases}\]
If we define \(A = \frac{1 + \beta – \xi}{2}\) and \(B = \beta + \theta \xi\), then for any scenario, the system has four stationary points on the boundary \(\left\{E_1, E_2, E_3, E_4\right\} = \left\{\left(0, 0\right), \left(0, \gamma\right), \left(1, 0\right), \left(\beta, 0\right)\right\}\). Let \(E_{5^*} = \left(A, A + \gamma\right)\), \(E_5 = \left(\bar{x}_1, \bar{x}_1 + \gamma\right)\) and \(E_6 = \left(\bar{x}_2, \bar{x}_2 + \gamma\right)\), 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, \[\beta \in \left(1 + \xi – 2 \sqrt{\xi(1 + \theta)}, 1 + \xi + 2 \sqrt{\xi(1 + \theta)}\right).\] (ii) If \(\beta = 1 + \xi \pm 2 \sqrt{\xi(1 + \theta)}\), 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 – \xi )^2 < 4 \xi \theta\). (iii) If \(\beta < 1 + \xi – 2 \sqrt{\xi(1 + \theta)}\) or \(\beta > 1 + \xi + 2 \sqrt{\xi(1 + \theta)}\), i-e, \(A^2 – B > 0\), then \(A < \sqrt{A^2 – B}\) if \(\theta \xi < – \beta\), 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 + \beta > \xi\) and \(\beta + \theta \xi > 0\). In strong Allee effect it is always true. In the weak Allee effect, this is true if \(1 > \xi – \beta > 0\) and \(\theta \xi > – \beta\). 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.
At fixed points, the respective jacobian matrices are given by \[\begin{aligned} J_1 &= \left[ \begin{array}{cc} 1 – h\dfrac{\beta}{\theta } & 0 \\ 0 & 1 + h \rho \\ \end{array} \right], & J_2 &= \left[ \begin{array}{cc} 1 – h\left(\xi -\dfrac{\beta }{\theta }\right) & 0 \\ h \rho & 1 – h \rho \\ \end{array} \right], \\ J_3 &= \left[ \begin{array}{cc} 1 + h\dfrac{\beta – 1}{1 + \theta} & – h\dfrac{\xi}{1 + \gamma} \\ 0 & 1 + h \rho \\ \end{array} \right], & J_4 &= \left[ \begin{array}{cc} 1 – h \dfrac{\beta (\beta – 1)}{\beta +\theta } & – h\dfrac{\beta \xi}{\beta + \gamma} \\ 0 & 1 + h \rho \\ \end{array} \right], \\ J_5 &= \left[ \begin{array}{cc} 1 + h \bar{x}_1 \left(D – C\right) & -h\bar{x}_1 D \\ h \rho & 1 – h\rho \\ \end{array} \right], & J_6 &= \left[ \begin{array}{cc} 1 + h \bar{x}_2 \left(D – C\right) & -h\bar{x}_2 D \\ h \rho & 1 – h\rho \\ \end{array} \right], \end{aligned}\] where \(C = – \frac{2\left(A – \bar{x}_i\right)}{\bar{x}_i + \theta}\) and \(D = \frac{\xi}{\bar{x}_i + \gamma}\), \(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 = Tr\left(J_i\right)\), \(det_i = \left\vert J_i\right\vert\) and \(p_i(z) = z^2 – t_i z + det_i\), \(i = 1, …, 6\). Suppose \(z_1\), \(z_2\) are the roots of \(p(z)\). Then:
\(\left\vert z_{1,2}\right\vert < 1\) if and only if \(p(-1) > 0\) and \(det < 1\).
\(\left\vert z_{1,2}\right\vert > 1\) if and only if \(p(-1) > 0\) and \(det > 1\).
\(\left\vert z_1\right\vert < 1\) and \(\left\vert z_2\right\vert > 1\) or \(\left\vert z_1\right\vert > 1\) and \(\left\vert z_2\right\vert < 1\) if and only if \(p(-1) < 0\).
\(z_1 = -1\) with \(\left\vert z_2\right\vert \neq 1\) iff \(p(-1) = 0\) and \(t \neq -2, 0\).
\(z_{1,2}\) are complex and \(\left\vert z_{1,2}\right\vert = 1\) if and only if \(4det – t^2 > 0\) and \(det = 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 \(\left\vert z_1\right\vert < 1\) and \(\left\vert z_2\right\vert < 1\), so sink is locally asymptotically stable.
(ii)source if \(\left\vert z_1\right\vert > 1\) and \(\left\vert z_2\right\vert > 1\), so source is locally unstable.
(iii)saddle if \(\left\vert z_1\right\vert < 1\) and \(\left\vert z_2\right\vert > 1\) (or \(\left\vert z_1\right\vert > 1\) and \(\left\vert z_2\right\vert < 1\)).
(iv)non-hyperbolic if either \(\left\vert z_1\right\vert = 1\) or \(\left\vert z_2\right\vert = 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
\(\beta > 1\), saddle if \(\beta < 1\) and non-hyperbolic with one
eigenvalue \(1\) and the other not on
the unit circle, if \(\beta = 1\), i-e,
it may undergo transcritical or fold bifurcation. Whereas, \(E_{4}\) is unstable if \(\beta < 1\), saddle if \(\beta > 1\) and non-hyperbolic with one
eigenvalue \(1\) and the other not on
the unit circle, if \(\beta = 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 \(\beta =\left(1 + \xi \pm 2 \sqrt{(1 + \theta) \xi }\right)\), then \(E_{5^*}\) may undergo transcritical or fold bifurcation with eigenvalues, \(z_1 = 1\) and \(\left\vert z_2\right\vert \neq 1\) if and only if \(h \neq \frac{2}{\rho – A D}\) and \(\rho \neq A D\).
Proof. Suppose, \(\beta =\left(1 + \xi \pm 2 \sqrt{(1 + \theta) \xi }\right)\). Then, for strong Allee effect, \(E_{5^*}\) will always exist, while for weak Allee effect, it will exist if additionally \((1 – \xi )^2 < 4 \xi \theta\). The jacobian of the system (1), at \(E_{5^*}\) is \[J^* = \left( \begin{array}{cc} 1 + h A D & – h A D \\ h \rho & 1 – h \rho \\ \end{array} \right)\] and the eigenvalues are \(z_1 = 1\) and \(z_2 = 1 + h (a d – \rho)\). Thus, \(\left\vert z_2\right\vert \neq 1\) if and only if \(h \neq \frac{2}{\rho – A D}\) and \(\rho \neq A D\). ◻
Theorem 2. For the positive interior fixed points \(E_5\) and \(E_6\), suppose \(1 + \beta > \xi\), \(\beta + \theta \xi > 0\) and \(\beta < 1 + \xi – 2 \sqrt{\xi(1 + \theta)}\) and let, \[\begin{aligned} h_{1} &= \frac{(C – D) \bar{x}_1 + \rho }{\rho C \bar{x}_1} -\frac{\sqrt{\left((C – D) \bar{x}_1 + \rho\right)^2 – 4\rho C \bar{x}_1}}{\rho C \bar{x}_1}, \\ h_{2} &= \frac{(C – D) \bar{x}_1 + \rho }{\rho C \bar{x}_1} +\frac{\sqrt{\left((C – D) \bar{x}_1 + \rho\right)^2 – 4\rho C \bar{x}_1}}{\rho C \bar{x}_1}. \end{aligned}\] Then, for both strong and weak Allee effects, the following holds true.
(i) If \(h \in \left(0, h_1\right)\), then the fixed point is a sink.
(ii) If \(h \in \left(h_2, \infty\right)\), then the fixed point is a source.
(iii) The fixed is saddle if \(h \in \left(h_1, h_2\right)\).
(iv) The fixed point is non-hyperbolic with eigenvalues \(z_1 = -1\) and \(\left\vert z_2\right\vert \neq 1\) if \(h = h_1\) or \(h = h_2\) and \(h \neq \frac{2}{\rho – \bar{x}_1 (D-C)}\) and \(h \neq \frac{4}{\rho -\bar{x}_1 (D-C)}\).
(v) If, \[\begin{aligned} % \nonumber % Remove numbering (before each equation) \rho &\in \bigl(\left(D – C\right)\bar{x}_1, \infty\bigr) \cap \left(\left(\sqrt{C}-\sqrt{D}\right)^2 \bar{x}_1, \left(\sqrt{C}+\sqrt{D}\right)^2 \bar{x}_1 \right) \end{aligned}\] then, the fixed point is non-hyperbolic with complex conjugate eigenvalues, \(\left\vert z_1\right\vert = 1 = \left\vert z_2\right\vert\) if and only if \(h = \frac{\rho – \bar{x}_1 (D-C)}{C \rho \bar{x}_1} < \frac{4}{\rho -\bar{x}_1 (D-C)}\).Proof. The characteristic polynomial of the system at the stationary points is, \[p(z) = z^2 – t z + det,\] where, \(t = 2 + h \left(\bar{x}_1 (D – C) – \rho\right)\) and \(det = 1 + h \left(\bar{x}_1 (D – C) – \rho + h \rho \bar{x}_1 C\right)\). Then, we can use lemma 1, since for any \(h > 0\), \[p(-1) = \rho C \bar{x}_1 h^2 – 2 ((C – D) \bar{x}_1 + \rho) h + 4,\] which shows that if we define \(h_{1,2}\) as given, then \[\begin{cases} p(-1) > 0, & \hspace{5mm} h \in \left(0, h_1\right) \cap \left(h_2, \infty\right); \\ p(-1) < 0, & \hspace{5mm} h \in \left(h_1, h_2\right); \\ p(-1) = 0, & \hspace{5mm} h = h_1 \text{ or } h = h_2. \end{cases}\] Also, \[\begin{cases} det > 1, & \hspace{5mm} h \in \left(\frac{\rho + \bar{x}_1 (C – D)}{\rho C \bar{x}_1}, \infty\right); \\ det < 1, & \hspace{5mm} h \in \left(0, \frac{\rho + \bar{x}_1 (C – D)}{\rho C \bar{x}_1}\right); \\ det = 1, & \hspace{5mm} h = \frac{\rho + \bar{x}_1 (C – D)}{\rho C \bar{x}_1}. \end{cases}\] Then, using Lemmas 1 and 2, we arrive at the desired results. ◻
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} \Omega_{PD} &= \left\{\left(\beta, \theta, \xi, \rho, \gamma, h_0\right) \in \mathbb{R}^6_+ : h_0 = h_1 \text{ or } h_2, 1 + \beta > \xi, \hspace{2mm} h \neq \frac{2}{\rho – x_1 (D – C)}, h \neq \frac{4}{\rho – x_1 (D – C)} \right\}, \\ \Omega_{NS} &= \left\{\left(\beta, \theta, \xi, \rho, \gamma, h_0\right) \in \mathbb{R}^6_+ : 1 + \beta > \xi, h_0 = \frac{\rho – \bar{x}_1 (D-C)}{C \rho \bar{x}_1} < \frac{4}{\rho -\bar{x}_1 (D-C)} \right\}. \end{aligned}\]
Let \(H\) be the parameter for the mapping (2), such that \(\left\vert H\right\vert \lll 1\), to obtain bifurcations. We can write the perturbed mapping as \[\begin{aligned} & \left[ \begin{array}{c} x \\ y \\ \end{array} \right] \to \left[ \begin{array}{c} x + \left(h + H\right) x\left(\dfrac{\left(1 – x\right)\left(x – \beta\right)}{x + \theta} – \xi \dfrac{y}{x + \gamma} \right) \\ y + \left(h + H\right)\rho y\left(1 – \frac{y}{x + \gamma}\right) \\ \end{array} \right]. \end{aligned}\]
Define \(X = x – \bar{x}_1\) and \(Y = y – \bar{y}_1\). Using the series expansion method, we can write \[\label{mod-sys}\tag{4} \left[ \begin{array}{c} X \\ Y \\ \end{array} \right] \to \left[ \begin{array}{cc} a_{11} & a_{12} \\ a_{21} & a_{22} \\ \end{array} \right] \left[ \begin{array}{c} X \\ Y \\ \end{array} \right] + \left[ \begin{array}{c} f\left(X, Y, H\right) \\ g\left(X, Y, H\right) \\ \end{array} \right],\] where \(a_{11} = 1 + h b_3\), \(a_{12} = h b_4\), \(a_{21} = h \rho\), \(a_{22} = 1 – h \rho\), \(f\left(X, Y, H\right) = 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\left(\left\vert X + Y + H\right\vert^4\right) ,\) \(g\left(X, Y, H\right) = – h b_7 X^2 + 2 h b_7 X Y + \rho X H – h b_7 Y^2 + h a_4 X^3- 2 h a_4 X^2 Y – a_4 X^2 H\rho Y H + 2 b_7 X Y H + h a_4 X Y^2 – b_7 Y^2 H + O\left(\left\vert X + Y + H\right\vert^4\right)\) and \(b_1 = \theta (\bar{x}_1 + \theta) 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 + \theta)^2 a_1 – 1) – \bar{y}_1^2 a_3,\) \(b_4 = – \bar{x}_1 \bar{y}_1^2 a_3,\) \(b_5 = – h (\theta a_1 + a_3),\) \(b_6 = h a_3,\) \(b_7 = \frac{\rho}{\bar{y}_1},\) \(a_1 = \frac{(1 + \theta) (\beta + \theta)}{(\bar{x}_1 + \theta)^4},\) \(a_2 = \xi \frac{\bar{x}_1}{\bar{y}_1},\) \(a_3 = \frac{\gamma \xi}{\bar{y}_1^3},\) \(a_4 = \frac{\rho}{\bar{y}_1^2}.\)
Assume that the parameters \(\left(\beta, \theta, \xi, \rho, \gamma, h_0\right) \in \Omega_{PD}\). 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 \[\left[ \begin{array}{c} X \\ Y \\ \end{array} \right] = T_{PD}\left[ \begin{array}{c} u \\ v \\ \end{array} \right],\] where \(T_{PD}\) is an invertible matrix defined as, \[T_{PD} = \left[ \begin{array}{cc} a_{12} & a_{12} \\ -1-a_{11} & z_2 – a_{22} \\ \end{array} \right].\] Using \(T_{PD}\), we can transform our system to the one given below. \[\begin{aligned} \label{PD-mod-sys} &\left[ \begin{array}{c} u \\ v \\ \end{array} \right] \to \left[ \begin{array}{cc} -1 & 0 \\ 0 & z_2 \\ \end{array} \right] \left[ \begin{array}{c} u \\ v \\ \end{array} \right] + \left[ \begin{array}{c} f_{PD}\left(H, u, v\right) \\ g_{PD}\left(H, u, v\right) \\ \end{array} \right], \end{aligned}\tag{5}\] where, \[\begin{aligned} f_{PD}\left(H, u, v\right) =& 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\left(\left\vert H + u + v\right\vert^4\right),\\ g_{PD}\left(H, u, v\right) = &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\left(\left\vert H + u + v\right\vert^4\right), \end{aligned}\] and \[\begin{aligned} d_1 =& – c_1 h \left(a_{12}^2 \left(b_1 c_2-a_4\right)-\left(a_{11}+1\right){}^2 b_7 – a_{12} \left(a_{11}+1\right) \left(b_2 c_2+2 b_7\right)\right), \\ d_2=& -c_1 h \left(2 a_{12}^2 \left(a_4-b_1 c_2\right)+ a_{12} \left(a_{11}+a_{22}-z_2+1\right) \left(b_2 c_2+2 b_7\right) + 2 \left(a_{11}+1\right) b_7 \left(a_{22}-z_2\right)\right),\\ d_3=& c_1 \left(a_{12} \left(b_3 c_2+\rho \right)- \left(1+a_{11}\right) \left(b_4 c_2-\rho \right)\right),\\ d_4 =& c_1 \left(a_{12} \left(b_3 c_2+\rho \right)+ \left(z_2-a_{22}\right) \left(b_4 c_2-\rho \right)\right),\\ d_5 =& -c_1 h \left(a_{12} \left(a_{22}-z_2\right) \left(b_2 c_2+2 b_7\right)+ a_{12}^2 \left(b_7-b_1 c_2\right)+b_7 \left(a_{22}-z_2\right){}^2\right),\\ d_6=& a_{12} c_1 \left(-a_{12} \left(a_{11}+1\right) \left(b_6 c_2-2 a_4 h\right)+ a_{12}^2 \left(a_4 h+b_5 c_2\right)+\left(a_{11}+1\right){}^2 a_4 h\right),\\ d_7 =& a_{12} c_1 \left(+3 a_{12}^2 \left(a_4 h+b_5 c_2\right)+ a_4 \left(a_{11}+1\right) h \left(a_{11}+2 a_{22}-2 z_2+1\right)\right.\\ &\left. – a_{12} \left(2 a_{11}+a_{22}-z_2+2\right) \left(b_6 c_2-2 a_4 h\right) \right),\\ d_8=& c_1 \left(a_{12}^2 \left(b_1 c_2-a_4\right)-\left(a_{11}+1\right){}^2 b_7- a_{12} \left(a_{11}+1\right) \left(b_2 c_2+2 b_7\right) \right),\\ d_9=& -c_1 \left(2 \left(a_{11}+1\right) b_7 \left(a_{22}-z_2\right) +2 a_{12}^2 \left(a_4-b_1 c_2\right)+ a_{12} \left(a_{11}+a_{22}-z_2+1\right) \left(b_2 c_2+2 b_7\right)\right),\\ d_{10}=& a_{12} c_1 \left(3 a_{12}^2 \left(a_4 h+b_5 c_2\right)- a_{12} \left(a_{11}+2 a_{22}-2 z_2+1\right) \left(b_6 c_2-2 a_4 h\right) \right.\\ &\left. +a_4 h \left(a_{22}-z_2\right) \left(2 a_{11}+a_{22}-z_2+2\right)\right),\\ d_{11} =& -c_1 \left(a_{12} \left(a_{22}-z_2\right) \left(b_2 c_2+2 b_7\right)+ a_{12}^2 \left(b_7-b_1 c_2\right)+b_7 \left(a_{22}-z_2\right){}^2\right),\\ d_{12}=& a_{12} c_1 \left(-a_{12} \left(a_{22}-z_2\right) \left(b_6 c_2-2 a_4 h\right)+ a_{12}^2 \left(a_4 h+b_5 c_2\right)+a_4 h \left(a_{22}-z_2\right){}^2\right),\\ d_{13}=& – h c_1 \left((1+a_{11}+a_{12})^2 b_7 + a_{12} (-a_{12} b_1+b_2+a_{11} b_2) c_3\right),\\ d_{14}=& – h c_1 \left(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)\right),\\ d_{15} =& c_1(a_{12} (b_3 c_3 + \rho) + (1+a_{11}) (-b_4 c_3 + \rho)),\\ d_{16} =& c_1((-a_{22}+z_2) (b_4 c_3-\rho)+a_{12} (b_3 c_3+\rho)),\\ d_{17} =& – c_1 h \left(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) \right),\\ d_{18} =& a_{12} c_1 \left(a_{12} (a_{12} b_5-(1+a_{11}) b_6) c_3+ (1+a_{11}+a_{12})^2 a_4 h\right),\\ d_{19} =& a_{12} c_1 \left((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) \right.\\ &\left.- a_{12} (b_6 c_3 – 2 a_4 h) (2 + 2 a_{11} + a_{22} – z_2)\right),\\ d_{20} =& – c_1 \left((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)\right),\\ d_{21} =& – c_1\left(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)\right), \\ d_{22} =& a_{12} c_1 \left(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)\right. \\ &\left. + a_4 h (a_{22} – z_2) (2 + 2 a_{11} + a_{22} – z_2)\right), \\ d_{23} =& – c_1 \left(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)\right),\\ d_{24} =& a_{12} c_1 \left(- 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) \right),\\ X =&a_{12}\left(u + v\right), \\ 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 \(\left(0, 0\right)\), define, \[\begin{aligned} M_{PD} &= \left\{\left(H, u, v\right) : v = G_{PD}\left(H, u\right), \left\vert u\right\vert < \delta_{1}, \left\vert H\right\vert < \delta_{2}, \delta_1, \delta_2 \in \mathbb{R}^+, G_{PD}(0, 0) = 0, DG_{PD}(0,0) = 0 \right\}, \end{aligned}\] where \(G_{PD}\left(H, u\right) = -\frac{d_{15}}{1 + z_2} H u + O\left(\left\vert H + u\right\vert^3\right)\). Then the following map illustrates the dynamics restricted to \(M_{PD}\). \[\begin{aligned} \label{PD-CenterManifold} F : u &\mapsto – u + d_1 u^2 + d_3 H u – \frac{d_4 d_{15}}{1 + z_2} H^2 u + \left( d_8 – \frac{d_2 d_{15}}{1 + z_2}\right) H u^2 + d_6 u^3 + O\left(\left\vert H + u\right\vert^4\right). \end{aligned}\tag{6}\]
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 &=& \left.\left(\frac{\partial^2 F}{\partial H \partial u} + \frac{1}{2}\frac{\partial F}{\partial H}\frac{\partial^2 F}{\partial u^2}\right)\right\vert_{\left(0, 0\right)} = d_3 \neq 0 ,\\ l_2 &=& \left.\left(\frac{1}{6}\frac{\partial^3 F}{\partial u^3} + \left(\frac{1}{2}\frac{\partial^2 F}{\partial u^2}\right)^2\right)\right\vert_{\left(0, 0\right)} = 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.
Let \(\left(\beta, \theta, \xi, \rho, \gamma, h_0\right) \in \Omega_{NS}\). Then, variation of \(h\) around \(h_0\), gives emergence to Neimark-Sacker bifurcation. Instead of (4), we can write the perturbed system as \[\label{mod-sys-NS}\tag{7} \left[ \begin{array}{c} X \\ Y \\ \end{array} \right] \to \left[ \begin{array}{cc} c_{11} & c_{12} \\ c_{21} & c_{22} \\ \end{array} \right] \left[ \begin{array}{c} X \\ Y \\ \end{array} \right] + \left[ \begin{array}{c} \bar{f}\left(X, Y\right) \\ \bar{g}\left(X, Y\right) \\ \end{array} \right],\] where \(c_{11} = 1 + \left(h+H\right) b_3\), \(c_{12} = \left(h+H\right) b_4\), \(c_{21} = \left(h+H\right) \rho\), \(c_{22} = 1 – \left(h+H\right) \rho\), \[\begin{aligned} % \nonumber % Remove numbering (before each equation) \bar{f}\left(X, Y\right) &= \left(h+H\right) b_1 X^2 + \left(h+H\right) b_2 X Y + b_5 X^3 + b_6 X^2 Y + O\left(\left\vert X + Y + H\right\vert^4\right), \\ \bar{g}\left(X, Y\right) &= – \left(h+H\right) b_7 X^2 + 2 \left(h+H\right) b_7 X Y – \left(h+H\right) b_7 Y^2 + h a_4 X^3 – 2 h a_4 X^2 Y\\ &\quad + h a_4 X Y^2 + O\left(\left\vert X + Y + H\right\vert^4\right), \end{aligned}\] and the above mentioned coefficients are defined at the start of this section. Let \(p_{NS}(z)\) be the characteristic polynomial of matrix in (9), given by \[\label{ch-pol-NS}\tag{8} p_{NS}(z) = z^2 – P(H) z + Q(H),\] where \[P(H) = c_{11} + c_{22} \hspace{1mm} \text{ and } \hspace{1mm} Q(H) = c_{11} c_{22} – c_{21}c_{12}.\]
Since, \(\left(\beta, \theta, \xi, \rho, \gamma, h\right) \in \Omega_{NS}\), the roots of (8) are complex conjugate \(z_1\), \(z_2\) and \(\left\vert z_1\right\vert = \left\vert z_2\right\vert = 1\) It follows that \(z_{1,2} = \frac{P(H)}{2} \pm \frac{\iota}{2} \sqrt{4 Q(H) – P(H)^2}\) and \(1 = \left\vert z_1\right\vert = \left\vert z_2\right\vert = \sqrt{Q(H)}\). Also, \(\left(\frac{d\left\vert z_{1,2}\right\vert}{d H}\right)_{H=0} \neq 0\) if and only if \(\frac{(\theta +1) (\beta +\theta )}{(\bar{x}_1 + \theta)^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, \(\left(\beta, \theta, \xi, \rho, \gamma, h\right) \in \Omega_{NS}\), \(P(H)^2 -4 Q(H) < 0\), which implies that \(P(H)^2 < 4Q(H)\) or \(P(H) \in \left(-2,2\right)\). Thus, \(P(0) \neq -2, 0, 2\). Finally, \[\begin{aligned} P(0) &= h \left[\bar{x} \left(\frac{(\theta +1) (\beta +\theta )}{\left(\bar{x}+\theta \right)^2}+\frac{\xi }{\bar{y}}-1\right)-\rho \right] + 2, \end{aligned}\] which implies that \(P(0) \neq 1\) if and only if \(h \neq \frac{1}{\rho -\bar{x} \left(\frac{(\theta +1) (\beta +\theta )}{\left(\bar{x}+\theta \right)^2}+\frac{\xi }{\bar{y}}-1\right)}\).
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{4Q(0) – P(0)^2}\) and let \[T_{NS} = \left( \begin{array}{cc} 0 & 1 \\ S & R \\ \end{array} \right).\] Here, \(T_{NS}\) is an invertible matrix. Consider the transformation \[\left( \begin{array}{c} X \\ Y \\ \end{array} \right)=T_{NS} \left( \begin{array}{c} u \\ v \\ \end{array} \right).\] Using the above given transformation, we can write system (4) as, \[\begin{aligned} \left[ \begin{array}{c} u \\ v \\ \end{array} \right] &\to \left[ \begin{array}{cc} R & -S \\ S & R \\ \end{array} \right] \left[ \begin{array}{c} u \\ v \\ \end{array} \right] + \left[ \begin{array}{c} f_{NS}\left(u, v\right) \\ g_{NS}\left(u, v\right) \\ \end{array} \right]. \end{aligned}\tag{9}\]
\[\begin{aligned} % \nonumber % Remove numbering (before each equation) f_{NS}\left(u, v\right) =& 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\left(\left\vert u + v\right\vert^4\right) ,\\ g_{NS}\left(u, v\right) =& e_{7} u v + e_{8} v^2 + e_{9} u v^2 + e_{10} v^3 + O\left(\left\vert u + v\right\vert^4\right), \end{aligned}\] where, \[\begin{aligned} % \nonumber % Remove numbering (before each equation) e_1 &= – \left(h + H\right) b_7 S, \\ e_2 &= -(h+H) (b_2R+2 b_7(R-1)), \\ e_3 &= -\frac{(h+H) \left(R (b_1+b_2R)+b_7(R-1)^2\right)}{S}, \\ e_4 &= h a_4 S, \\ e_5 &= 2 a_4h (R-1)-b_6R, \\ e_6 &= \frac{a_4h (R-1)^2-R (b_5+b_6R)}{S}, \\ e_7 &= (h+H) b_2S, \\ e_8 &= (h+H) (b_1+b_2R), \\ e_9 &= b_6S, \\ e_{10} &= b_5+b_6R,\\ X &= v, \hspace{2mm} 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{aligned} \zeta_{NS} &= \left\{ \left[-\operatorname{Re} \left(\varpi_{20}\varpi_{11}\frac{\left(1 – 2z_{1}\right)z_{2}^{2}}{1 – z_{1}}\right) – \frac{1}{2} \left\vert\varpi_{11}\right\vert^{2} – \left\vert\varpi _{02}\right\vert^{2} + \operatorname{Re}(\varpi_{21}z_{2}) \right]\right\}_{H=0}, \end{aligned}\tag{10}\] where, \[\begin{aligned} % \nonumber % Remove numbering (before each equation) \varpi_{20} &= \frac{1}{8} \left[ \frac{\partial^2 f_{NS}}{\partial u^2} – \frac{\partial^2 f_{NS}}{\partial v^2} + 2\frac{\partial^2 g_{NS}}{\partial u \partial v} + \iota \left(\frac{\partial^2 g_{NS}}{\partial u^2} – \frac{\partial^2 g_{NS}}{\partial v^2} – 2\frac{\partial^2 f_{NS}}{\partial u \partial v}\right) \right] ,\\ \varpi_{11} &= \frac{1}{4} \left[ \frac{\partial^2 f_{NS}}{\partial u^2} + \frac{\partial^2 f_{NS}}{\partial v^2} + \iota \left(\frac{\partial^2 g_{NS}}{\partial u^2} + \frac{\partial^2 g_{NS}}{\partial v^2}\right) \right], \\ \varpi_{02} &= \frac{1}{8} \left[ \frac{\partial^2 f_{NS}}{\partial u^2} – \frac{\partial^2 f_{NS}}{\partial v^2} – 2\frac{\partial^2 g_{NS}}{\partial u \partial v} + \iota \left(\frac{\partial^2 g_{NS}}{\partial u^2} – \frac{\partial^2 g_{NS}}{\partial v^2} + 2\frac{\partial^2 f_{NS}}{\partial u \partial v}\right) \right],\\ \varpi_{21} &= \frac{1}{16} \left[\frac{\partial^3 f_{NS}}{\partial u^3} + \frac{\partial^3 f_{NS}}{\partial u \partial v^2} + \frac{\partial^3 g_{NS}}{\partial u^2 \partial v} + \frac{\partial^3 g_{NS}}{\partial v^3}+ \iota \left(\frac{\partial^3 g_{NS}}{\partial u^3} + \frac{\partial^3 g_{NS}}{\partial u \partial v^2} – \frac{\partial^3 f_{NS}}{\partial u^2 \partial v} – \frac{\partial^3 f_{NS}}{\partial v^3}\right) \right]. \end{aligned}\] Based on the prior analysis, we have the following result [41].
Theorem 4. If \(\frac{(\theta +1) (\beta +\theta )}{(\bar{x}_1 + \theta)^2} \neq 1\) and \(\zeta_{NS}\) 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 \(\left(\bar{x_1}, \bar{y_1}\right)\). Moreover, an attracting (resp., repelling) invariant closed curve bifurcate from the positive fixed point if \(\zeta_{NS} < 0\) (resp., \(\zeta_{NS} > 0\)).
In this section, we will use specific numerical values of parameters \(\left(\beta, \theta, \xi, \rho, \gamma\right)\) 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.
Example 1. Select \(\left(\beta, \theta, \xi, \rho, \gamma\right)\) \(=\) \(\left(0.1, 0.01, 0.05, 0.05, 0.3\right)\), \(h \in [2.4, 3.345]\) and the initial value \((x_0, y_0) = (1.2, 2)\). The jacobian matrix is, \[\left( \begin{array}{cc} 1 – 0.790244 h & -0.0379371 h \\ 0.05 h & 1 – 0.05 h \\ \end{array} \right),\] 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 \left(-1,1\right)\). 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 \(\left(0.943479, 1.24348\right)\) 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.
Example 2. Let us assume that for \(h \in [2.1, 2.34]\) and the initial value \((x_0, y_0) = \left(0.98, 0.6\right)\), we have the following parameter values, \(\left(\beta, \theta, \xi, \rho, \gamma\right)\) \(=\) \(\left(0.15, 0.06, 0.05, 1.1, 0.1\right)\). The jacobian matrix is, \[\left( \begin{array}{cc} 1 – 0.681574 h & -0.0451768 h \\ 1.1 h & 1 – 1.1 h \\ \end{array} \right),\] and \(z_1 = -0.985179 – 0.171532 \iota\) and \(z_2 = -0.985179 + 0.171532 \iota\) are the corresponding eigenvalues, for \(h = h_0 = 2.2285674489338163\). Both eigenvalues have modulus equal to \(1\). Also, \(h_0 \in \left(0, \frac{4}{\rho – x_1 (D – C)}\right) = \left(0, 2.24521\right)\), \(\left(\frac{d\left|z_{1}\right|}{d H}\right)_{H=0} = \left(\frac{d\left|z_{2}\right|}{d H}\right)_{H=0} = 1.78157 \neq 0\) and \(\zeta_{NS} = -7834.72 \neq 0\). Thus, the fixed point exhibits Neimark-Sacker bifurcation, at the equilibrium point \(\left(0.936652, 1.03665\right)\), as can be seen in Figure 5. Moreover, for \(h > h_0\) we have \(\zeta_{NS} < 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)\).
Example 3. Select \(\left(\beta, \theta, \xi, \rho, \gamma\right)\) \(=\) \(\left(-0.02, -\frac{\beta}{\xi}, 0.09, 0.01, 1.5\right)\), \(h \in [2.9, 3.92]\) and the initial value \((x_0, y_0) = (0.6, 2)\). The jacobian matrix is, \[\left( \begin{array}{cc} 1 – 0.678663 h & -0.0335146 h \\ 0.01 h & 1 – 0.01 h \\ \end{array} \right),\] 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 \left(-1,1\right)\). 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 \(\left(0.943479, 1.24348\right)\) 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.
Example 4. Let us assume that for \(h \in [3.23, 3.2]\) and the initial value \((x_0, y_0) = \left(0.6, 0.4\right)\), we have the following parameter values, \(\left(\beta, \theta, \xi, \rho, \gamma\right)\) \(=\) \(\left(-0.08, 0.06, 0.1, 0.9, 0.01\right)\). The jacobian matrix is, \[\left( \begin{array}{cc} 1 – 0.316267 h & -0.0987952 h \\ 0.9 h & 1 – 0.9 h \\ \end{array} \right),\] and \(z_1 = -0.355266 – 0.136097\iota\) and \(z_2 = -0.355266 + 0.136097\iota\) are the corresponding eigenvalues, for \(h = h_0 = 3.25592\). Both eigenvalues have modulus equal to \(1\). Also, \(h_0 \in \left(0, \frac{4}{\rho – x_1 (D – C)}\right) = \left(0, 3.28875\right)\), \(\left(\frac{d\left|z_{1}\right|}{d H}\right)_{H=0} = \left(\frac{d\left|z_{2}\right|}{d H}\right)_{H=0} = 0.832494 \neq 0\) and \(\zeta_{NS} = -2237.02 \neq 0\). Thus, the fixed point exhibits Neimark-Sacker bifurcation, at the equilibrium point \(\left(0.82, 0.83\right)\), as can be seen in Figure 9. Moreover, for \(h > h_0\) we have \(\zeta_{NS} < 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)\).
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.
The authors declare no conflict of interest.