Dynamical complexities and chaos control in a Ricker type predator-prey model with additive Allee effect

1.   Introduction

The intricate ecosystems found in the natural world are mostly the result of predator-prey interactions. Effective mathematical models are developed to address the problems in predator-prey interaction. A basic model for two species communities was formulated by Lotka in 1925 and Volterra in 1931. The predator's prey feeding rate is one essential component of predator-prey interaction. In literature, many interaction functions have been modeled, for example, Holling type Ⅰ, type Ⅱ or type Ⅲ ^[1], ratio-dependent ^[2], Crowley-Martin ^[3] and Beddington-DeAngelis models ^[4,5]. In recent years, the predator-prey model of Holling type Ⅱ with different biological situations has been considered in the literature. For instance, in ^[6], the authors considered such a model with prey refuge and showed that the model parameters have crucial importance in the coexistence, stability and oscillation of their considered model. In reality, a predator consumes prey at different rates, depending on the species. Leslie and Gower ^[7] proposed a predator-prey model in which the predator's growth function is given by the ratio of the sizes of predator and prey. Differential equations describe the evolution of systems in continuous and discrete time. Further, discrete models are more realistic if the populations are small or if births and deaths occur at discrete times or within certain time intervals, such as a generation ^[8,9].

To model real ecosystems, the term "Allee effect'' is unavoidable and meaningful. Thus, many researchers have paid much attention to studying the dynamics of predator-prey models with the Allee effect ^[10,11]. There are two types of Allee effects: weak and strong. In the weak Allee effect, for the low population size, the per capita growth rate is positive and there is no threshold to grow ^[12]. But, in the case of the strong Allee effect, the per capita growth is negative and there is a population threshold introduced by the Allee effect, so the population must surpass this threshold to grow ^[13]. The Allee effect can be introduced in different ways, namely, multiplicative ^[14] and additive ^[13]. Discrete-time predator-prey model with Allee effects was studied in ^[15,16], and some more recent developments in predator-prey models can be found in ^[17,18,19]. Based on the above discussions, the Allee effect term is incorporated into prey growth for this study.

The study of stability analysis, bifurcation behavior, namely, flip bifurcation and Neimark-Sacker bifurcation, and control of the chaos of the discrete-time model has long been an issue of significant interest for researchers ^[20,21]. Among the mathematics, physics and engineering communities, chaos control has been extensively studied ^[22], whereas the problem of controlling chaos in ecology has been given in ^[24]. In some nonlinear systems, there is a major route to chaos via period doubling. In ^[25], the authors proposed a new hybrid control method to control the chaotic orbits by using state feedback and parameter perturbation for the period-doubling bifurcation. The chaotic attractor can be converted into one of the possible attracting time-periodic motions by giving small perturbations to the system parameters ^[26]. Recently, the work related to controlling chaos in the predator-prey model has been of greater interest to researchers. For instance, a predator-prey model with the Allee effect and cannibalism in the discrete-time model has been found in the literature ^[27], where the authors showed that the considered system undergoes flip and Neimark-Sacker bifurcations. Further, they discussed the various methods to control the chaos. A discrete-time classical Lotka-Volterra model obtained by applying the method of a piecewise constant argument for differential equations has been proposed in ^[28]. The authors studied the stability, bifurcation and chaotic behavior of their proposed model. Also, they implemented hybrid control and the Ott-Grebogi-Yorke (OGY) method to control bifurcation and chaotic behavior.

Since it was an interesting topic for the researchers, we are motivated by the above literature. To the best of our knowledge, the study of stability, bifurcation and chaos control analysis for a discrete-time prey-predator model with the additive Allee effect has not been studied. In this article, we aim to study the discrete-time Holling type Ⅱ predator-prey model with the additive Allee effect in prey growth. The main highlights of this paper are listed as follows:

1) We use the method of piecewise constant arguments for the differential equations to obtain the discrete-time model. The model is a good representation of the populations' interaction, which has non-overlapping generations.

2) We study the impact of the Allee effect in terms of the existence and local stability properties of all positive equilibrium points of the discrete model.

3) We drive the conditions for the existence of flip and Neimark-Sacker bifurcations near the interior equilibrium point for the discrete model by taking the Allee parameter as a bifurcation parameter. Further, we discuss the properties of both bifurcations with the help of the center manifold (CM) theorem and normal form theory.

4) We utilize the state and OGY feedback and hybrid control methods to control the bifurcating and chaotic behavior of the discrete model. We provide extensive numerical simulations to show the rich dynamics, including periodic windows, invariant circles and chaos in the discrete model.

The rest of the article is arranged as follows: The discrete-time Holling type Ⅱ predator-prey model with the additive type Allee effect is given, and the existence and local stability of equilibrium points are discussed in Section 2. Furthermore, an analysis of the model's dynamics is presented, including the identification of key parameters and their influence on the system's behavior. In Section 3, the bifurcation behavior of the proposed model is thoroughly examined. Two types of bifurcations, such as flip and Neimark-Sacker bifurcations, are investigated to gain a comprehensive understanding of the system's response to parameter changes. In Section 4, the state and OGY feedback and hybrid control techniques for the considered predator-prey model are discussed. The effectiveness of these control approaches is highlighted, emphasizing their ability to regulate and stabilize system dynamics. To ensure the validity and reliability of the mathematical results, extensive numerical simulations are presented in Section 5. Lastly, in Section 6, a comprehensive conclusion and a brief discussion are presented.

2.   Model formation

First, we describe the continuous-time Holling type Ⅱ prey-predator model with the additive Allee effect ^[29], as follows:

where $H(t)$ and $P(t)$ are the population sizes of prey and predator respectively at time $t$; $r$ and $k$ represent the growth rate and carrying capacity of prey, respectively; $\gamma$ and $d$ represent the consumption rate of the predator on prey and per capita death rate of the predator, respectively. The most common type of functional response is Holling type Ⅱ, which is based on the idea that, if prey species is low, the rate of predation is proportional to prey density and is of the form $\frac{ a_1 H(t)}{1+a_1 a_2 H(t)}$, where $a_1$ and $a_2$ represent the predator capture rate and handling time per prey, respectively.

In general, the term $\frac{\alpha}{\beta+H(t)}$ is known as the additive Allee effect, where $\alpha$ and $\beta$ represent the severity of the Allee effect and the prey population size at which fitness is half of its maximum value, respectively. Many studies have been conducted to understand the dynamics of prey-predator system models with the additive Allee effect ^[30,31,32].

To reduce the complexity, and, for our convenience, after redefining the parameters, the model (2.1) takes the form

where $H(t)\leq0$, $P(t)\leq 0$, $a = r$, $b = \frac{r}{k}$, $c = \frac{a_1}{a_1 a_2}$, $e = \frac{1}{a_1 a_2}$ and $f = \frac{\gamma a_1}{a_1 a_2}$. Note that, in the sense of ^[33], if $\alpha < a \beta$ ($\alpha > a \beta$), the Allee effect is the weak (strong) one.

Next, based on the appropriate modification of overlapping generations, one can get the difference equations for modeling a population with non-overlapping generations. Euler approximations, non-standard finite difference schemes ^[34] and piecewise constant arguments ^[35,36] are the most common methods to obtain discrete versions from the continuous models. In this way, differential equations with piecewise constant arguments have been useful. Let us assume that the populations have no overlap between the successive generations and the population growth occurs in discrete steps $t\in[n, n+1)$, $n = 0, 1, 2, 3, \dots$. Also, consider that the variables and constants in (2.2) change in regular time intervals and incorporate this idea; then the corresponding discrete-time model for (2.2) is obtained by the method of piecewise constant arguments for differential equations, as follows:

where $H(n+1)$ and $P(n+1)$ denote the populations in generation $n+1$ that are related to the sizes $H(n)$ and $P(n)$ of the populations in the preceding generation $n$ and $a$, $b$, $c$, $d$, $e$, $f$, $\alpha$ and $\beta$ are all positive constants.

Note, in the absence of a predator $P$ and the Allee effect, the above model (2.3) becomes a one-dimensional model similar to the Ricker model ^[37], namely,

This system represents the relationship between the current and previous population size $H$. In a logistic map, the population's rate of per capita growth falls as the population size approaches a limit set by the availability of resources. On the Ricker map, population growth is almost exponential; however, as population size increases, the instantaneous growth rate declines linearly; and, eventually, population size reaches a plateau and oscillates around a mean.

The main aim of this paper is to study the local stability, bifurcation behavior and various chaos control analyses for the discrete time prey-predator model (2.3) with the additive Allee effect in the case of the weak one in $\mathbb{R} = \left\{ (H, P)|H > 0, P > 0\right\}$.

2.1.   Equilibria

We need to solve the following equations to find the equilibrium points of the model (2.3):

By direct substitution, from (2.4), we have three positive equilibrium points: The trivial equilibrium point $(0, 0)$, the predator-free equilibrium point $(\bar{H}, 0)$, where $\bar{H}$ is the positive root of the equation

and the interior equilibrium point $(H^*, P^*)$, which is given by

Then from ^[29], we have the following lemma.

Lemma 1. Consider $\alpha < a \beta$; then, the model (2.3) has a unique boundary equilibrium $(\bar{H}, 0) = \left(\frac{(a-b \beta)+\sqrt{(a-b \beta)^2+4 b(a \beta+\alpha)}}{2 b}, 0\right)$. Also, the model (2.3) has a unique interior equilibrium $(H^*, P^*) = \left(\frac{d e}{f-d}, \frac{(a-b H^*)(\beta+H^*)(e+H^*)-\alpha (e+H^*)}{c(\beta+H^*)}\right)$ if $\alpha < (a-b H^*)(\beta+H^*)$ and $\max\left(\beta, \frac{d e}{f-d}\right) < \frac{a}{b} < \beta+\frac{d e}{f-d}$.

2.2.   Local stability analysis

Next, the Jacobian matrix for the model (2.3) at arbitrary equilibrium point $(H, P)$ is calculated to investigate the local stability property:

where $A_1 = \exp \left(a-b H- \frac{\alpha}{\beta+H}- \frac{c P}{e+H}\right)$ and $A_2 = \exp \left(\frac{f H}{e+H}-d\right)$.

The Jacobian matrix of (2.3) derived at $(0, 0)$ is

and we have the Jacobian matrix at $(\bar{H}, 0)$ as

where $a_1 = \left(1-b \bar{H}+\frac{\alpha \bar{H}}{(\beta+\bar{H})^2}\right)\exp \left(a-b \bar{H}-\frac{\alpha}{\beta+\bar{H}}\right)$, $a_2 = -\left(\frac{c \bar{H}}{e+\bar{H}}\right)\exp \left(a-b \bar{H}-\frac{\alpha}{\beta+\bar{H}}\right)$ and $a_3 = \exp \left(\frac{f \bar{H}}{e+\bar{H}}-d\right)$.

The Jacobian matrix of (2.3) evaluated at $(H^*, P^*)$ is

The characteristic equation of $J^*$ is given by

where

From ^[38], let $\lambda_1$ and $\lambda_2$ be the eigenvalues of the Jacobian matrix (2.6) for some arbitrary equilibrium point $(H, P)$. We recall some topological classifications of the equilibrium points. The equilibrium point $(H, P)$ is a sink (locally asymptotically stable) if $|\lambda_1| < 1$ and $|\lambda_2| < 1$. $(H, P)$ is a source (locally unstable) if $|\lambda_1| > 1$ and $|\lambda_2| > 1$. $(H, P)$ is a saddle if $|\lambda_1| > 1$ and $|\lambda_2| < 1$ (or $|\lambda_1| < 1$ and $|\lambda_2| > 1$). And, $(H, P)$ is non-hyperbolic if either $|\lambda_1| = 1$ or $|\lambda_2| = 1$. As we know, the eigenvalue of the Jacobian matrix plays an important role in establishing the stability properties of the equilibrium points. Using the eigenvalues, we discuss the stability of the equilibrium points in the following lemmas:

Lemma 2. The trivial equilibrium $(0, 0)$ is always a saddle and the predator-free equilibrium $(\bar{H}, 0)$ is a

(1) sink if $|a_1| < 1$ and $|a_3| < 1$;

(2) source if $|a_1| > 1$ and $|a_3| > 1$;

(3) saddle if $|a_1| > 1$ and $|a_3| < 1$ (or $|a_1| < 1$ and $|a_3| > 1$).

Lemma 3. The coexisting equilibrium point

(1) $(H^*, P^*)$ is a sink if $D < 1$ and $|T| < D+1$;

(2) $(H^*, P^*)$ is a source if $D > 1$ and $|T| < D+1$ or $|T| > D+1$;

(3) $(H^*, P^*)$ is a saddle if $0 < |T|+D+1 < 2 |T|$;

(4) $(H^*, P^*)$ is non-hyperbolic if $|T| = |D+1|$, or $D = 1$ and $|T|\leq 2$.

3.   Bifurcations

Previously, we discussed the local stability properties of the equilibrium point $(H^*, P^*)$. Now, we study the bifurcation behavior, namely, flip and Neimark-Sacker bifurcations, of the model (2.3) at $(H^*, P^*)$ by varying the parameter $\alpha$ and keeping all model parameters fixed. It should be noted that further analyses can hold for other parameters also. Also, stability properties of the existing period-two orbit and the invariant closed curve are discussed with the help of the CM theorem and normal form theory as in ^[38,39]. The conditions for the model (2.3) to undergo both types of bifurcation are discussed below.

We know that the model can undergo a flip bifurcation if one of the eigenvalues of $J^*$ is $-1$ and the other not $1$ or $-1$. Hence, we assume one eigenvalue as $-1$, then, from (2.10),

where $\alpha = \alpha_f$ is the critical point of $\alpha$ that satisfies (3.1), and it is necessary for the flip bifurcation near the equilibrium point $(H^*, P^*)$ for the model (2.3).

Let $\Omega_F = \left\{ (a, b, c, d, e, f, \alpha, \beta):\alpha = \alpha_f, a, b, c, d, e, f, \beta > 0\right\}$; at the critical parameter value $\alpha = \alpha_f$, the model exhibits flip bifurcation at $(H^*, P^*)$, when $\alpha = \alpha_f$ varies in the neighborhood of $\Omega_F$.

Also, the model (2.3) admits the Neimark-Sacker bifurcation if the Jacobian matrix $J^*$ has complex conjugate eigenvalues with a modulus value of one, and also satisfies the following conditions

which gives

Now $A_2(\alpha) = 0$ gives

Let $\Omega_{NS} = \left\{ (a, b, c, d, e, f, \alpha, \beta):\alpha = \alpha_h, A_1(\alpha) < 0, a, b, c, d, e, f, \beta > 0\right\}$; at the critical parameter value $\alpha = \alpha_h$, the model (2.3) can exhibit Neimark-Sacker bifurcation near the equilibrium point $(H^*, P^*)$, when $\alpha = \alpha_h$ changes in the neighborhood of $\Omega_{NS}$. Then, we conclude the above results in the following theorem:

Theorem 1. (i) If $\Omega_F = \left\{ (a, b, c, d, e, f, \alpha, \beta):\alpha = \alpha_f, a, b, c, d, e, f, \beta > 0\right\}$ exists and satisfies (3.1) at critical parameter value $\alpha = \alpha_f$, then the model (2.3) undergoes flip bifurcation.

(ii) If $\Omega_{NS} = \left\{ (a, b, c, d, e, f, \alpha, \beta):\alpha = \alpha_h, A_1(\alpha) < 0, a, b, c, d, e, f, \beta > 0\right\}$ exists and satisfies (3.2) at critical parameter value $\alpha = \alpha_h$, then the model (2.3) undergoes Neimark-Sacker bifurcation.

3.1.   Flip bifurcation

Now we investigate the possible flip bifurcation of the model (2.3) at $(H^*, P^*)$. From the above discussion the model (2.3) undergoes flip bifurcation at $\alpha = \alpha_f$, where $\alpha$ varies in $\Omega_{F}$. Given a perturbation $|\alpha_1| \ll 1$ of $\alpha_f$, then perturbation of model (2.3) is described as

Next by shifting $(H^*, P^*)$ to the origin of (2.3) by using the transform $u_n = H_n-H^*$ and $z_n = P_n-P^*$, we have

where

Let us assume that the eigenvalues are $\lambda_1 = -1$ and $\lambda_2 = 3-b H^*+\frac{c H^* P^*}{(e+H^*)^2}+\frac{\alpha_f H^*}{(\beta +H^*)^2}$ for the matrix $J$ with $|\lambda_1| = 1$, $|\lambda_2|\neq 1$.

Next, we construct the non-singular matrix $L$ as follows

and we use the translation $\left(\begin{array}{l}u_n\\v_n\end{array}\right) = L\left(\begin{array}{l}U_n\\V_n\end{array}\right)$; then, (3.4) can be written as

where

and

Now, let us assume $G^c$ to be the CM; then, by using the CM theorem, we approximate the CM $G^c$ of (3.5) at the origin for small changes in $\alpha_1 = 0$:

By substituting (3.5) on both sides of $V_n = h(U_n, \alpha_1)$, we have

where

Accordingly, on the CM $G^c$ at the origin, we have

where

Moreover, the map confined to the CM $G^c(0, 0)$ takes the form

where

Finally, from ^[40], we define $\Delta_1$ and $\Delta_2$ as follows:

Therefore, we have the following findings about flip bifurcation from the aforementioned study.

Theorem 2. If $\Delta_1\neq0$ and $\Delta_2\neq0$, the model (2.3) exhibits a flip bifurcation at $(H^*, P^*)$ while changing the parameter $\Delta$ nearby $\Delta_1$. Moreover, if $\Delta_2 > 0$ (or $\Delta_2 < 0$) then the existing period two orbits from $(H^*, P^*)$ are stable (or unstable).

3.2.   Neimark-Sacker bifurcation

Now we analyze the properties of possible Neimark-Sacker bifurcation around $(H^*, P^*)$ for model (2.3) if, for instance, (3.2) holds for some $\alpha_h$. Given a perturbation $|\alpha_2| \ll 1$ of $\alpha_h$, perturbation of model (2.3) is described as

Let us use the transforms $u_n = H_n-H^*$ and $v_n = P_n-P^*$ and shift $(H^*, P^*)$ to $(0, 0)$; the model (3.9) takes the form

Then, the Taylor expansion of (3.10) at the origin, up to order three, is

where

The characteristic polynomial equation associated with the linearized system of (3.10) at the origin:

where

with $A_1 = \exp \left(a-b H^*-\frac{(\alpha_h+\alpha_2)}{\beta+H^*}-\frac{c P^*}{e+H^*}\right)$ and $A_2 = \exp \left(\frac{f H^*}{e+H^*}-d\right)$. Now, the roots of (3.12) are expressed as a pair of complex conjugates:

Since $(a, b, c, d, e, f, \alpha, \beta)\in \Omega_{NS}$, we have that $|\lambda_{1, 2}| = \sqrt{q_2(\alpha_2)}$ and

Further, we assume that $q_1(0) = -2+b H^*-\frac{\alpha_h H^*}{(\beta+H^*)^2}-\frac{c H^* P^*}{(e+H^*)^2}\neq0, -1$ and (2.6) implies that $q_1(0)\neq\pm\; 2, 0, -1$, which means that $\lambda^k_1, \lambda^k_2\neq 1$ for $k = 1, 2, 3, 4$ when $\alpha_2 = 0$. We only require that $q_1(0)\neq 0, 1$, which we can attain if it satisfies

Letting $\alpha_2 = 0$, $\xi = -\frac{q_1(0)}{2}$ and $\theta = \frac{\sqrt{4 q_2(0)-q_1^2(0)}}{2}$, we construct the non-singular matrix

and we use the translation $\left(\begin{array}{c}u_n\\v_n\end{array}\right) = L\left(\begin{array}{c}U_n\\V_n\end{array}\right)$; thus, the model (3.9) takes the form

where

Next, we require the non zero quantity $a^*$ to ensure that (3.9) admits Neimark-Sacker NS bifurcation.

where

Finally, from ^[38], we can state the following findings:

Theorem 3. If (3.13) and (3.14) hold and the quantity $a^*$ is non zero then the model (3.1) admits Neimark-Sacker bifurcation at $(H^*, P^*)$ when $\alpha_h$ changes in the neighborhood of $\Omega_{NS}$. Additionally, the quantity $a^* < 0$ (or resp. $a^* > 0$); then, the stable (or resp. unstable) invariant closed curve from $(H^*, P^*)$ starts to bifurcate.

4.   Chaos control

The theory of chaos control and bifurcation is known to be an essential area of current research. In particular, it has significant uses in biological and engineering sciences. Its features have emerged in ecological models. The nature of discrete-time population models is mostly unstable and chaotic in contrast to the continuous-time models. It is known that the logistic difference equation (logistic map) and the Ricker model ^[37] are the discrete analogs of the logistic differential equation ^[9]. The logistic map is derived from the logistic differential equation by using Euler approximation, taking the finite difference quotient as $\Delta t = 1$, i.e., the time step is one generation. The sequence of periodic orbits starts oscillate from stable point like period 2, period 4, period 8, .., cycles and finally enter into the chaotic regime for many difference equations; this collection of such simple difference equations have been reported in ecology ^[41]. The Neimark-Sacker bifurcation is similar to the Hopf bifurcation in the continuous system. Similarly, there has been a lot of research work on the presence of complex and chaotic nature in the two-dimensional discrete prey-predator system, this chaotic dynamics absent in the continuous time system ^[42]. Further, it is important to introduce the appropriate methods for controlling chaos to prevent unexpected scenarios.

First, we use the state feedback control method as in ^[43,44] to control the chaotic system. For the model (2.3) we consider the following corresponding controlled model:

where $u(H_n, P_n) = h_1(H_n-H^*)+h_2(P_n-P^*)$ is the feedback controlling force with feedback gains $h_1$ and $h_2$, and $(H^*, P^*)$ is the unique positive interior equilibrium point of (4.1). Consider the controlled model (4.1) whose Jacobian matrix at $(H^*, P^*)$ is given as

where

Then we have the characteristic polynomial of $J(H^*, P^*)$ as

Let $\lambda_1$ and $\lambda_2$ be the eigenvalues of (4.3); then, we obtain

Then, the lines of marginal stability are derived by solving $\lambda_1 = \pm 1$ and $\lambda_1 \lambda_2 = 1$, assuming that $|\lambda_{1, 2}| < < 1$. If $\lambda_1 \lambda_2 = 1$, using the second part of (4.4), we have

Next, if $\lambda_1 = 1$, and using (4.4), we get

Also, if $\lambda_1 = -1$, and using (4.4), we obtain

Then the triangular region enclosed by lines $L_1$, $L_2$, $L_3$ have stable eigenvalues.

Next, we study a chaos controlling technique based on the OGY control method as in ^[26]. By taking $\alpha$ as the control parameter, we rewrite the model (2.3) as follows:

Further, $\alpha$ must lie in some interval $|\alpha-\alpha_0| < \delta$ with $\delta > 0$, and $\alpha_0$ denotes the nominal value, for which the model (2.3) has chaotic dynamics. Now, we apply the state feedback control strategy to shift the trajectory to the desired state. Let $(H^*, P^*)$ be the unstable equilibrium point of the model (2.3) in the chaotic state due to Neimark-Sacker bifurcation; then, model (2.3) can be approximated near $(H^*, P^*)$ by the following linear map:

where

It is easy to see that model (4.8) is controllable provided that the following matrix

has rank 2. Furthermore, assume that $\frac{-H^*}{\beta+H^*}\neq 0$ and $\frac{\partial f_2(H^*, P^*, \alpha_0)}{\partial H_n}\neq 0$; then, model (4.8) is controllable. Next, we assume that $[\alpha-\alpha_0] = -K \left[\begin{array}{c}H_n-H^*\\ P_n-P^*\end{array}\right]$, where $K = \left[\begin{array}{cc}s_1 & s_2\end{array}\right]$; then, model (4.8) can be written as

Then, the controller model is given by

Moreover, equilibrium point $(H^*, P^*)$ is locally asymptotically stable if and only if both eigenvalues of the matrix $A-B K$, say, $\lambda_1$ and $\lambda_2$, lie in an open unit disk. Then, the matrix $A-B K$ can be written as

where

Then, the characteristic polynomial of (4.13) is given by

Then, the lines of marginal stability can be obtained, that is,

Then, stable eigenvalues lie within the triangular region in the $s_1, s_2$-plane bounded by the straight lines $L_1, L_2, L_3$ for the model's parametric values.

Next, we apply a hybrid control feedback methodology ^[21,25,45] for the control of flip and Neimark-Sacker bifurcations of the model near the equilibrium point $(H^*, P^*)$; then, the controlled model can be written as

where $0 < \epsilon < 1$ is the controlled strategy of a combination of both feedback control and parameter perturbation. The Jacobian matrix evaluated for the model (4.18) at $(H^*, P^*)$ is given by

Note, one can select appropriate values for $\epsilon$ to ensure that all eigenvalues of the above matrix satisfy $|\lambda_{1, 2}| < 1$.

5.   Numerical simulations

In this section, we provide a numerical simulation to ensure our theoretical analysis. The numerical iteration used to validate our results is similar to the processes described in ^{[27,28,36,44,46]}. The mathematical software Matlab (2018a) has been used to plot the numerical graphs. The computational approach provides fascinating details about how ecological systems function. Analytical approaches are exact and wide, but they only operate with the most fundamental models. The parameter values used here are only for illustrative purposes; neither experiments nor field research was used to obtain them ^[29,30,31]. Here, all parameters are considered positive values for which the model (2.3) has at least one positive interior equilibrium point (coexistence of both populations). Also, the chosen parameter values for the considered model (2.3) are biologically feasible.

Next, the conditions for the occurrence of two possible bifurcation behaviors of the model (2.3) were examined by varying the Allee parameter, namely, flip and Neimark-Sacker bifurcation verified numerically. It is also shown with the help of bifurcation diagrams. For this, we take two sets of parameter values assuming $\alpha < a \beta$, and described them in the following cases:

Case (i): Show the model (2.3) exhibits flip bifurcation with the fixed parametric values of $a = 5.5$, $b = 5.5$, $c = 20$, $d = 0.1$, $e = 10$, $f = 2$, $\beta = 0.5$, while varying $\alpha = (0, 2]$. We obtained the unique interior positive equilibrium points $(H^*, P^*)$ for the various values of $\alpha$ given in Table 1. At the critical point $\alpha = \alpha_f = 1.5833188$, the model (2.3) has the unique equilibrium point $(H^*, P^*) = (0.526316, 0.559233)$, and (2.10) becomes

which shows that, $\lambda_1, \lambda_2 = 0.949529, -1$ are the eigenvalues of (2.9), Then, it satisfies the condition. It ensures the occurrence of flip bifurcation near the interior equilibrium point for the model (2.3). Also, for $\alpha = 0.20$, the model (2.3) shows chaotic behavior, and by further increasing, it shows period-16 at $\alpha = 0.37$, period-8 at $\alpha = 0.43$, period-4 at $\alpha = 0.50$, period-2 at $\alpha = 1.25$, and stable equilibrium point at $\alpha = 1.70$. The different states of the model (2.3) are plotted for different values of $\alpha$ in the time plot and phase portrait in Figures 1 and 2, respectively. It indicates that the model becomes stable from chaos by the reverse period-doubling phenomenon. So, the existence of flip bifurcation for the model (2.3) is shown in Figure 3(a). Also, the existence of chaos is confirmed by plotting the corresponding largest Lyapunov exponent in Figure 3(b). Further, from (3.7) and (3.8), we obtain $\Delta_1 = 0.415148\neq 0$, $\Delta_2 = -2.22965\neq 0$ and the property of the existing flip bifurcation which is illustrated in Theorem 2.

Case (ii): Show that the model (2.3) exhibits Neimark-Sacker bifurcation with the fixed parametric values of $a = 3.71$, $b = 4.8$, $c = 10$, $d = 1$, $e = 10$, $f = 29$, $\beta = 0.5$, and varying values of $\alpha = (0, 0.5]$. We obtained the unique interior positive equilibrium points $(H^*, P^*)$ for the various values of $\alpha$ given in Table 2. At the critical point $\alpha = \alpha_h = 0.4135275$, the model (2.3) has the unique equilibrium point $(H^*, P^*) = (0.357143, 1.56731)$ and the characteristic Eq (2.10) becomes

which shows that $\lambda_1, \lambda_2 = 0.269458\pm 0.963012 i$ are the eigenvalues of (2.9). Then, it satisfies the condition (3.2). It ensures the occurrence of Neimark-Sacker bifurcation near the interior equilibrium point for the model (2.3). Also, for $\alpha = 0.015$, the model (2.3) shows chaotic orbit. Further, for $\alpha = 0.045$, the result is 13-period; for $\alpha = 0.048$, it is 13 invariant circles; for $\alpha = 0.130$, it is 9-period; for $\alpha = 0.410$, it is a closed invariant circle and there is a asymptotically stable equilibrium point at $\alpha = 0.450$. The different nature of the model around the interior equilibrium point for the model (2.3) is shown in the time plots and phase portraits in Figures 4 and 5, correspondingly. The figures show the model undergoes Neimark-Sacker bifurcation. The bifurcation diagram and its largest Lyapunov exponent are shown in Figure 6(a), (b), correspondingly. Since $a^* = -0.151225$ from (3.16), its results are explained in Theorem 3.

In order to show the effect of Allee parameter $\alpha$ on the growth rate parameter $a$ for the values in both Cases (i) and (ii), we have plotted the largest Lyapunov exponents in Figure 7(a), (b). This figures clearly shows the chaotic and stable dynamics of the proposed model.

5.1.   Chaos control

For the particular choice of parameter values $a = 3.71$, $b = 4.8$, $c = 10$, $d = 1$, $e = 10$, $f = 29$, $\beta = 0.5$, and $\alpha = 0.015$, we get the positive unstable equilibrium point $(H^*, P^*) = (0.357143, 2.0488647)$ for the model (2.3) and for its phase portrait (see Figure 5(a)). So, we chose to use the state feedback control method and want to shift the unstable equilibrium point into stable dynamics. For this, taking $\alpha_0 = 0.015$ and the corresponding controlled model yields

where $u(H_n, P_n) = u(H_n, P_n) = h_1(H_n-H^*)+h_2(P_n-P^*)$, and $h_1$ and $h_2$ are feedback gains. Furthermore, the Jacobian matrix of the controlled model is given as follows

Then, we obtain the characteristic equation as

Moreover, the lines of marginal stability for the model (5.1) are computed as

Then, the triangular region enclosed by the lines $L_1, L_2, L_3$ has stable eigenvalues for the controlled model (5.1), as shown in Figure 8(a). Clearly, if $h_2 = -0.3$, then $(H^*, P^*) = (0.357143, 2.0488647)$ is locally stable if and only if $h_2\in[-0.642, 0.4853]$. Taking $h_1 = -0.3$ and $h_2 = 0.3$, then the time-series plots for the controlled model (5.1) are shown in Figure 8(b), (c).

Second, for $a = 3.71$, $b = 4.8$, $c = 10$, $d = 1$, $e = 10$, $f = 29$, $\beta = 0.5$ and $\alpha = 0.015$, the model (2.3) yields $(H^*, P^*) = (0.357143, 2.0488647)$ and its phase portrait is shown in Figure 5(a). So, we chose to use the OGY feedback control method to move the strange chaotic dynamics into the stable dynamics. For this, taking $\alpha = 0.015$ and the corresponding controlled model yields

where $s_1$ and $s_2$ are feedback gains. Furthermore, in this case, we have

and the Jacobian matrix of the controlled model is given as follows

Then we obtain the characteristic equation

Moreover, the lines of marginal stability for the model (5.2) are computed as

Then the triangular region enclosed by the lines $L_1, L_2, L_3$ has stable eigenvalues for the controlled model (5.2) as shown in Figure 9(a). Clearly, if $s_1 = 0.1$, then $(H^*, P^*) = (0.357143, 2.0488647)$ is locally stable if and only if $s_2\in[0.55128, 0.82758]$. Taking $s_1 = 0.1$ and $s_2 = 0.7$, then the time-series plots for the controlled model (5.2) are depicted in Figure 9(b), (c).

Finally, we discuss the hybrid control strategy for some parametric values. In this case, the controlled model can be written as

Moreover, the Jacobian matrix of the model (5.2) at $(H^*, P^*) = (0.526316, 1.26863)$ is given by

The eigenvalues of the above matrix lie inside the open disk if and only if $0 < \epsilon < 0.773477$. For $\epsilon\in(0, 1]$ and $(H^*, P^*) = (0.526316, 1.26863)$, the bifurcation diagram for the controlled model (5.3) are depicted in Figure 10(a), (b).

Remark 1. The model (2.3) is the discrete counterpart of the continuous-time model (2.1) investigated in ^[29]. The method of using a piecewise constant argument for the differential equations has been is adopted to obtain the discrete time model, similar to those in ^[36,46]. In the absence of predators, the resultant model transforms into the classical Ricker stock recruitment model ^[37]. It is worth mentioning here that this type of Ricker-type model (2.3) has not been investigated previously. For instance, in ^[47], the authors studied a basic discrete predator model with the effect of using Ricker map for the prey growth by replacing the classical logistic map; they showed the existence of flip and Neimark-Sacker bifurcations. A Ricker-type predator-prey model with hunting cooperation was investigated by the authors of ^[48], who studied the existence of discrete Hopf bifurcation. The authors of ^[29] demonstrated that the model (2.1) experiences Hopf bifurcation, while the model (2.3) undergoes flip and Neimark-Sacker bifurcations and also exhibits chaotic dynamics. The direction and stability properties of both bifurcations can be discussed by using the results form ^[38,39]. Also, the bifurcation and chaos control analyses have been carried out using the methods described in ^[27,28,44].

6.   Discussion and conclusions

This paper deals with the discrete-time predator-prey model of Holling Ⅱ-type interaction with a weak additive Allee effect. In the first step of this study, we derived the discrete-time model from the corresponding continuous model in ^[29] by using the method of piecewise constant arguments for the differential equations ^[35]. This discrete form of Lotka-Volterra equations is due to May ^[9]. The model (2.3) without the Allee effect or predator population has been reduced to a one-dimensional model similar to the model in ^[37]. The discrete model (2.3) considered in this present work is a good representation of population interaction, as it has a non-overlapping generation, i.e., the birth and death rates occur in certain time intervals. From a biological point of view, the positive equilibrium points for the continuous-time model described in ^[29] and the discrete-time model (2.3) in this present study are the same, as we showed in Lemma 1. For model (2.3), we showed that the density of the prey population remains unchanged, and that the density of the predator population decreases by increasing the Allee parameter $\alpha$; see Tables 1 and 2. The local stability of a coexisting equilibrium $(H^*, P^*)$ is described via Lemma 3, which ensures the long-time survival of both species with the impact of $\alpha$. In the contrasting continuous-time model ^[29] of (2.3), for a small value the of Allee parameter, the model exhibits a periodic solution and becomes locally asymptotically stable via Hopf bifurcation, when the Allee parameter crosses some critical value. Moreover, this is similar to the Neimark-Sacker bifurcation for the discrete model (2.3) in this present study.

Next, we derived the conditions for the existence of flip and Neimark-Sacker bifurcations for model (2.3) around the positive interior equilibrium point by taking the Allee parameter as a bifurcation parameter. We showed the properties of both bifurcations with the help of the CM theorem and normal form theory in Theorems 2 and 3. In the simulation part, we verified all the analytical results numerically. We showed that model (2.3) undergoes flip bifurcation at $\alpha_h$ in Case (i) and Neimark-Sacker bifurcation at $\alpha_h$ in Case (ii), for a different set of parameter values. We notice that model (2.3) has strange dynamics for smaller values of $\alpha$; see Figures 2(a) and 5(a). The strange behavior of the model (2.3) for both cases is shown clearly in the one-parameter bifurcation diagrams in Figures 3(a) and 6(a). We notice that the proposed model reduces complex dynamical behavior, and we say that the sensitivity of the population dynamics to initial conditions has been reduced in a sense by the Allee effect. Further, we verified this sensitivity to the initial conditions by plotting the largest Lyapunov exponents in Figures 3(b) and 6(b). In this study, we successfully implemented state and OGY feedback and hybrid control methods by constructing the corresponding controlled models (5.1), (5.2) and (5.3). The lines of stable margins for the controlled models using state feedback control and the OGY feedback control method have been derived and plotted in Figures 8(a) and 9(a). The choice of feedback value inside the stability triangle helps to shift the desired chaotic orbit into the area of stable dynamics, see Figures 8(b), (c) and 9(b), (c). In the case of the hybrid control method, for the control parameter $0 < \epsilon < 0.773477$, the model (5.3) is stable, and the one-parameter bifurcation diagram is shown in Figure 10. A different possibility for showing the existence of chaos in real populations has been provided in the theory of chaos control ^[26]. The need for this control supports the adoption of the OGY technique for the three fish food chain, i.e., stabilization for the desired periodic orbit. This control condition goes against empirical attempts for lakes and rivers ^[23].

This study concludes that increasing the Allee parameter, which leads to converting the strange attractor into a proper pattern, i.e., period-16, period-8, period-4, period-2, period-13, period-9 and a closed invariant circle, means that the population oscillates periodically and finally reaches stable dynamics (long-term coexistence of both species) in the considered model. Also, this periodic window and chaotic orbits cannot occur in the counter continuous-time model. This stabilization takes place via both the flip and the Neimark-Sacker bifurcations. In one word, the Allee effect plays a positive role in stabilizing the proposed model, i.e., keeping the population persistent for a long time. However, there will be more complexity in predicting the future population size under chaotic dynamics because a small perturbation in the initial population size can cause massive divergence in population size. For smaller values of the Allee parameter, the presence of bifurcation and chaos in the proposed model can cause both species to have a higher risk of extinction due to unpredictability.

However, it is important to study the discrete-time predator-prey model with the Allee effect in both prey and predator populations. It is interesting to note that one may study the model's relationship to other ecological phenomena like prey refuge, cannibalism and fear impact by appropriately adjusting the strength of the Allee effect simultaneously. Also, considering the other functional responses such as ratio-dependent, Crowley-Martin and Beddington-DeAngelis functional responses, may lead to complex models and cause strange behaviors; we will leave this for future research.

Use of AI tools declaration

The authors declare that they have not used artificial intelligence tools in the creation of this article.

Acknowledgments

This work was partially supported by the Center for Nonlinear Systems, Chennai Institute of Technology, India, with funding number CIT/CNS/2023/RP/016.

Conflict of interest

Prof. Dimplekumar Chalishajar is the Guest Editor of special issue "Recent advances in differential and partial differential equations and its applications'' for AIMS Mathematics. Prof. Dimplekumar Chalishajar was not involved in the editorial review or the decision to publish this article.

All authors declare no conflict of interest this paper.