Bifurcation analysis and chaos control of a discrete-time prey-predator model with fear factor

1.   Introduction

In nature, various organisms are divided into different levels according to their different physiological characteristics and food sources, and there are various connections between populations at different levels. The prey-predator process is the most basic and important process in population dynamics. Under normal circumstances, on the one hand, predators can directly affect the population number of prey, and on the other, they can also affect the population change of predators themselves. The two relationships are very complex ^[1,2].

Mathematical modeling is an important method of scientific research. Although mathematical modeling can not completely solve the specific problems in biology, it plays an important role in describing the change law of biological population. For example, in population dynamics, the continuous model can describe the population growth law when the population number is relatively large or the generations overlap ^{[3,4,5,6,7,8]}. However, when the generations of the population do not overlap each other and their growth is discontinuous and step-by-step, the discrete model can more truly describe the change law of the population than the continuous model.

In 1976, May showed that first-order difference equations have a series of surprising dynamic behaviors in biological and economic contexts ^[9]. Thus, in recent years more and more literatures pay attention to the dynamical behavior of discrete-time systems ^{[10,11,12,13,14,15,16,17]}. Cheng et al. ^[18] studied the influence of Allee effect on the dynamic behavior of discrete predator-prey model, and analyzed the asymptotic behavior and bifurcation structure at the positive equilibrium point. He et al. ^[19] used the central manifold theorem and bifurcation theory to show that discrete predator-prey systems undergo flip bifurcation and Neimark-Sacker bifurcation in the interior of $R_{+}^{2}$. In addition, they had stabilized the chaotic orbits at an unstable fixed point using the feedback control method.

In 2016, Wang et al.^[20] showed through experiments that prey fear of predators would lead to lower birth rate of prey, and $F(k, y) = \frac{1}{1+ky}$ was used to represent the fear factor. Here, the parameter $k$ reflects the level of fear which drives anti-predator behavior of the prey.

Based on the discussion above, in this paper, we consider the following discrete model

where $r$, $k$, $K$, $a$, $b$ and $s$ are positive, $r$ and $s$ are the intrinsic growth rates of the prey $x$ and predator $y$ populations, respectively. $K$ is the carrying capacity for $x$. The functional response $d(x, y) = \frac{bx}{x+a}$ depending on $x$ only. The constant $b$ is the maximum of the predation rate when the predator will not or cannot kill more prey even when the latter is available. The constant $a$ refers to some value of the prey population beyond which the predators attacking capability begins to saturate ^[21,22].

The rest of this paper is organized as follows. In Section 2, the existence and stability of fixed points are analyzed. Section 3 discusses the existence of flip bifurcation and Neimark-Sacker bifurcation. In Section 4, chaos is controlled to an unstable fixed point using the feedback control method. In Section 5, we perform numerical simulations which include the bifurcation diagrams and the phase portraits. Finally, in Section 6, we will analyze and summarize our conclusions.

2.   The existence and stability of fixed points

In this section, we will calculate the fixed points of system (1.1) and give the conditions for the asymptotic stability of the fixed point. To find the fixed points of system (1.1), we assume that:

By solving the following system

we can get the following proposition.

Proposition 1.

$(A_{1})$ System (1.1) has a boundary equilibrium point $E_{0}(K, 0)$;

$(A_{2})$ System (1.1) has a unique positive equilibrium point $E^{*}(x^{*}, y^{*})$, where

and

The Jacobian matrix $J(x, y)$ corresponding to system (1.1) at point $(x, y)$ is as follows:

Assume that $\lambda_{1}$ and $\lambda_{2}$ are two roots of the characteristic equation of the Jacobian matrix $J|_{(x, y)}$, and we have the following definition and conclusions.

Definition 2.1. ^[10] The fixed point (x, y) is called

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

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

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

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

The Jacobian matrix at $E_{0}(K, 0)$ is:

Then, $\lambda_{1} = 1+s, \lambda_{2} = 1-r$. Thus, the following proposition holds.

Proposition 2. The characteristic root at the boundary equilibrium point $E_{0}(K, 0)$ is $\lambda_{1} = 1+s, \lambda_{2} = 1-r$, then

$(A_{1})$ $E_{0}(K, 0)$ is a saddle point, if $0 < r < 2$;

$(A_{2})$ $E_{0}(K, 0)$ is a source point, if $r > 2$;

$(A_{3})$ $E_{0}(K, 0)$ is a non-hyperbolic point, if $r = 2$.

Proof. It can be seen from (2.3) that the two characteristic roots of system (1.1) at the boundary equilibrium point are $\lambda_{1} = 1+s, \lambda_{2} = 1-r$. Since $s > 0$ and $r > 0$, then $|\lambda_{1}| > 1$. Thus from Definition 2.1, when $|\lambda_{2}| > 1$, $E_{0}(K, 0)$ is the source point; when $|\lambda_{2}| < 1$, $E_{0}(K, 0)$ is the saddle point. It can be known by calculation that when $0 < r < 2$, $E_{0}(K, 0)$ is the saddle point; when $r > 2$, $E_{0}(K, 0)$ is the source point; when $r = 2$, the boundary equilibrium point $E_{0}(K, 0)$ is non-hyperbolic.

$J(x, y)$ evaluated at the positive equilibrium point $E^{*}(x^{*}, y^{*})$ is

Then the characteristic equation corresponding to $J(E^{*})$ is

where

Then

Lemma 2.1. ^[10] Assume that $F(\lambda) = \lambda^{2}-A\lambda+B$, and $F(1) > 0$ with $\lambda_{1}, \lambda_{2}$ are roots of $F(\lambda) = 0$. Then the following results hold:

$(A_{1})$ $|\lambda_{1}| < 1$ and $|\lambda_{2}| < 1$ if and only if $F(-1) > 0$ and $B < 1$;

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

$(A_{3})$ $|\lambda_{1}| > 1$ and $|\lambda_{2}| > 1$ if and only if $F(-1) > 0$ and $B > 1$;

$(A_{4})$ $\lambda_{1} = 1$ and $\lambda_{2}\neq1$ if and only if $F(-1) = 0$ and $B\neq0, 2$;

$(A_{5})$ $\lambda_{1}$ and $\lambda_{2}$ are complex and $|\lambda_{1}| = 1$ and $|\lambda_{2}| = 1$ if and only if $A^{2}-4B < 0$ and $B = 1$.

The following proposition shows the local dynamics of the unique positive equilibrium point $(x^{*}, y^{*})$ from Lemma 2.1.

Proposition 3. Let $E^{*}(x^{*}, y^{*})$ be the unique positive equilibrium point of system (1.1), and the following propositions hold:

$(A_{1})$ $E^{*}(x^{*}, y^{*})$ is a sink point and it is locally asymptotically stable if and only if

$(A_{2})$ $E^{*}(x^{*}, y^{*})$ is a source point and it is locally unstable if and only if

$(A_{3})$ $E^{*}(x^{*}, y^{*})$ is a saddle point if and only if

$(A_{4})$ $E^{*}(x^{*}, y^{*})$ is non-hyperbolic if one of the following conditions holds:

or

Proof. According to Lemma 2.1, $E^{*}(x^{*}, y^{*})$ is a sink point if and only if $F(1) > 0$, $F(-1) > 0$ and $B < 1$, it can be obtained by calculation $p < |1+q|$ and $q < 1$. Therefore, Proposition 3 $(A_{1})$ holds. Similarly, Proposition 3 $(A_{2})$, $(A_{3})$ and $(A_{4})$ can be established.

3.   Bifurcation analysis

In this section, we will use the central manifold theorem and bifurcation theorem ^[23,24,25] to discuss the existence conditions and related conclusions of flip bifurcation and Neimark-Sacker bifurcation of system (1.1) at the fixed points $E_{0}(K, 0)$ and $E^{*}(x^{*}, y^{*})$.

3.1.   Bifurcation at $E_{0}(K, 0)$

The eigenvalues of $J(E_{0})$ are $\lambda_{1} = 1+s$ and $\lambda_{2} = 1-r$. According to Proposition 2, $E_{0}(K, 0)$ is non-hyperbolic if $r = 2$. If $r$ takes that value, we have $\lambda_{1} = 1+s$ and $\lambda_{2} = -1$ with $|\lambda_{1}|\neq1$. Now, let

System (1.1) around steady state $E_{0}(K, 0)$ admits a flip bifurcation when parameters vary in a small neighborhood of F.B.1. Take parameter $s$ as the bifurcation parameter. Then system (1.1) can be transformed into

where $|s^{*}|\ll 1$ is a small perturbation parameter.

Let $X = x-K, \; Y = y$, then we get

where

$A_{11} = 1-r$, $A_{12} = -\frac{bK}{K+a}$, $A_{13} = -\frac{r}{K}$, $A_{14} = rk-\frac{ab}{(K+a)^{2}}$,

$A_{15} = A_{16} = 0$, $A_{17} = \frac{rk}{K}+\frac{ab}{(K+a)^{3}}$, $A_{18} = -rk^{2}$,

$A_{21} = A_{23} = A_{24} = A_{26} = A_{27} = 0$, $A_{22} = 1+s$, $A_{25} = -\frac{s}{K}$,

$A_{28} = \frac{s}{K^{2}}$, $B_{1} = B_{3} = B_{5} = 0$, $B_{2} = 1$, $B_{4} = -\frac{1}{K}$.

We construct a nonsingular matrix

and use the following transformation for (3.2)

Then (3.2) can be changed into

where

According to the central manifold theorem, there exists a center manifold $W^{c}(0)$:

and satisfy

where $P$ and $s^{*}$ sufficiently small.

It can be obtained by calculation

Hence, the map restricted to the center manifold $W^{c}(0)$ is given by

where

Next, we define the following two nonzero real numbers:

Therefore, based on the above analysis, we can conclude the following:

Theorem 3.1. If $h_{1}\neq0, \; h_{2}\neq0$, then system (1.1) undergoes a flip bifurcation at the boundary equilibrium point $E_{0}(K, 0)$ when parameters vary in a small neighborhood of F.B.1. And when $h_{2} > 0$ (respectively, $h_{2} < 0$), system (1.1) bifurcates from the fixed point to a 2-periodic stable orbit(respectively, unstable).

3.2.   Bifurcation at $E^{*}(x^{*}, y^{*})$

3.2.1.   Flip bifurcation

We consider the following set:

Firstly, the dynamic analysis of system (1.1) is analyzed when the parameters change in the small field of F.B.2. Select parameter $(a, b, s, k, r, K)\in$ F.B.2 and consider the following system:

where $s$ is the bifurcation parameter, $|s^{*}|\ll 1$ is a small perturbation parameter.

Let $u = x-x^{*}, \; v = y-y^{*}$. Then we get

where

$C_{11} = 1+\frac{r}{1+kx^{*}}(1-\frac{2x^{*}}{K})-\frac{abx^{*}}{(x^{*}+a)^{2}}$, $C_{12} = -\frac{rkx^{*}}{(1+kx^{*})^{2}}(1-\frac{x^{*}}{K})-\frac{bx^{*}}{x^{*}+a}$,

$C_{13} = -\frac{r}{K(1+kx^{*})}+\frac{abx^{*}}{(x^{*}+a)^{3}}$, $C_{14} = -\frac{rk}{(1+kx^{*})^{2}}(1-\frac{2x^{*}}{K})-\frac{ab}{(x^{*}+a)^{2}}$,

$C_{15} = \frac{rk^{2}x^{*}}{(1+kx^{*})^{3}}(1-\frac{x^{*}}{K})$, $C_{16} = -\frac{abx^{*}}{(1+kx^{*})^{3}}(1-\frac{x^{*}}{K})$,

$C_{17} = \frac{rk}{K(1+kx^{*})^{2}}+\frac{ab}{(x^{*}+a)^{3}}$, $C_{18} = \frac{rk^{2}}{(1+kx^{*})^{3}}(1-\frac{2x^{*}}{K})$,

$C_{21} = s$, $C_{22} = 1-s$, $C_{23} = C_{24} = -\frac{s}{x^{*}}$, $C_{25} = \frac{2s}{x^{*}}$, $C_{26} = C_{28} = \frac{s}{x^{*2}}$,

$C_{27} = -\frac{2s}{x^{*2}}$, $D_{1} = 1$, $D_{2} = -1$, $D_{3} = D_{4} = -\frac{1}{x^{*}}$, $D_{5} = \frac{2}{x^{*}}$.

We construct a nonsingular matrix

and use the translation

then (3.9) can be changed into

where

According to the central manifold theorem, there exists a center manifold $W^{c}(0)$:

and satisfy

where $\widetilde{P}$ and $s^{*}$ sufficiently small.

It can be obtained by calculation

Hence, the map restricted to the center manifold $W_{1}^{c}(0)$ is given by

where

Remember

If parameters $\beta_{1}$ and $\beta_{1}$ are not 0, system (3.13) has flip bifurcation. And we can get the following theorem:

Theorem 3.2. If $\beta_{1}\neq0, \; \beta_{2}\neq0$, then system (1.1) undergoes a flip bifurcation at the positive equilibrium point $E^{*}(x^{*}, y^{*})$ when parameters vary in a small neighborhood of F.B.2. And when $\beta_{2} > 0$ (respectively, $\beta_{2} < 0$), system (1.1) bifurcates from the fixed point $E^{*}(x^{*}, y^{*})$ to a 2-periodic stable orbit(respectively, unstable).

3.2.2.   Neimark-Sacker bifurcation

Next, we study the Neimark-Sacker bifurcation of positive equilibrium point of system (1.1). The characteristic polynomial (2.4) of Jacobian matrix of linearized system of (1.1) about positive equilibrium point $(x^{*}, y^{*})$ can be rewritten as:

where

For convenience, we assume

Then

Let

Then, the dynamic analysis of system (1.1) is analyzed when the parameters change in the small field of $N.S$. Select parameter $(a, b, s, k, r, K)\in N.S$ and consider the following systems:

where $|s^{*}|\ll 1$ is a small perturbation parameter. And we choose $s$ as the bifurcation parameter.

Let $u = x-x^{*}, \; v = y-y^{*}$, then we get

where $C_{11}, C_{12}, C_{13}, C_{14}, C_{15}, C_{16}, C_{17}, C_{18}, C_{21}, C_{22}, C_{23}, C_{24}, C_{25}, C_{26}, C_{27}, C_{28}$ are given in (3.14) by substituting $s$ for $s+s^{*}$.

The characteristic equation of system (3.14) at $(u, v) = (0, 0)$ is as follows:

where

Since parameters $(a, b, r, s, k, K)\in N.S$, the roots of the characteristic equation are

and we have

In addition, it is required that $s^{*} = 0, \lambda^{j}_{1, 2} \neq 1 (j = 1, 2, 3, 4)$ which is equivalent to $M(0)\neq$ -2, 0, 1, 2. Because $(a, b, r, s, k, K)\in N.S$, thus $M(0)\neq$ -2, 2. We only require $M(0) \neq 0, 1,$ so that

Therefore, eigenvalues $\lambda_{1}, \lambda_{2}$ of the fixed point (0,0) of system (3.14) do not lay in the intersection of the unit circle with the coordinate axes when $s^{*} = 0$ and the condition (3.15) holds.

Let $\rho = -\frac{M(0)}{2}, \; \omega = \frac{\sqrt{4N(0)-M^{2}(0)}}{2}$, we use the following transformation:

and system (3.14) becomes into

where

System (3.16) undergoes the Neimark-Sacker bifurcation if the following quantity is not zero

where

Through some complicated calculations, we get

If $\mathfrak{L}\neq0$, Neimark-Sacker bifurcation will occur in system (1.1), and the following theorem holds:

Theorem 3.3. System (1.1) undergoes a Neimark-Sacker bifurcation at the positive equilibrium point $E^{*}(x^{*}, y^{*})$ if conditions in (3.15) are satisfied and $\mathfrak{L}\neq0$ in (3.17). Moreover, if $\mathfrak{L} < 0$(resp., $\mathfrak{L} > 0$), an attracting (resp., repelling) invariant closed curve bifurcates from the steady state for $s > s^{*}$ (resp., $s < s^{*}$).

4.   Chaos control

In this section, we will use the feedback control method ^[26,27,28] to control the chaos of system (1.1). Specifically, a feedback control term is added to system (1.1) to stabilize the chaotic orbit of system (1.1) at the equilibrium point. Thus, system (1.1) becomes the following form:

where

is feedback controlling force, $r_{1}$ and $r_{2}$ are feedback gains, and $(x^{*}, y^{*})$ the unique positive equilibrium point of system (1.1). Moreover $f(x^{*}, y^{*}) = x^{*}$ and $g(x^{*}, y^{*}) = y^{*}$.

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

where

$C_{11} = 1+\frac{r}{1+kx^{*}}(1-\frac{2x^{*}}{K})-\frac{abx^{*}}{(x^{*}+a)^{2}}$, $C_{12} = -\frac{rkx^{*}}{(1+kx^{*})^{2}}(1-\frac{x^{*}}{K})-\frac{bx^{*}}{x^{*}+a}$, $C_{21} = s$, $C_{22} = 1-s$.

Thus, the characteristic equation corresponding to $J^{'}(x^{*}, y^{*})$ is:

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

and

In order to make the absolute values of $\lambda_{1}$ and $\lambda_{2}$ less than 1, we assume that $\lambda_{1} = \pm1$ and $\lambda_{1}\lambda_{2} = 1$ hold.

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

Assume that $\lambda_{1} = 1$, and we get

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

Then stable eigenvalues lie within the triangular region bounded by the straight lines $L_{1}, \; L_{2}, \; L_{3}$. Therefore, when the control parameters $r_{1}$ and $r_{1}$ take values in the triangular region, system (4.1) will not produce chaos.

5.   Numerical simulations

In this section, we will use numerical simulation to verify the previous theoretical results and show the dynamic behavior of the discrete system (1.1) at the positive equilibrium point $E^{*}(x^{*}, y^{*})$. This is due to the presence of prey only at the boundary equilibrium point $E_{0}(K, 0)$ and the extinction of predators.

5.1.   Dynamical without fear factor ($k=0$)

Firstly, we assume that the fear factor $k = 0$ and take s as the bifurcation parameter to analyze the dynamic behavior of system (1.1) at the positive equilibrium point. We consider the parameter values as

and the initial value is taken as $(x_{0}, y_{0}) = (2.5, 2)$. Figure 1 is the Neimark-Sacker bifurcation diagram of system (1.1) at the positive equilibrium point. In Figure 3, system (1.1) undergoes Neimark-Sacker bifurcation when the parameter values above are taken. The critical value of $s = 1.7802$ for the bifurcation to occur can be calculated. Combined with the maximum Lyapunov exponents diagram (Figure 4(a)), when $s < 1.7802$, system (1.1) is in equilibrium, and when $s > 1.7802$, the phase diagram corresponding to system (1.1) appears closed track, and thus the periodic solution appears. However, when $s$ continues to increase, the value of the maximum Lyapunov exponents corresponding to system (1.1) is greater than 0, and thus chaos will occur, i.e., the solution of system (1.1) is arbitrarily periodic (Figure 3).

Now, we consider the parameter values as

and $k = 0$. In Figure 2, system (1.1) undergoes a period-doubling bifurcation (flip bifurcation), and it is stable when $s < 2.1142$ and when $s > 2.1142$, system (1.1) oscillates with periods of $2, 2^{2}, 2^{3}, \cdot\cdot\cdot$. It can be seen from Figure 4(b) that when the bifurcation parameters $s$ continue to increase, chaos will occur in system (1.1).

5.2.   Dynamical with fear factor ($k > 0$)

We assume that fear factor $k > 0$ and make bifurcation diagrams and maximum Lyapunov exponents diagrams (Figure 5) for $k = 0.5$ and $k = 10$ respectively with $s$ as the bifurcation parameter on the basis of (4.1). In Figures 1 and 5, when $k > 0$, the dynamic behavior of system (1.1) will change significantly. When $k = 0.5$, system (1.1) changes from Neimark-Sacker bifurcation to period-doubling bifurcation (Figure 5(a), (b)), and will produce chaos (Figure 5(c)). However, when $k = 10$, system (1.1) will produce not only period-doubling bifurcation and chaos (Figure 5(d)–(f)), but also the equilibrium point be lowered. So it can be concluded that fear factor will have a significant impact on the dynamic behavior of system (1.1).

Here, we consider the parameter values as

and fear factor $k > 0$. At this time, the bifurcation phenomenon of system (1.1) will not occur. In Figure 6, when the fear factor $k$ increases, the density of predators and preys will continue to decrease and tend to 0, which leads to the collapse of the population system and the extinction of predators and prey.

Furthermore, when the parameter value is $r = 3.5, \; b = 1.5, \; a = 1, \; K = 5, \; s = 0.5$, $k$ is bifurcation parameter. Figure 7 shows the occurrence of Neimark-Sacker bifurcation as the value of parameter $k$ varies. According to the calculation, the critical value for Neimark-Sacker bifurcation can be determined as $k_{0} = 0.1865$. As shown in the $x-y$ space, when $k < k_{0}$, the stable stationary state is stable. Moreover, the maximum Lyapunov exponents are plotted in Figure 7(c), Thus Neimark-Sacker bifurcation occurs at $k = k_{0} = 0.1865$ and [0.1865, 0.7] behaves like chaotic region as $k$ varies in this interval.

5.3.   Controlling chaos

Next, we will conduct numerical simulation of chaos control. Parameter values are fixed as $r = 2.2, \; K = 6, \; k = 10, \; b = 0.2, \; a = 4$, and the initial value is $(x_{0}, y_{0}) = (2.5, 2)$. In Figure 5(f) when the bifurcation parameter $s = 2.8$, system (1.1) will produce chaos. When the feedback gains are $r_{1} = 2$ and $r_{2} = -0.8$, Figure 7(a), (b) show that a chaotic trajectory is stabilized at the fixed point (2.0443, 2.0443). In Figures 8 and 9, when the parameters $r_{1}$ and $r_{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 fixed point and become an asymptotically stable state.

6.   Conclusions

Research shows that the dynamic behavior of discrete systems is richer and more complex than that of continuous systems ^[9,10]. Therefore, based on the previous research work, this paper studies the dynamic behavior and nonlinear characteristics of a class of discrete predator-prey systems with the fear effect. Based on the findings of the research, we can obtain the following mathematical and ecological results:

(1) System (1.1) has two fixed points, in which the only stable fixed point is positive, which reflects the stable coexistence of predators and prey.

(2) System (1.1) has flip bifurcation at the boundary equilibrium point, and flip bifurcation and Neimark-Sacker bifurcation at the positive fixed point. It can be found from Figures 1, 2 and 4 that when $k = 0$, the Neimark-Sacker bifurcation at the positive equilibrium point will produce chaos, and the flip bifurcation will also produce chaos. We can also find the orbits of periods 2, 4, and 8 periodic windows.

(3) When fear $k$ is larger, the number of both predators and prey decreases. It is worth noting that the cost of fear cannot induce the extinction of predators but the extinction of prey. And the system will change from Neimark-Sacker bifurcation to flip bifurcation when $k$ increases (see Figures 1 and 5). Therefore, fear $k$ is very important in analyzing the change of population size.

Acknowledgments

This work was supported by the National Natural Science Foundation of China(Grant Nos. 12161079, 12171192) and Natural Science Foundation of Gansu Province(No.20JR10RA086).

Conflict of interest

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