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

: In this paper, we investigate the complex dynamics of a classical discrete-time prey-predator system with the cost of anti-predator behaviors. We ﬁrst give the existence and stability of ﬁxed points of this system. And by using the central manifold theorem and bifurcation theory, we prove that the system will experience ﬂip bifurcation and Neimark-Sacker bifurcation at the equilibrium points. Furthermore, we illustrate the bifurcation phenomenon and chaos characteristics via numerical simulations. The results may enrich the dynamics of the prey-predator systems.


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) = 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 x t+1 = x t + rx t 1+ky t (1 − x t K ) − bx t y t x t +a , y t+1 = y t + sy t (1 − y t x t ), (1.1) 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) = 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.

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 The Jacobian matrix J(x, y) corresponding to system (1.1) at point (x, y) is as follows: (2.2) Assume that λ 1 and λ 2 are two roots of the characteristic equation of the Jacobian matrix J| (x,y) , and we have the following definition and conclusions.

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 * ).

Bifurcation at E 0 (K, 0)
The eigenvalues of J(E 0 ) are 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 * | 1 is a small perturbation parameter. where 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 0, h 2 0, 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).

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) ∈ F.B.2 and consider the following system: where s is the bifurcation parameter, |s * | 1 is a small perturbation parameter. where x * .We construct a nonsingular matrix and use the translation .9) 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 1 (0) is given by where Remember If parameters β 1 and β 1 are not 0, system (3.13) has flip bifurcation.And we can get the following theorem: Theorem 3.2.If β 1 0, β 2 0, 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 β 2 > 0 (respectively, β 2 < 0), system (1.1) bifurcates from the fixed point E * (x * , y * ) to a 2-periodic stable orbit(respectively, unstable).
Therefore, eigenvalues λ 1 , λ 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.
, 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 L 0, 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 L 0 in (3.17).Moreover, if L < 0(resp., L > 0), an attracting (resp., repelling) invariant closed curve bifurcates from the steady state for s > s * (resp., s < s * ).

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 Thus, the characteristic equation corresponding to J (x * , y * ) is: Let λ 1 and λ 2 be the eigenvalues of characteristic equation (4.3), then and In order to make the absolute values of λ 1 and λ 2 less than 1, we assume that λ 1 = ±1 and λ 1 λ 2 = 1 hold.
Assume that λ 1 λ 2 = 1, and we have Assume that λ 1 = 1, and we get Assume that λ 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.

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.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).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 , • • •.
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.

Mathematical Biosciences and Engineering
Volume 19, Issue 7, 6659-6679.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

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.

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.

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 r = 2.2, b = 0.2, a = 4, K = 6, (5.1) (a) Bifurcation diagram for x t .(b) Bifurcation diagram for y t .

Figure 1 .
Figure 1.The bifurcation diagram of system(1.1)corresponding to the bifurcation parameter s.

Figure 5 .
Figure 5.The bifurcation diagrams and maximum Lyapunov exponents of system (1.1) for different values of k.
(a) Bifurcation diagram for x t .(b) Bifurcation diagram for y t .

Figure 7 .
Figure 7.The bifurcation diagrams and maximum Lyapunov exponents of system (1.1) corresponding to the bifurcation parameter k.
* (x * , y * ) is a source point and it is locally unstable if and only if p < |1 + q| and q > 1;