Stability analysis and chaos control in a discrete predator-prey system with Allee effect, fear effect, and refuge

1.   Introduction

In the natural world, the mutual relationship between predators and prey represents one of the most fundamental and critically important types of ecological interactions within population dynamics. Predators exert direct influence on prey population sizes and simultaneously affect their own population dynamics, resulting in complex interplays between them. A fundamental continuous population model was devised in 1926 by American biologist Lotka and Italian mathematician Volterra to explain the relationship between predators and prey ^[1]. Subsequent to this foundational work, numerous experts conducted in-depth investigations into the predator-prey model ^[2,3,4]. Nevertheless, due to the model's limited ability to address diverse real-world scenarios and complexities, scholars have developed multifaceted modifications to provide a more realistic interpretation and enhance the understanding of ecosystems. To present a more precise depiction of ecological interactions, several ecological concepts have been incorporated into the predator-prey system. These include functional responses, shelters, harvesting, fear, and the Allee effect ^[5,6,7].

In the natural world, numerous species of animals, plants, and insects frequently exhibit reduced individual fitness at lower population densities. The ecologist W.C. Allee succinctly summarized this phenomenon as the Allee effect ^[8]. The Allee effect can be attributed to various environmental factors, including difficulties in mating at low densities, genetic inbreeding, reduced anti-predator defense capabilities, depression, and social hindrances at low densities, among others ^[9,10,11]. The Allee effect is dichotomized into two distinct categories based on the magnitude of density dependence at low population densities: The strong Allee effect and the weak Allee effect. In the case of the strong Allee effect, the population growth rate becomes negative when the population falls below the Allee effect threshold. Conversely, a weak Allee effect is characterized by a per capita population growth rate at low density that is lower than that at high density, though it remains positive. For endangered species, Allee effects are more likely to occur, underscoring the significance of Allee effects in the management of endangered species conservation, population development, utilization, and species introduction. Recent years have witnessed attention being increasingly dedicated to comprehending the impact of the Allee effect on biological populations ^[12,13,14], revealing its significant influence on the dynamics of ecological systems.

In the investigation of predator-prey systems, the predominant focus in most studies has revolved around the direct predation exerted by predators on prey populations, given the ease of directly observing such behavior within ecosystems. In recent years, scholars have gradually become aware that prey, upon sensing the presence of predators, exhibit fear responses. This fear can have diverse impacts on various aspects of prey populations, including changes in habitat use, foraging behavior, and vigilance, as well as different psychological changes ^[15,16,17]. On certain occasions, the effects of this fear may be even more significant than those of direct predation ^[18,19,20]. To date, numerous scholars have conducted research on predator-prey models that incorporate the effects of fear ^[21,22,23].

The efficiency of predators in prey capture is intricately tied to the strategies employed by the prey to evade predation. In situations of danger, prey organisms typically seek refuge in areas with lower predation risks, ecologically referred to as prey refuges ^[24,25,26]. The successful evasion of predation is achieved by certain prey populations through the exploration of refuges, thereby effectively mitigating the pressure of being captured by predators. In 1946, Crombic introduced the concept of refuge into the prey-predator model through experimental investigations ^[27]. Both theoretical analysis and experimental studies suggest that this type of system is more realistic and holds greater research value ^[28,29,30].

When examining populations characterized by non-overlapping generations or low numbers, mathematical expressions often present as discrete models. Despite certain variables exhibiting continuous changes in practical applications, the recording of quantity changes often occurs within specific time intervals during data processing. Therefore, utilizing discrete systems to analyze the dynamic behavior of biological populations has immense practical significance. For example, Hadeler and Gerstmann ^[31] showed that for the Rosenzweig-MacArthur model, the discrete-time version of the model has similar equilibrium points as the continuous-time version, but the discrete-time version of the model exhibits more complex phenomena such as period doubling and chaos, which are phenomena that are not present in the continuous-time model. In turn, we can conclude that discrete systems are better suited to describe the richer patterns and chaotic behavior of nonlinear dynamics. Numerous mathematical methods exist for the discretization of continuous systems, including the forward Euler method, the use of piecewise constant arguments, and the exponential Euler difference scheme ^{[32,33,34,35,36]}. Among these, the forward Euler method is widely employed due to its ease of comprehension and computational simplicity. Inspired by the aforementioned literature, we chose to adopt the forward Euler method to discretize continuous systems.

Our focus on understanding system stability and controlling chaos has been driven by the desire to clarify the intricate dynamics between predators and prey in ecosystems, as well as to probe the complexity and stability of community dynamics. Stability analysis offers insights into system resilience under varying environmental conditions, shedding light on how environmental factors influence community structure and evolution. Conversely, the emergence of chaotic behavior can disrupt system predictability and stability, posing risks to ecosystem equilibrium and biodiversity. By investigating stability and chaos control in predator-prey models, we contribute to the broader understanding of ecosystem dynamics, offering insights and methodologies that are crucial for sustainable ecosystem management and conservation.

In 2014, Rana et al. conducted a study on a discrete model of predator-prey systems that incorporate an Allee effect and prey refuge. The findings revealed that moderate levels of prey refuge can stabilize the dynamics of the system through the Allee effect. Additionally, it was observed that the population remained stable at moderate refuge parameter levels, while chaotic oscillations of prey occurred at relatively low and high refuge levels ^[37]. In 2020, Ma et al. investigated a discrete model of predator-prey dynamics involving fear and refuge. The study identified flip branching at immobile points where only the predator is present, as well as Neimark-Sacker branching at positive immobile points where both prey and predator coexist ^[38]. In 2023, Hong and Zhang explored a discrete predator-prey system that incorporates an Allee effect and a fear factor. The investigation demonstrated that the system undergoes flip bifurcation and Neimark-Sacker bifurcation at the sole positive immobile point ^[39]. To the best of our knowledge, there is a dearth of studies focusing on the stability, bifurcation, and chaos control analysis of discrete-time predator-prey models that incorporate a strong Allee effect, the fear effect, and prey refuge. This article serves to contribute to this area of study and address existing gaps in the literature. This study brings to light several significant contributions, notably, those summarize below:

● Utilizing the forward Euler method, we have derived a discrete-time model that aptly represents interactions among non-overlapping populations across generations.

● Stability analysis of the discrete model was conducted from the perspective of the fixed points inherent in the system.

● Investigating internal fixed points, various types of bifurcations were explored by introducing $h$ as the bifurcation parameter, leveraging the center manifold theorem and bifurcation theory. Notably, the proposed discrete system exhibited flip and Neimark-Sacker bifurcations.

● In addressing system-induced chaos, chaos control was implemented in the discrete model through state feedback control, OGY feedback control, and hybrid control methods, effectively managing the chaotic behavior of the system.

The rest of the paper is organized as follows. In Section 2, we formulate a discrete biodynamic system by using the forward Euler method. Subsequently, in Section 3, we analyze the existence and stability of fixed points within the discrete system. Section 4 focuses on demonstrating the possibility of flip bifurcation and Neimark-Sacker bifurcation in the discrete system, considering appropriate parameter values. In Section 5, chaos is effectively controlled through the application of state feedback control, OGY feedback control, and hybrid control methods. The practicality of the main theoretical results is validated through numerical simulations in Section 6. Finally, the paper concludes with a brief summary.

2.   Model

First, we describe the continuous-time prey-predator model with the Allee effect, a fear effect, and a prey refuge ^[40], as follows:

where all coefficients are positive constants; $u$ and $v$ are the prey and predator population densities, respectively; $r$ represents the intrinsic growth rate of the prey; $a$ indicates the predator's rate of predation on prey; $k$ is the maximum prey carrying capacity of the environment; $\alpha$ is the transformation efficiency of a predator as a result of preying on its prey; ${{m}_{0}}$ is the natural mortality rate of predators; $0 < {{\theta }_{0}} < k$ indicates a strong Allee effect; $F\left({{f}_{0}}, v \right) = \frac{1}{1+{{f}_{0}}v}$ is the fear function; $0 < \eta < 1$ is the prey refuge constant.

For simplicity, we initially rescale system (2.1) by introducing the following parameters: $N = \frac{u}{k}, P = \frac{av}{rk}, \tau = a\alpha kt, \varepsilon = \frac{a\alpha }{r}, \theta = \frac{{{\theta }_{0}}}{k}, f = \frac{rk{{f}_{0}}}{a}, m = \frac{{{m}_{0}}}{a\alpha k}$. Rewriting $\tau$ as $t$, the simplified system (2.1) becomes as follows:

where $0 < \theta < 1$ and $0 < m < 1$.

Implementing the forward Euler scheme to system (2.2) yields the discrete model under scrutiny in this investigation:

where $h$ is the step size. The system of difference equation (2.3) can be written in mapping form as follows:

3.   The existence and stability of fixed points

In this section, the presence of fixed points for the system (2.4) is investigated. Subsequently, an investigation into the stability of these fixed points is conducted, employing either the characteristic polynomial or the eigenvalues of the Jacobian matrix assessed at these fixed points.

The subsequent lemmas aid in examination of the stability of fixed points.

Lemma 3.1. ^[41] Let $F(\lambda) = {{\lambda }^{2}}+p\lambda +q$ be the characteristic equation of the Jacobian matrix at $(x, y)$ and ${{\lambda }_{1}}$, ${{\lambda }_{2}}$ be solutions of $F(\lambda) = 0$; then, $(x, y)$ is a

(1) sink iff $|{{\lambda }_{1}}| < 1$ and $|{{\lambda }_{2}}| < 1$,

(2) saddle point iff $|{{\lambda }_{1}}| < 1$ and $\left| {{\lambda }_{2}} \right| > 1$ (or $\left| {{\lambda }_{1}} \right| > 1$ and $\left| {{\lambda }_{2}} \right| < 1$),

(3) source iff $|{{\lambda }_{1}}| > 1$ and $\left| {{\lambda }_{2}} \right| > 1$,

(4) non-hyperbolic point iff $|{{\lambda }_{1}}| = 1$ or $|{{\lambda }_{2}}| = 1$.

Lemma 3.2. ^[41] Let $F(\lambda) = {{\lambda }^{2}}+p\lambda +q$, where $p$, $q$ are constants. Suppose that $F(1) > 0$ and ${{\lambda }_{1}}$, ${{\lambda }_{2}}$ are roots of $F(\lambda) = 0$; then,

(1) $|{{\lambda }_{1}}| < 1$ and $|{{\lambda }_{2}}| < 1$ iff $F(-1) > 0$ and $F(0) < 1$,

(2) $|{{\lambda }_{1}}| < 1$ and $\left| {{\lambda }_{2}} \right| > 1$ (or $\left| {{\lambda }_{2}} \right| < 1$ and $\left| {{\lambda }_{1}} \right| > 1$) iff $F(-1) < 0$,

(3) $|{{\lambda }_{1}}| > 1$ and $\left| {{\lambda }_{2}} \right| > 1$ iff $F(-1) > 0$ and $F(0) > 1$,

(4) ${{\lambda }_{1}} = -1$ and $\left| {{\lambda }_{2}} \right|\ne 1$ iff $F(-1) = 0$ and $p\ne 0, 2$,

(5) ${{\lambda }_{1}}$ and ${{\lambda }_{2}}$ are complex and $\left| {{\lambda }_{1}} \right| = \left| {{\lambda }_{2}} \right| = 1$ iff ${{p}^{2}}-4q < 0$ and $F(0) = 1$.

3.1.   Existence of fixed points

In this subsection, we will compute the fixed points of system (2.4) and establish the conditions for the asymptotic stability of these fixed points. To determine the fixed points of system (2.4), we need to create the following system of equations:

Solving the algebraic equation, we get that system (2.4) has the following fixed points:

Theorem 3.1. For system (2.4), the following statements hold true:

(1) A trivial fixed point ${{E}_{0}} = (0, 0)$,

(2) there are two boundary fixed points ${{E}_{1}} = (1, 0)$ and ${{E}_{2}} = (\theta, 0)$,

(3) there is a unique positive fixed point ${{E}^{*}} = ({{N}^{*}}, {{P}^{*}})$, where ${{N}^{*}} = \frac{m}{1-\eta }$, ${{P}^{*}} = \frac{-1+\sqrt{1+4fH}}{2f}$ and ${H^{*}} = \frac{(1-\eta -m)(m-\theta (1-\eta))}{{{\left(1-\eta \right)}^{3}}}, 1-\frac{m}{\theta } < \eta < 1-m$.

3.2.   Stability of fixed points

Theorem 3.2. The trivial fixed point ${{E}_{0}} = (0, 0)$ is a

(1) sink if $0 < h < \min \{\frac{2\varepsilon }{\theta }, \frac{2}{m}\}$,

(2) source if $\max \{\frac{2\varepsilon }{\theta }, \frac{2}{m}\} < h$,

(3) saddle if $\frac{2\varepsilon }{\theta } < h < \frac{2}{m}$ or $\frac{2}{m} < h < \frac{2\varepsilon }{\theta }$,

(4) non-hyperbolic point if $h = \frac{2}{m}$ or $h = \frac{2\varepsilon }{\theta }$.

Proof. The Jacobian matrix computed at ${E}_{0}$ is given by

The eigenvalues of $J({{E}_{0}})$ are ${{\lambda }_{1}} = 1-\frac{h\theta }{\varepsilon }$ and ${{\lambda }_{2}} = 1-hm$. By calculating this, we can get the following:

□

Remark 1. In the context of system (2.1) ^[40], the fixed point $(0, 0)$ consistently maintains stability. Within the framework of system (2.4), it manifests complex dynamics. From this it can be seen that discrete systems have more complex dynamics than continuous systems.

Theorem 3.3. The boundary fixed point ${{E}_{1}} = (1, 0)$ is a

(1) sink if $\left\{ \begin{aligned} & m > 1-\eta, \\ & 0 < h < \min \{\frac{2\varepsilon }{1-\theta }, \frac{2}{m+\eta -1}\}, \\ \end{aligned} \right.$

(2) source if $\left\{ \begin{aligned} & m > 1-\eta, \\ & h > \max \{\frac{2\varepsilon }{1-\theta }, \frac{2}{m+\eta -1}\}, \\ \end{aligned} \right.$ or $\left\{ \begin{aligned} & m < 1-\eta, \\ & h > \frac{2\varepsilon }{1-\theta }, \\ \end{aligned} \right.$

(3) saddle if any one of the conditions listed below is met:

(a) $\left\{ \begin{aligned} & m > 1-\eta, \\ & \frac{2\varepsilon }{1-\theta } < h < \frac{2}{m+\eta -1}, \\ \end{aligned} \right.$

(b) $\left\{ \begin{aligned} & m > 1-\eta, \\ & \frac{2}{m+\eta -1} < h < \frac{2\varepsilon }{1-\theta }, \\ \end{aligned} \right.$

(c) $\left\{ \begin{aligned} & m < 1-\eta, \\ & 0 < h < \frac{2\varepsilon }{1-\theta }, \\ \end{aligned} \right.$

(4) non-hyperbolic point if any one of the conditions listed below is met:

(a) $h = \frac{2\varepsilon }{1-\theta }$,

(b) $\left\{ \begin{aligned} & m > 1-\eta, \\ & h = \frac{2}{m+\eta -1}, \\ \end{aligned} \right.$

(c) $m = 1-\eta$.

Proof. The Jacobian matrix computed at ${{E}_{1}} = (1, 0)$ is given by

The eigenvalues are ${{\lambda }_{1}} = 1+\frac{h\left(\theta -1 \right)}{\varepsilon }$ and ${{\lambda }_{2}} = 1+h\left(1-\eta -m \right)$. By calculating this, we can get the following:

and

□

Theorem 3.4. The boundary fixed point ${{E}_{2}} = (\theta, 0)$ is a

(1) source if $\left\{ \begin{aligned} & m > \theta (1-\eta), \\ & h > \frac{2}{m-\theta (1-\eta)}, \\ \end{aligned} \right.$ or $m < \theta (1-\eta)$,

(2) saddle if $\left\{ \begin{aligned} & m > \theta (1-\eta), \\ & 0 < h < \frac{2}{m-\theta (1-\eta)}, \\ \end{aligned} \right.$

(3) non-hyperbolic point if $m = \theta (1-\eta)$ or $\left\{ \begin{aligned} & m > \theta (1-\eta), \\ & h = \frac{2}{m-\theta (1-\eta)}. \\ \end{aligned} \right.$

Proof. The Jacobian matrix computed at ${{E}_{2}} = (\theta, 0)$ is given by

The eigenvalues are ${{\lambda }_{1}} = 1+\frac{h\theta }{\varepsilon }\left(1-\theta \right)$ and ${{\lambda }_{2}} = 1+h\left(\theta \left(1-\eta \right)-m \right)$. By calculating this, we can get the following: $\left| 1+\frac{h\theta }{\varepsilon }\left(1-\theta \right) \right| > 1$ is constant, and

□

Theorem 3.5. The following conditions hold for the unique positive fixed point ${{E}^{*}}$:

(1) ${{E}^{*}} = ({{N}^{*}}, {{P}^{*}})$ is a sink if any one of the conditions listed below is met:

(a) $-2\sqrt{S}\le A < 0$ and $0 < h < -\frac{A}{S}$,

(b) $A < -2\sqrt{S}$ and $0 < h < \frac{-A-\sqrt{{{A}^{2}}-4S}}{S}$,

(2) ${{E}^{*}} = ({{N}^{*}}, {{P}^{*}})$ is a source if any one of the conditions listed below is met:

(a) $-2\sqrt{S}\le A < 0$ and $h > -\frac{A}{S}$,

(b) $A < -2\sqrt{S}$ and $h > \frac{-A+\sqrt{{{A}^{2}}-4S}}{S}$,

(c) $A\ge 0$,

(3) ${{E}^{*}} = ({{N}^{*}}, {{P}^{*}})$ is a saddle if $A < -2\sqrt{S}$ and $\frac{-A-\sqrt{{{A}^{2}}-4S}}{S} < h < \frac{-A+\sqrt{{{A}^{2}}-4S}}{S}$,

(4) ${{E}^{*}} = ({{N}^{*}}, {{P}^{*}})$ is non-hyperbolic if any one of the conditions listed below is met:

(a) $-2\sqrt{S} < A < 0$ and $h = -\frac{A}{S}$,

(b) $A < -2\sqrt{S}$, $h = \frac{-A\pm \sqrt{{{A}^{2}}-4S}}{S}$ and $h\ne -\frac{2}{A}, h\ne -\frac{4}{A}$,

where $A = \frac{m\left(\left(\theta +1 \right)\left(1-\eta \right)-2m \right)}{\varepsilon {{\left(1-\eta \right)}^{2}}\left(f{{P}^{*}}+1 \right)}$, $S = \frac{m}{\varepsilon }\left(\frac{(1-\eta -m)\left(m-\theta \left(1-\eta \right) \right)f}{{{\left(1-\eta \right)}^{2}}{{\left(f{{P}^{*}}+1 \right)}^{2}}}+1-\eta \right){{P}^{*}}$.

Proof. The Jacobian matrix computed at the point ${{E}^{*}} = ({{N}^{*}}, {{P}^{*}})$ is as follows:

The characteristic polynomial of $J({{E}^{*}})$ is given by

where

Then, we get

It is clear that $F\left(1 \right) = S{{h}^{2}} > 0$.

Thus, we can conclude with the result. □

4.   Bifurcation analysis

4.1.   Flip bifurcation

Here, we choose $h$ as the bifurcation parameter. According to the analysis in Section 3, it is determined that the model can undergo a flip bifurcation when $h$ varies within a small vicinity of ${{L}_{1}}$ or ${{L}_{2}}$, where

We exclusively examine the flip bifurcation phenomenon when the parameter $h$ undergoes variations within a small neighborhood of ${{L}_{1}}:h = \frac{-A-\sqrt{{{A}^{2}}-4S}}{S}$. Analogous reasoning can be extended to the case of ${{L}_{2}}:h = \frac{-A+\sqrt{{{A}^{2}}-4S}}{S}$. Now, we consider a perturbed system derived from system (2.4):

where $h$ represents the bifurcation parameter and $\left| {{h}^{*}} \right|\ll 1$ serves as a small perturbation parameter.

Let $x = N-{{N}^{*}}, y = P-{{P}^{*}}$. Then, we get the following:

where

We construct a nonsingular matrix:

and we use the following transformation: $\left(\begin{aligned} & x \\ & y \\ \end{aligned} \right) = T\left(\begin{aligned} & X \\ & Y \\ \end{aligned} \right)$. Then, system (4.3) will become as follows:

where

In accordance with the center manifold theorem, there is the existence of a center manifold denoted as ${{W}^{c}}\left(0, 0 \right)$:

and

Through calculation, we obtain

The system (2.4) restricted to the center manifold, is given by

where

Now, regarding the occurrence of flip bifurcation, it is crucial that the two quantities, ${\alpha }_{1}$ and ${\alpha }_{2}$, are both non-zero, where

As a result of the aforementioned analysis, we summarize the following lemma on the existence and direction of the flip bifurcation.

Lemma 4.1. If ${{\alpha }_{1}}\ne0$, ${{\alpha }_{2}}\ne 0$, then system (2.4) undergoes a flip bifurcation at the positive equilibrium point ${{E}^{*}} = ({{N}^{*}}, {{P}^{*}})$ when $h$ in a small neighborhood of ${{L}_{1}}$. And when ${{\alpha }_{2}} > 0$ (respectively, ${{\alpha }_{2}} < 0$), system (2.4) bifurcates from the fixed point ${{E}^{*}} = ({{N}^{*}}, {{P}^{*}})$ to a 2-periodic stable orbit (respectively, unstable).

4.2.   Neimark-Sacker bifurcation

Neimark-Sacker bifurcation manifests when the eigenvalues of the characteristic equation form complex conjugate pairs with a unit modulus. Based on the preceding discussion, we know that the model can undergo a Neimark-Sacker bifurcation when $h$ varies within a small neighborhood of $L_{3}$, where

Now, we examine a perturbation system derived from system (2.4):

where $h$ represents the bifurcation parameter, and $\left| {\bar{h}} \right|\ll 1$ serves as a small perturbation parameter. Let $x = N-{{N}^{*}}$ and $y = P-{{P}^{*}}$:

where ${{a}_{ij}}\left(i = 1, 2, j = 1, 2, \cdots 9 \right)$ are given in (4.4) by changing $h$ to $h+\bar{h}$. The characteristic equation corresponding to the linear part of the system (4.13) at the fixed point (0, 0) is given by

where $p\left({\bar{h}} \right) = -2-A\left(h+\bar{h} \right)$, $q\left({\bar{h}} \right) = 1+A\left(h+\bar{h} \right)+S{{\left(h+\bar{h} \right)}^{2}}$.

The roots of (4.14) are complex and given by

When $\bar{h} = 0$ and $h\in l$, by computation, we obtain $q\left(0 \right) = 1, {{\lambda }_{1}}\left(0 \right)$, ${{\lambda }_{2}}\left(0 \right) = 1+\frac{Ah}{2}\pm \frac{h\sqrt{4S-{{A}^{2}}}}{2}i$, and $\left| {{\lambda }_{1}}\left(0 \right) \right| = \left| {{\lambda }_{2}}\left(0 \right) \right| = 1$.

It is clear from the above discussion that

In addition, it is required that when $\bar{h} = 0$, ${{\lambda }_{1}}^{n}, {{\lambda }_{2}}^{n}\ne 1\left(n = 1, 2, 3.4 \right)$, which is equivalent to $p\left(0 \right)\ne 0, 1, 2, -2$. Because $h\in {{L}_{3}}$, we get that $p\left(0 \right)\ne 2, -2$. Additionally, we necessitate that $p(0) \ne 0$ and $1$, which gives

Let ${{\lambda }_{1}}\left(0 \right), {{\lambda }_{2}}\left(0 \right) = \rho \pm \omega i$. To bring the linear component of (4.13) into its canonical form at $\bar{h} = 0$, we apply the following transformation:

So $x = \omega u+\left({{a}_{11}}-\rho \right)v$ and $y = {{a}_{12}}v$. Under the transformation (4.17), the system (4.13) will become as follows:

where

Then, we get

Next, a non zero real number is defined as follows:

where

Based on the aforementioned calculations, we establish the lemma regarding the presence and direction of Neimark-Sacker bifurcation.

Lemma 4.2. If the condition (4.16) holds and $\tilde{l}\ne 0$, then system (2.4) undergoes a Neimark-Sacker bifurcation at the fixed point ${{E}^{*}} = ({{N}^{*}}, {{P}^{*}})$ when the parameter $h$ varies in a small neighborhood of ${{L}_{3}}$. Moreover, if $\tilde{l} < 0$ (resp., $\tilde{l} > 0$), then an attracting (resp., repelling) invariant closed curve bifurcates from the fixed point for $\bar{h} > 0$ (resp., $\bar{h} < 0$).

5.   Chaos control

Chaos is often a detrimental phenomenon for biological systems; they often disrupt the ecological balance of biological population systems and directly affect one's long-term projections of population growth. By implementing appropriate control policies, not only can the size of ecological populations be protected, human beings can also lay a good foundation for the sustainable exploitation of ecological resources. Therefore, three control methods are used in this section to provide an effective control of the chaos generated by system (2.4).

5.1.   State feedback control

In this section, the state feedback control method ^[42] will be employed to regulate the chaos exhibited by system (2.4). To achieve this, a feedback control term is introduced into system (2.4):

where

denotes the feedback controlling force, while ${r}_{1}$ and ${r}_{2}$ represent the feedback gains. The point $({{N}^{*}}, {{P}^{*}})$ is the unique positive equilibrium point of system (2.4).

The Jacobian matrix at the equilibrium point $({{N}^{*}}, {{P}^{*}})$ for system (5.1) is given by

where

Therefore, the characteristic equation associated with $J({{N}^{*}}, {{P}^{*}})$ is given by

Let ${{\lambda }_{1}}$ and ${{\lambda }_{2}}$ represent the eigenvalues of the characteristic Eq (5.5). Then,

Then, we can determine the lines of marginal stability:

Theorem 5.1. If ${{r}_{1}}$, ${{r}_{1}}$ lie in a triangular region bounded by lines ${{l}_{1}}$, ${{l}_{2}}$, and ${{l}_{3}}$, it can be concluded that the system (5.1) is stable.

5.2.   OGY control method

In this subsection, we investigate a technique for chaos control that relies on the OGY method, as outlined in ^[43]. By taking $m$ as the control parameter, we rewrite the system (2.4) as follows:

Furthermore, let us assume that $m\in \left({{m}_{0}}-{{\delta }_{0}}, {{m}_{0}}+{{\delta }_{0}} \right)$ with ${{\delta }_{0}} > 0$ and ${{m}_{0}}$ denotes the nominal value of $m$. Then, model (5.9) can be approximated near $({{N}^{*}}, {{P}^{*}})$ by applying the following linear map:

where

The controllability of model (5.9) is evident when considering the matrix

and ensuring that it has rank 2.

Next, we assume that $\left(m-{{m}_{0}} \right) = -K\left(\begin{matrix} {{N}_{n}}-{{N}^{*}} \\ {{P}_{n}}-{{P}^{*}} \\ \end{matrix} \right)$, where $K = \left[\begin{matrix} {{k}_{1}} & {{k}_{2}} \\ \end{matrix} \right]$. Then, model (5.9) can be written as

Then, the controller model is given by

Furthermore, the equilibrium point $({{N}^{*}}, {{P}^{*}})$ achieves local asymptotic stability if and only if both eigenvalues of the matrix $J-MK$, denoted as ${{\lambda }_{1}}$ and ${{\lambda }_{2}}$, reside within an open unit disk. The matrix $J-MK$ can be expressed as follows:

where

Then, the characteristic polynomial is given by

Subsequently, we can derive the lines of marginal stability, denoted as follows:

Theorem 5.2. If ${{k}_{1}}$, ${{k}_{1}}$ lie in a triangular region bounded by lines ${{m}_{1}}$, ${{m}_{2}}$, and ${{m}_{3}}$, it can be concluded that the system (5.13) is stable.

5.3.   Hybrid control

In this subsection, we employ the hybrid control technique ^[44] to control chaos. The controlled system discussed herein corresponds to system (2.4), as depicted below:

where $0 < \zeta < 1$. This is a control strategy that combines feedback control and parameter adjustment. The Jacobian matrix at the equilibrium point $({{N}^{*}}, {{P}^{*}})$ for the system given by Eq (5.19) is as follows:

where

According to the Jury condition, the equilibrium point of the system is stable if and only if the following conditions are satisfied:

Theorem 5.3. If $\left| {{\sigma }_{11}}+{{\sigma }_{22}} \right| < 1+{{\sigma }_{22}}{{\sigma }_{11}}-{{\sigma }_{21}}{{\sigma }_{12}} < 2$ can hold, it can be concluded that the system (5.19) is stable.

6.   Numerical simulations

The objective of this section is to illustrate the feasibility of the main theoretical results and to depict the intricate dynamical behaviors of system (2.4). To achieve this, we present bifurcation diagrams and phase portraits through numerical simulations. Exploring the impact of various biological effects on discrete predator-prey systems is crucial. In this section, we will respectively investigate the effects of the fear effect, refuges, and the Allee effect on the system.

6.1.   Flip bifurcation

We chose the parameter values as follows:

The flip bifurcation value is computed as $h = \frac{-A-\sqrt{{{A}^{2}}-4S}}{S} = 0.8085$, and the interior fixed point is computed as ${{E}^{*}} = \left(0.9859, 0.00956 \right)$. The eigenvalues of the Jacobian matrix are ${{\lambda }_{1}} = -1$ and ${{\lambda }_{2}} = 0.991554931907570$, as well as ${{\alpha }_{1}} = -1.235833646$ and ${{\alpha }_{2}} = 11.9477919$, affirming that system (2.4) undergoes flip bifurcation at ${{E}^{*}}$ as the bifurcation parameter $h$ crosses $h = 0.8085$. Referring to Figure 1, we see that the bifurcation values obtained in Figure 1(a) align with the conclusions from Lemma 4.1, and it is evident that the prey-eating population experiences a flip bifurcation, eventually leading to chaos. Furthermore, Figure 1(b) shows the maximum Lyapunov exponents to confirm the chaotic behavior of system (2.4). When we compare Figure 1(b) to Figure 1(a), we can observe that the prey population enters a chaotic state when the maximum Lyapunov exponent is greater than 0.

6.2.   Neimark-Sacker bifurcation

We chose the parameter values as follows:

The Neimark-Sacker bifurcation value is computed as $h = -\frac{A}{S} = 0.6226263491$ and the interior fixed point of system (2.4) is identified as ${{E}^{*}} = \left(0.6185, 0.1737 \right)$. The eigenvalues of the Jacobian matrix are given by ${{\lambda }_{1, 2}} = 0.9775042841\pm 0.2109155626i$ with $\left| {{\lambda }_{1}} \right| = \left| {{\lambda }_{2}} \right| = 1$, as well as $d = 0.03613036280 > 0$, confirming the occurrence of Neimark-Sacker bifurcation at ${{E}^{*}}$. Additionally, the value of the first Lyapunov exponent is computed as $\tilde{l} = -0.00080098583420$, verifying the validity of Lemma 4.2. The system's two-dimensional bifurcation diagram is shown in Figure 2(a). As $h$ increases, we can see that the system's fixed point shifts from a stable to a chaotic state. The maximum Lyapunov exponent map corresponding to Figure 2(a) is given in Figure 2(b); we observe that, when the maximum Lyapunov exponent exceeds 0, the prey population enters a chaotic state. Some classical phase diagrams are given in Figure 3, which clearly show how smooth invariant circles bifurcate from stable fixed points ${{E}^{*}} = \left(0.6185, 0.1737 \right)$ and eventually develop into chaotic sets.

Remark 2. In the context of the system described in ^[39], it was found that the system remained non-chaotic as parameters changed during the Neimark-Sacker bifurcation. However, our system in this study demonstrated chaos in the system as parameter $h$ changed Neimark-Sacker bifurcation. Hence, we conclude that the addition of prey refuge increases the complexity of the biological system.

6.3.   Chaos control

We chose the parameter values as follows:

By evaluating Figure 3, we can get that the variables $N$ and $P$ in the system (2.4) are in a chaotic state when the parameter is (6.3); the corresponding time series plot is as follows; see Figure 4.

6.3.1.   State feedback control

We consider the same parameter values as in (6.3). In Figure 5, the triangular region, as determined by Theorem 5.1, confines the parameters ${{r}_{1}}$ and ${{r}_{2}}$. Within this region, the chaotic behavior generated by the system (2.4) is effectively controlled, leading to asymptotic convergence near the fixed point ${{E}^{*}} = \left(0.6185, 0.1737 \right)$.

For the case that the feedback gains are ${{r}_{1}} = 1.27$ and ${{r}_{2}} = 1.1$ and the controller is initiated at the 500th iteration of the system, Figure 6 illustrates the control effect, wherein a chaotic trajectory is successfully stabilized at the fixed point ${{E}^{*}} = \left(0.6185, 0.1737 \right)$. Consequently, it can be inferred that the utilization of the feedback control approach proves effective in the mitigation of bifurcation and chaos.

6.3.2.   OGY control method

We consider the same parameter values as in (6.3). In Figure 7, the triangular region, as determined by Theorem 5.2, confines the parameters ${{k}_{1}}$ and ${{k}_{2}}$. Within this region, the chaotic behavior generated by the system (2.4) is effectively controlled, leading to asymptotic convergence near the fixed point ${{E}^{*}} = \left(0.6185, 0.1737 \right)$.

The values ${{k}_{1}} = 0.2$ and ${{k}_{2}} = -3$ were taken in the stabilization region and the controller was initiated at the 500th iteration of the system; Figure 8 illustrates the control effect, wherein a chaotic trajectory is successfully stabilized at the fixed point ${{E}^{*}} = \left(0.6185, 0.1737 \right)$. Consequently, it can be inferred that the utilization of the OGY control method approach proves effective in the mitigation of bifurcation and chaos.

6.3.3.   Hybrid control

For the controlled system (5.19), when we consider the same parameter values as in (6.3), according to Theorem 5.3, the stability range of the hybrid control system is $0 < \xi < 0.831662531$. Taking $\xi = 0.4$, Figure 9 documents the bifurcation of the controlled system; by comparison with Figure 2(a), it can be seen that the Neimark-Sacker bifurcation of the system is delayed. By employing small values of the control parameter $\xi$, the Neimark-Sacker bifurcation can be postponed over a broader range of $h$. Consequently, it can be inferred that the utilization of the hybrid control approach proves effective in the mitigation of bifurcation and chaos.

6.4.   The impact of the fear effect on the system

We chose the parameter values as follows:

Neimark-Sacker bifurcation value is calculated as $h = -\frac{A}{S} = 0.6775878938$, and the internal immobility point is calculated as ${{E}^{*}} = \left(0.6186, 0.1928\right)$. See Figure 10. Comparison with Figure 2(a) shows that the Neimark-Sacker bifurcation of the system was delayed from $h = 0.6226263491$ to $h = 0.6775878938$ when the fear effect constant was varied from $1$ to $0.3$. This means that the change of biological populations in the system will not be too drastic. This helps to maintain biodiversity in the ecosystem so that various biological populations can coexist in a relatively stable manner.

6.5.   The impact of the refuge on the system

We chose the parameter values as follows:

The Neimark-Sacker bifurcation value is calculated as $h = -\frac{A}{S} = 0.4954418794$, and the internal immobility point is calculated as ${{E}^{*}} = \left(0.6061, 0.1718\right)$. See Figure 11. Comparison with Figure 2(a) shows that the Neimark-Sacker bifurcation of the system was advanced from $h = 0.6226263491$ to $h = 0.4954418794$ when the refuge constant was varied from $0.03$ to $0.01$. This means that ecosystems are more susceptible to external disturbances, thereby increasing the risk of abrupt changes in the system. This is an important warning sign for environmental conservation and resource managers that ecosystem change and stability need to be more closely monitored.

6.6.   The impact of the Allee effect on the system

We chose the parameter values as follows:

The Neimark-Sacker bifurcation value is calculated as $h = -\frac{A}{S} = 0.7678598057$, and the internal immobility point is calculated as ${{E}^{*}} = \left(0.6186, 0.1882\right)$. See Figure 12. In terms of comparison with Figure 2(a), Figure 12 shows that the Neimark-Sacker bifurcation of the system was delayed from $h = 0.6226263491$ to $h = 0.7678598057$ when the Allee effect constant was varied from $0.1$ to $0.05$. This means that the change of biological populations in the system will not be too drastic. This helps to maintain biodiversity in the ecosystem so that various biological populations can coexist in a relatively stable manner.

7.   Conclusions

This study investigated a discrete prey-predator system that incorporates a fear factor, a strong Allee effect, and prey refuge. Employing the forward Euler scheme on the system (2.2) yielded the discrete system (2.4). Notably, the fixed points of the discrete system (2.4) align with those of its continuous counterpart (2.2). However, despite this similarity, the dynamic behaviors of systems (2.2) and (2.4) diverge significantly. We performed a comprehensive examination of the local stability from the perspective of the fixed points, delving meticulously into the specifics of local bifurcations that occur at the interior fixed point. Our findings revealed that system (2.4) undergoes both flip and Neimark-Sacker bifurcations. Under the conditions of system (2.2), the fixed point $(0, 0)$ consistently maintains stability. However, within the framework of system (2.4), it manifests complex dynamics. Furthermore, the topological classification of fixed points in the system (2.4) is contingent upon the chosen step size, denoted as $h$. Our numerical simulations demonstrated that flip and Neimark-Sacker bifurcations manifest when a large step size is employed in Euler's method. In contrast, the dynamics of the system (2.2) remain independent of $h$. Consequently, we have compelling grounds to assert that the dynamic behavior of the system (2.4) is more intricate than that of the system (2.2).

We also discuss the impacts of the fear effect, refuge, and the Allee effect on the system. Decreasing the fear or Allee effect parameters delays the Neimark-Sacker bifurcation and results in an increase in the population densities of prey and predators at the stable fixed point. Conversely, reducing prey refuge size parameters advances the Neimark-Sacker bifurcation and leads to a decrease in the population densities of prey and predators at the stable fixed point. Hence, all three effects significantly influence the system.

In addition to investigating the dynamic behavior of discrete predator-prey systems under various biological effects, it is imperative to understand the biological significance of chaos control. Effectively regulating chaos in ecological systems holds profound implications for ecosystem stability and biodiversity conservation. Chaotic dynamics in predator-prey interactions can lead to unpredictable population fluctuations and destabilize ecosystems, potentially threatening the survival of species and disrupting ecological balance. By applying these three chaos control methods, i.e., state feedback control, OGY feedback control, and hybrid control, we effectively managed the chaotic behavior and bifurcations generated by the system (2.4). This intervention allowed us to mitigate the adverse effects of chaos and bifurcations, consequently enhancing ecosystem resilience. Furthermore, this approach not only advances our comprehension of ecological processes, it also offers invaluable insights for sustainable ecosystem management and biodiversity preservation.

To validate and illustrate our theoretical analysis, we conducted numerical simulations. The outcomes of these simulations are presented in the forms of bifurcation diagrams, phase portraits, and time-series plots, providing a comprehensive visualization of the system's dynamic behavior as discussed in the theoretical framework. Furthermore, in the numerical simulations, we observed that feedback control methods and the OGY control method are simple and intuitive, allowing for the control of chaotic systems by adjusting parameters or applying external forcing. However, the abundance of parameters and their interdependence necessitate meticulous calibration. In contrast, hybrid control, which does not require precise system models, combines the advantages of parameter perturbation control and state feedback control, thus enhancing adaptability to chaotic systems.

In future work, the model can be discretized by using different discretization methods. In addition, new discretization parameters could be chosen for the study in order to obtain specific biological properties to explore the influence of different ecological effects on the predator-prey model.

Appendix A

Here are some useful lemmas for novice authors to include in their articles.

Lemma 7.1. (Jacobian matrix ^[45]) Let ${\mathbf f}({\mathbf x}) = [f_1({\mathbf x}), f_2({\mathbf x}), ..., f_m({\mathbf x})]^T$ be a vector-valued function of $n$ variables $x_1, x_2, ..., x_n$, where ${\mathbf x} = [x_1, x_2, ..., x_n]^T$ is an $n$-dimensional vector and $f_1, f_2, ..., f_m$ are real-valued functions. If each component function $f_1, f_2, ..., f_m$ has continuous first-order partial derivatives at some point, then the Jacobian matrix $J$ is an $m \times n$ matrix, where the element in the $i$-th row and $j$-th column is the partial derivative of the component function $f_i$ with respect to the variable $x_j$, given by

Assume the following discrete system:

where $f:G\to {{R}^{n}}$ is a continuous smooth map from the ${{R}^{n}}$ open region $G$ to ${{R}^{n}}$. The mapping (7.1) can be transformed by a linear transformation into the following form

where $B$ has eigenvalues on the unit disk and all $C$ eigenvalues are either inside or outside of the circle.

Lemma 7.2. ^[46] The localized center manifold theorem for map (7.2) can be expressed as follows:

The center manifold provides an effective way to reduce the dimensionality of the system, thereby simplifying the description and analysis of the system. By studying the dynamic behavior of the system on the center manifold, we can better understand the stability and dynamical characteristics of the system near fixed points.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (No. 61074005) and the Talent Project of the High Education of Liaoning Province (LR2012005), as well as by the Youth Research Project of the Educational Department of Liaoning Province (JYTQN2023284) and the Science and Technology Project of the Science of Technology Department of Liaoning Province (2023BSBA241).

Conflict of interest

The authors declare no conflict of interest.