Dynamics of a nonlinear discrete predator-prey system with fear effect

1.   Introduction

Populations in nature are rarely isolated, and intensively interact with others in the biological community. All kinds of organisms are divided into different levels according to their physiological characteristics, food sources, etc. and different levels of populations have a variety of connections. The predator-prey process is the most fundamental, important and universal process in the study of population dynamic behavior. Many researchers have studied the dynamic behavior of many prey-predator systems in ecology and behavioral phenomena between species ^{[1,2,3,4,5,6,7,8,9,10,11]}. Some authors have also explored the complexity, stability, and conditional requirements for spatial pattern formation in prey-predator systems ^[12,13,14].

Numerous studies have shown that discrete-time systems are more suitable than the continuous system of small populations, and provide valid evidence for these ^{[15,16,17,18]}. Cheng et al. ^[19] studied a discrete-time ratio-dependent prey-predator system with Allee effect, and obtained the model with logistic growth function that have somewhat similar bifurcation structures. In the past few years, a large number of research studies have indicated that discrete prey-predator systems have more abundant dynamic behaviors than continuous systems, such as chaos. Scientists have also analyzed the corresponding dynamic behaviors between populations by numerical simulation ^{[20,21,22,23,24,25,26,27,28,29,30]}.

Qamar Din ^[31] studied the following discrete-time system:

where $r, K, \beta, \gamma, a, b, c, n$ and $d$ are greater than zero, $r$ is the intrinsic growth rate of the prey $u$ population, $K$ denotes environmental carrying capacity of the prey $u$ in a particular habitat, $\beta$ and $a$ represent the maximum value of the per capita reduction rate of the prey $u$ and predator $v$, respectively. $\gamma$ and $c$ indicate the extent to which the environment provides protection to prey and predator. $b$ denotes the quality of food that the prey provides for conversion into predator births, and $d$ measures the death rate of the predator. $n$ stands for time.

In 2016, Wang et al. ^[21] showed, through experiments, that prey$'$s fear of predators would lead to a decrease in the birth rate of prey, and $F(k, v) = \frac{1}{1+kv}$ was used to denote the fear factor. Here, $k$ reflects the degree of fear that drives prey anti-predator behavior. In the past, many researchers have only studied the effects of direct killing, no matter how they improve the predator-prey model. In this paper, we combine fear (indirect effects) and investigate the effects of fear on population dynamics.

To study the effects of fear on population dynamics, on the basis of system (1.1), we introduce the fear factor $F(k, v) = \frac{1}{1+kv}$ and the growth rate $\alpha$ of the predator $v$ population, and consider the discrete-time predator-prey system:

This article is organized as follows: in Section 2, the existence, stability and boundedness of the system at different equilibrium points are analyzed. In Section 3, we discuss the specific conditions for the existence of Neimark-Sacker bifurcation and flip bifurcation. In Section 4, chaos is controlled by the feedback control method. In Section 5, we carry out numerical simulations, including the bifurcation diagrams, phase portraits and solution diagrams. Finally, a brief conclusion is given in the last section.

2.   The properties of equilibrium points

In this section, we consider the discrete-time system (1.2) in the closed first quadrant $\mathbb{R}_{+}^{2}$ of the $(u, v)$ plane. We study the existence, stability and boundedness of the equilibrium points by the eigenvalues for the Jacobian matrix of (1.2) at the equilibrium points.

2.1.   Existence and stability

To obtain the equilibrium points of (1.2), we calculate the following equations:

Through calculation, the following results can be gained directly:

Proposition 1. (i) For all parameter values, system (1.2) has two equilibrium points $H_{0} = (0, 0), \; H_{1} = (K, 0)$;

(ii) If $\alpha > d$, then system (1.2) has a boundary positive equilibrium point $H_{2} = (0, \frac{(\alpha-d)c}{a});$

(iii) System (1.2) has a unique positive equilibrium point $H_{3} = (u^{*}, v^{*}) = (u^{*}, \frac{1}{a}(\alpha-d)(bu^{*}+c))$, where $\alpha > d$, and $u^{*}$ is the only positive solution to the quadratic equation of one variabie

where

$C_{0} = \frac{Kk\beta}{a^{2}}(\alpha-d)^{2}b^{2}+r,$

$C_{1} = \frac{2bcKk\beta}{a^{2}}(\alpha-d)^{2}+\frac{bK\beta}{a}(\alpha-d)+r\gamma-rK,$

$C_{2} = \frac{Kk\beta}{a^{2}}(\alpha-d)^{2}c^{2}+\frac{cK\beta}{a}(\alpha-d)-rK\gamma.$

Definition 1. ^[11] Suppose that $\lambda_{1}$ and $\lambda_{2}$ are two roots of the characteristic equation $F(\lambda) = \lambda^{2}+M\lambda+N = 0$, where $M$ and $N$ are constants. Then equilibrium point (u, v) is called

(i) sink if $|\lambda_{1}| < 1$ and $|\lambda_{2}| < 1$, and it is locally asymptotically stable;

(ii) source or repeller if $|\lambda_{1}| > 1$ and $|\lambda_{2}| > 1$, and it is locally unstable;

(iii) saddle if either ($|\lambda_{1}| < 1$ and $|\lambda_{2}| > 1)$ or ($|\lambda_{1}| > 1$ and $|\lambda_{2}| < 1$);

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

The Jacobian matrix for equilibrium point $H_{0}(0, 0)$ is:

Then, $\lambda_{1} = e^{r}, \; \lambda_{2} = e^{\alpha-d}$. Thus, the following proposition holds.

Proposition 2. Equilibrium point $H_{0}(0, 0)$ is

(i) source and it is locally unstable if $\alpha > d$;

(ii) saddle if $\alpha < d$.

Proof. According to (2.1), the two eigenvalues of (1.2) at the equilibrium point $H_{0}(0, 0)$ are $\lambda_{1} = e^{r}, \; \lambda_{2} = e^{\alpha-d}$. If $\alpha-d > 0$ and $r > 0$, then $\left|\lambda_{1}\right| > 1$, $|\lambda_{2}| > 1$. Thus from Definition 1, $H_{0} = (0, 0)$ is a source. If $\alpha < d$, then $0 < |\lambda_{2}| < 1$. $H_{0} = (0, 0)$ is a saddle. This completes the proof. □

For equilibrium point $H_{1} = (K, 0)$, the Jacobian matrix is described as follows:

the corresponding characteristic roots are $\lambda_{1} = 1-r, \; \lambda_{2} = e^{\alpha-d}$. Thus, the following proposition holds.

Proposition 3. The eigenvalues at the boundary equilibrium point $H_{1} = (K, 0)$ are $\lambda_{1} = 1-r, \; \lambda_{2} = e^{\alpha-d}$, then

(i) $H_{1} = (K, 0)$ is sink, if $0 < r < 2$ and $\alpha-d < 0$;

(ii) $H_{1} = (K, 0)$ is saddle, if one of the following conditions is true:

    (ii-1) $\alpha-d > 0$ and $0 < r < 2;$

    (ii-2) $\alpha-d < 0$ and $r > 2;$

(iii) $H_{1} = (K, 0)$ is non-hyperbolic, if either $r = 2$ or $\alpha-d = 0;$

(iv) $H_{1} = (K, 0)$ is source, if $r > 2$ and $\alpha-d > 0.$

Proof. According to (2.2), the two eigenvalues of (1.2) at the boundary equilibrium point are $\lambda_{1} = 1-r, \; \lambda_{2} = e^{\alpha-d}$. If $\alpha-d$ are greater than zero, then $|\lambda_{2}| > 1$. Thus from Definition 1, when $\left|\lambda_{1}\right| < 1$, then $0 < r < 2$. Thus, $H_{1} = (K, 0)$ is a saddle. Similarly, when $\left|\lambda_{1}\right| > 1$, then $r > 2$, $H_{1} = (K, 0)$ is a source. Similarly, we can prove (i), (iii) by the same way. □

For equilibrium point $H_{2} = (0, \frac{(\alpha-d)c}{a})$, the Jacobian matrix is evaluated as follows:

the corresponding characteristic roots are $\lambda_{1} = \exp[\frac{ar}{a+kc(\alpha-d)}-\frac{\beta(\alpha-d)c}{a\gamma}], \; \lambda_{2} = 1-(\alpha-d)$. Thus, the following proposition holds.

Proposition 4. The eigenvalues of $Q_{H_{2}}$ are $\lambda_{1} = \exp[\frac{ar}{a+kc(\alpha-d)}-\frac{\beta(\alpha-d)c}{a\gamma}]$ and $\lambda_{2} = 1-(\alpha-d)$, then

(i) $H_{2} = (0, \frac{(\alpha-d)c}{a})$ is sink if $0 < \alpha-d < 2$ and $0 < r < \frac{\beta c(\alpha-d)[a+kc(\alpha-d)]}{a^{2}\gamma};$

(ii) $H_{2} = (0, \frac{(\alpha-d)c}{a})$ is source if $\alpha-d > 2$ and $r > \frac{\beta c(\alpha-d)[a+kc(\alpha-d)]}{a^{2}\gamma};$

(iii) $H_{2} = (0, \frac{(\alpha-d)c}{a})$ is saddle if one of the following conditions is true:

    (iii-1) $\alpha-d = 2$ and $0 < r < \frac{\beta c(\alpha-d)[a+kc(\alpha-d)]}{a^{2}\gamma};$

    (iii-2) $0 < \alpha-d < 2$ and $r > \frac{\beta c(\alpha-d)[a+kc(\alpha-d)]}{a^{2}\gamma};$

(iv) $H_{2} = (0, \frac{(\alpha-d)c}{a})$ is non-hyperbolic if either $\alpha-d = 2$ or $r = \frac{\beta c(\alpha-d)[a+kc(\alpha-d)]}{a^{2}\gamma}.$

Proof. (i) According to (2.3), the two eigenvalues of (1.2) at the boundary equilibrium point $H_{2}$ are $\lambda_{1} = \exp[\frac{ar}{a+kc(\alpha-d)}-\frac{\beta(\alpha-d)c}{a\gamma}], \; \lambda_{2} = 1-(\alpha-d)$. $H_{2} = (0, \frac{(\alpha-d)c}{a})$ is sink if and only if $\left|\lambda_{1}\right| < 1$ and $\left|\lambda_{2}\right| < 1$. When $\left|\lambda_{1}\right| < 1$, then $0 < \alpha-d < 2$. When $\left|\lambda_{2}\right| < 1$, then $0 < r < \frac{\beta c(\alpha-d)[a+kc(\alpha-d)]}{a^{2}\gamma}$. In conclusion, $H_{2} = (0, \frac{(\alpha-d)c}{a})$ is sink if $0 < \alpha-d < 2$ and $0 < r < \frac{\beta c(\alpha-d)[a+kc(\alpha-d)]}{a^{2}\gamma};$ Similarly, Proposition 4 (ii)–(iv) can be proved. □

Lemma 1. ^[11] Suppose that $F(\lambda) = \lambda^{2}-M\lambda+N,$ and $F(1) > 0$, $\lambda_{1}$ and $\lambda_{2}$ are roots of $F(\lambda) = 0$. Then the following results hold true:

(i) $\left|\lambda_{1}\right| < 1$ and $\left|\lambda_{2}\right| < 1$ if and only if $F(-1) > 0$ and $N < 1;$

(ii) $\left|\lambda_{1}\right| < 1$ and $\left|\lambda_{2}\right| > 1$ (or $\left|\lambda_{1}\right| > 1$ and $\left|\lambda_{2}\right| < 1$) if and only if $F(-1) < 0;$

(iii) $\left|\lambda_{1}\right| > 1$ and $\left|\lambda_{2}\right| > 1$ if and only if $F(-1) > 0$ and $N > 1;$

(iv) $\lambda_{1} = -1$ and $\left|\lambda_{2}\right|\neq1$ if and only if $F(-1) = 0$ and $N\neq0, 2;$

(v) $\lambda_{1}$ and $\lambda_{2}$ are complex and $\left|\lambda_{1}\right| = \left|\lambda_{2}\right| = 1$ if and only if $M^{2}-4N < 0$ and $N = 1.$

The Jacobian matrix $Q(u^{*}, v^{*})$ relevant system (1.2) at the positive equilibrium point $H_{3}(u^{*}, v^{*})$ is as follows:

Then the characteristic equation related to $Q_{H_{3}}$ is

where

$m(u^{*}, v^{*}) = trQ_{H_{3}} = 2-\frac{ru^{*}}{(1+kv^{*})K}+\frac{\beta u^{*}v^{*}}{(u^{*}+\gamma)^{2}}-(\alpha-d)$

and

$n(u^{*}, v^{*}) = detQ_{H_{3}} = [1-\frac{ru^{*}}{(1+kv^{*})K}+\frac{\beta u^{*}v^{*}}{(u^{*}+\gamma)^{2}}][1-(\alpha-d)]+[(1-\frac{u^{*}}{K})\frac{rku^{*}}{(1+kv^{*})^{2}}+\frac{\beta u^{*}}{u^{*}+\gamma}]\frac{abu^{*2}}{(bu^{*}+c)^{2}}.$

Therefore

Using Lemma 1, we have the following proposition:

Proposition 5. Let $H_{3}(u^{*}, v^{*})$ be the unique positive equilibrium point of system (1.2), then the following propositions hold:

(i) $H_{3}(u^{*}, v^{*})$ is sink if and only if $|m(u^{*}, v^{*})| < 1+n(u^{*}, v^{*}) < 2;$

(ii) $H_{3}(u^{*}, v^{*})$ is saddle if and only if $m^{2}(u^{*}, v^{*}) > 4n(u^{*}, v^{*})$ and $|m(u^{*}, v^{*})| > |1+n(u^{*}, v^{*})|$;

(iii) $H_{3}(u^{*}, v^{*})$ is source if and only if $|n(u^{*}, v^{*})| > 1$ and $|m(u^{*}, v^{*})| < |1+n(u^{*}, v^{*})|$;

(iv) $H_{3}(u^{*}, v^{*})$ is non-hyperbolic if and only if $|m(u^{*}, v^{*})| = |1+n(u^{*}, v^{*})|$ or $n(u^{*}, v^{*}) = 1$ and $|m(u^{*}, v^{*})|\leq2.$

Proof. (i) According to Lemma 1, $H_{3}(u^{*}, v^{*})$ is a sink point if and only if $F(1) > 0, F(-1) > 0$ and $N < 1$, it can be acquired by calculation $|m(u^{*}, v^{*})| < 1+n(u^{*}, v^{*}) < 2.$ Consequently, Proposition 5 (i) holds. Similarly, Proposition 5 (ii)–(iv) can be established. □

2.2.   Boundedness

Lemma 2. ^[20] Assume that $u_{t}$ satisfies $u_{0} > 0$, and $u_{t+1}\leq u_{t} \exp[A(1-B u_{t})]$ for $t\in[t_{1}, \infty)$, where $B$ is a positive constant. Then $\lim\limits_{t \to \infty} \sup\ u_{t} \leq \frac{1}{A B} \exp (A-1)$.

Theorem 1. Every positive solution $\left\{\left(u_{n}, v_{n}\right)\right\}$ of system (1.2) is uniformly bounded.

Proof. Suppose that $\left\{\left(u_{n}, v_{n}\right)\right\}$ be an arbitrary positive solution corresponding to system (1.2). Then, by the first part of (1.2), it is known

for all $n = 0, 1, 2, \cdots$. Suppose that $u_{0} > 0$, then according to Lemma 2, we gain

From the second part of (1.2), we acquire

Assume that $v_{0} > 0$, then using Lemma 2, we gain

That is to say that $\lim\limits_{n \to \infty} \sup\ (u_{n}, v_{n}) \leq M,$ where $M = \max\left\{M_{1}, M_{2}\right\}$. This completes the proof. □

3.   Bifurcation analysis

3.1.   Flip bifurcation

The characteristic equation related to system (1.2) at the positive interior equilibrium point $H_{3}$ is

where

Assume that $m^{2}(u^{*}, v^{*}) > 4n(u^{*}, v^{*})$, that is,

and $m(u^{*}, v^{*})+n(u^{*}, v^{*}) = -1,$ that is to say

Then eigenvalue of $F(\lambda) = 0$ are $\lambda_{1} = -1$ and $\lambda_{2} = 2+\Phi+\Theta-\frac{ru^{*}}{(1+kv^{*})K}$. The condition $\left|\lambda_{2}\right|\neq1$ indicates that

Consider the following set

Based on the above analysis, we can obtain that when the parameters change on set $A_{1}$, system (1.2) will occur flip bifurcation at $H_{3}(u^{*}, v^{*}).$

We consider the following system

here $(a, b, c, d, K, r_{1}, k, \alpha, \beta, \gamma)\in A_{1}$.

Consider a perturbation corresponding to system (3.5) as follows:

where $\overline{r}$ is a small perturbation parameter and $\left|\overline{r}\right|\ll1.$

Let $p = u-u^{*}$ and $q = v-v^{*}.$ Then we gain

where

We construct a nonsingular matrix $D_{1}$ and translate it as follows:

where

Taking $D_{1}^{-1}$ on both sides of Eq (3.8), we obtain

where

Applying the center manifold theorem $W^{c}(0)$ of system (3.9) at the trivial equilibrium point (0, 0) in a limited field of $\overline{r} = 0.$ Then there exists a center manifold $W^{c}(0)$ as follows:

and satisfies

and we have

Therefore, we consider the map restricted to the center manifold $W^{c}(0)$ as below:

where

According to flip bifurcation, we define the following two nonzero real numbers $\delta_{1}$ and $\delta_{2}$, where

From the above analysis, we get the following theorem:

Theorem 2. If $\delta_{1}\neq0, \delta_{2}\neq0$, then system (1.2) passes through a flip bifurcation at the fixed point $H_{3}(u^{*}, v^{*})$ when the parameter $r$ alters in the small region of $r_{1}$. In addition, if $\delta_{2} > 0$ (resp., $\delta_{2} < 0$), then the period-two orbits that bifurcate from fixed point $H_{3}(u^{*}, v^{*})$ are stable (resp., unstable).

3.2.   Neimark-Sacker bifurcation

When the parameters change on set $A_{2}$, system (1.2) will undergo Neimark-Sacker bifurcation at the unique positive interior equilibrium point $H_{3}(u^{*}, v^{*})$, where

Consider a perturbation related to system (1.2) as follows:

where $\overline{r}$ is a limited perturbation parameter and $\left|\overline{r}\right|\ll1$.

The characteristic equation of system (3.10) at $H_{3}(u^{*}, v^{*})$ is as follows:

where

and

Since parameters $(a, b, c, d, K, \alpha, \beta, \gamma, k, r_{2})\in A_{2}$, the characteristic values of system (3.10) at $H_{3}(u^{*}, v^{*})$ are a pair of complex conjugate numbers $\lambda$ and $\bar{\lambda}$ with $|\lambda| = |\bar{\lambda}| = 1$ as follows

Therefore we have

When $\overline{r}$ changes in limited field of $\overline{r} = 0,$ then $\lambda, \; \bar{\lambda} = x\pm iy,$ where

In addition, Neimark-Sacker bifurcation requires that $\overline{r} = 0$, $\lambda^{z}, \; \bar{\lambda}^{z}\neq1$ ($z$ = 1, 2, 3, 4), which is equivalent to $m(0)\neq$ -2, 0, -1, 2. Because parameters $(a, b, c, d, K, \alpha, \beta, \gamma, k, r_{2})\in A_{2}$, therefore $m(0)\neq$ -2, 2. We only require $m(0)\neq$ 0, -1, so that

Let $p = u-u^{*}$ and $q = v-v^{*}.$

After the transformation of the equilibrium point $H_{3}(u^{*}, v^{*})$ of system (3.10) to the origin, we have

where

and $S_{11}, S_{12}, S_{13}, S_{14}, S_{15}, T_{11}, T_{12}, T_{13}, T_{14}, S_{21}, S_{22}, S_{23}, S_{24}, S_{25}, T_{21}, T_{22}, T_{23}, T_{24},$ are given in (3.7) by substituting $r_{1}$ for $r_{2}+\overline{r}$.

Besides that, we analyse the normal form of system (3.12) when $\overline{r} = 0.$

Consider the translation as follows:

where

Taking $D_{2}^{-1}$ on both sides of system (3.12), we acquire

where

System (1.2) occurs the Neimark-Sacker bifurcation if the following quantity $\vartheta$ is not zero,

where

If $\vartheta\neq0$, Neimark-Sacker bifurcation will undergo in system (1.2), and the following theorem holds:

Theorem 3. System (1.2) undergoes a Neimark-Sacker bifurcation at the positive equilibrium point $H_{3}(u^{*}, v^{*})$ if conditions (3.11) are satisfied and $\vartheta\neq0$. In addition, if $\vartheta > 0$ (resp., $\vartheta < 0$), then an repelling (resp., attracting) invariant closed curve bifurcates from fixed point $H_{3}(u^{*}, v^{*})$ for $r < r_{2}$ (resp., $r > r_{2}$).

4.   Chaos control

In this section, we will adopt the feedback control method ^[26,27,28] to stabilize the chaotic orbit at an unstable equilibrium point by adding a feedback control term to the system (1.2). Therefore, system (1.2) makes the following form:

where $x(u_{n}, v_{n}) = h_{1}(u_{n}-u^{*}) + h_{2}(v_{n}-v^{*})$ is feedback controlling force, $h_{1}$ and $h_{2}$ are feedback gains, and $(u^{*}, v^{*})$ is the unique positive equilibrium point of (1.2). Furthermore, $f(u^{*}, v^{*}) = u^{*}$, and $g(u^{*}, v^{*}) = v^{*}$.

The Jacobian matrix of system (4.1) at positive equilibrium point $(u^{*}, v^{*})$ is as follows:

where

Thus, the characteristic equation related to $J(u^{*}, v^{*})$ is:

Let $\lambda_{1}$ and $\lambda_{2}$ be the eigenvalues of characteristic equation (4.2), then

Next, we must solve equations $\lambda_{1} = \pm1$ and $\lambda_{1}\lambda_{2} = 1$ to gain the critical stability line. At the same time, it also ensures that the absolute value $\lambda_{1}$ and $\lambda_{2}$ are less than one.

Suppose that $\lambda_{1}\lambda_{2} = 1$, then we gain

Assume that $\lambda_{1} = 1$, then we have

Assume that $\lambda_{1} = -1$, then we obtain

Thus, the stable eigenvalues lie within the triangular area with the boundaries of the straight lines $L_{1}, L_{2}, L_{3}$. In addition, when the control parameters $h_{1}$ and $h_{2}$ take values in the triangular region, system (4.1) will not generate chaos.

5.   Numerical simulations

In this section, we draw the bifurcation diagrams, phase portraits, solution of the figures and maximum Lyapunov exponents for system (1.2) to verify the above theoretical analysis and show the new interesting complex dynamical behaviors and the stability of the predator-prey system at the equilibrium point by using numerical simulations.

5.1.   System (1.2) without fear factor ($k$ = 0)

First, in Figure 1, we consider that the fear factor $k = 0$ and take $r$ as the bifurcation parameter to discuss the dynamic behavior of (1.2) at $H_{3}(u^{*}, v^{*})$. We consider the parameter values as $(a, b, c, d, \alpha, \beta, \gamma, K) = (1.8, 2.8, 3.5, 0.6, 1.27, 1.1, 2.7, 1.5)\in A_{1}$ with the initial value of $(u_{0}, v_{0})$ = (2, 1) and $r \in$ [2.8, 4.6]. Flip bifurcation emerges from the unique positive equilibrium point and loses its stability as $r$ goes through a critical value $r = 3.102$, and it is stable when $r < 3.102$, and when $r > 3.102$, system (1.2) oscillates with periods of $2, 2^{2}, 2^{3}, \cdots.$ It can be obtained from Figure 1(c) and Figure 2(a–c) that chaos will happen in system (1.2) as the bifurcation parameters $r$ continue to increase.

In Figure 3, taking $(a, b, c, d, \alpha, \beta, \gamma, K) = (1.7, 1.6, 3.1, 0.01, 1.2, 2.6, 1.2, 3.5)\in A_{2}$ with the initial value of $(u_{0}, v_{0})$ = (1, 2) and $r \in$ [4.6, 6.2]. Neimark-Sacker bifurcation emerges from the unique positive equilibrium point and loses its stability as $r$ goes through a critical value $r = 5.05$. It can be seen from Figure 3(a) that when $r < 4.7,$ the equilibrium point of (1.2) with respect to prey does not exist. Not only does the system have an attractive-invariant loop and periodic solutions, but it also exhibits dynamical chaos as the bifurcation parameters $r$ continue to increase. Figure 3(c) is the maximum Lyapunov exponent diagram related to Figure 3(a, b). It can be seen from the MLE that system (1.2) will appear chaotic. By observing the Figure 4(a–c), it can be found that when $r > 5.05$, a limit cycle, a periodic window and chaos appear in system (1.2).

5.2.   System (1.2) with fear factor $(k > 0)$

In Figure 5, we consider that the fear factor $k > 0$. Taking $(a, b, c, d, \alpha, \beta, \gamma, k, K) = (2, 2.8, 3.5, 0.6, 1.1, 1.1, 3, 0.5, 2)\in A_{1}$ with the initial value of $(u_{0}, v_{0})$ = (2, 1) and $r\in$ [5.2, 7]. Flip bifurcation appears from the unique positive equilibrium point and loses its stability as $r$ goes through a critical value $r = 5.63$, and it is stable when $r < 5.63$ and when $r > 5.63$, system (1.2) oscillates with periods of $2, 2^{2}, 2^{3}, \cdots.$ It can be acquired from Figure 5(c) that chaos will happen in system (1.2) as the bifurcation parameters $r$ continue to increase.

In Figure 6, taking $(a, b, c, d, \alpha, \beta, \gamma, k, K) = (1.7, 1.6, 3.1, 0.01, 1.2, 2.6, 1.2, 0.3, 3.5)\in A_{2}$ with the initial value of $(u_{0}, v_{0})$ = (1, 2) and $r \in$ ^[8,11]. Neimark-Sacker bifurcation emerges from the unique positive equilibrium point and loses its stability as $r$ goes through a critical value $r = 8.63$. We notice that the equilibrium point of (1.2) is stable for $r < 8.63$, loses its stability at $r = 8.63$ and not only a limit cycle but also periodic solution emerge when the bifurcation parameter $r > 8.63$. Other than that, the value of the MLE related to (1.2) is greater than 0 as $r$ continues to increase, and thus chaos will occur, i.e., the solution of (1.2) is arbitrarily periodic.

In Figure 7, taking $(a, b, c, d, \alpha, \beta, \gamma, r, K) = (1.5, 2.8, 3.5, 0.6, 1, 1.1, 3, 2, 1.5)$ with the initial value of $(u_{0}, v_{0})$ = (2, 1) and $k \in$ [0, 6], $k$ is a bifurcation parameter. At this time, the bifurcation phenomenon of (1.2) will not occur. The population density of prey and predator will continue to decrease and tend to 0 with the increase of fear factor $k$. It is important to note that the cost of fear does not lead to the extinction of predators, but rather to the extinction of prey.

5.3.   Controlling chaos

In Figure 8, when the parameter value is $(a, b, c, d, \alpha, \beta, \gamma, k, r, K) = (1.7, 1.6, 3.1, 0.01, 1.2, 2.6, 1.2, 0.3, 10.8, 3.5)$ with the initial value of $(u_{0}, v_{0})$ = (1, 2). In Figure 6(c) when the bifurcation parameter $r = 10.8$, system (1.2) will produce chaos. When the $h_{1}$ and $h_{2}$ are controlled in the triangular region surrounded by three straight lines $L_{1}$, $L_{2}$, and $L_{3}$, the chaos generated by system (4.1) will be controlled near the equilibrium point and become an asymptotically stable state.

6.   Conclusions

Studies have shown that discrete systems have richer and more complex dynamic behaviors than continuous systems. Hence, on the basis of previous research work, this paper discusses the stability, bifurcation and chaos control of a nonlinear discrete prey-predator system with fear effect. Based on the results of the study, we can draw the following conclusions:

(a) System (1.2) has four equilibrium points, where the stable equilibrium point is positive, and depicts the coexistence of prey and predators.

(b) System (1.2) has flip bifurcation and Neimark-Sacker bifurcation happen at the positive interior equilibrium point when $r$ alters in $A_{1}$ and $A_{2}$ small fields. (see Figures 1, 3, 5, 6). We can also observe the orbits of periods 2, 4, and 8 periodic windows of flip bifurcation.

(c) When $k = 0$, system (1.2) at the positive equilibrium point will generate the Neimark-Sacker bifurcation, flip bifurcation and chaos as the bifurcation parameters $r$ continue to increase.

(d) When the fear parameter $k$ is greater, both predators and prey populations decrease. It is important to note that the cost of fear does not lead to the extinction of predators, but rather to the extinction of prey. (see Figure 7).

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 National Natural Science Foundation of China (No.12161079), the Northwest Normal University Graduate Research Grant Project (No.2022KYZZ-S114), and the Gansu Province Innovation Star Project (No.2023CXZX-325).

Conflict of interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.