Bifurcations and dynamics of a discrete predator–prey model of ricker type

A discrete-time predator–prey model is investigated in this paper. In considered model, the population is assumed to follow the model suggested by Ricker 1954. Existence and stability of equilibria are studied. Numerical simulations reveal that, depending on the parameters, the system has complicated and rich dynamics and can exhibit complex patterns. Also the bifurcation diagrams are presented.

where X stands for the number of prey and Y stands for the number of predators. Furthermore, the parameters a, b, c, and d are all positive. It has been demonstrated that in the absence of predators, the number of prey increases exponentially, while the number of predators declines exponentially in the absence of a prey population. The terms bXY and d XY describe predator-prey interactions that are beneficial to predators but deadly to prey. A good update for this model is to replace the term a X by the logistic growth a X(1 − X ) where population's per capita growth rate decreases as population size approaches a maximum imposed by limited resources. Stability, bifurcations and chaos control of the improved system studied by [37]. However, we are quite interested in studying the effect of using Ricker map r X n e 1− Xn k instead of the logistic map to describe the dynamics of preys. In Ricker map, the population growth is virtually exponential; but, as population size grows, the instantaneous growth rate decreases linearly; and finally, population size plateaus and fluctuates around a mean.
The outline of this paper is as follows: The construction of the model is discussed in Sect. 2. Section 3 investigates the existence of equilibria and local stability of the model (5) with various topological types. We explore the existence of bifurcations around equilibria in Sect. 4. Numerical simulations using MATLAB are applied in Sect. 5 to support the theoretical analysis and visualize the newly observed complex dynamics of the system. Section 6 contains the conclusion and discussion.

Derivation of the model
This paper deals with the study of stability and bifurcations of a discrete-time predatorprey model. Let X n and Y n represent the population density of prey and predator at time n, respectively. We impose the following assumptions to formulate the difference equations which describe the required system.

Assumption 1
In the absence of predator Y , the population of prey X would grow in a classical discrete Ricker population map manner.
where r > 0 is per-capita birth or growth rate in absence of predators r is the prey finite rate of increase in the absence of predators and k > 0 is the carrying capacity parameter.

Assumption 2
The population of prey X will drop in the presence of the predator Y . If we suppose that the interaction between the prey and the predator has the mass action principle response due to the prey's defensive ability, then Eq. (1) can be expressed as where b > 0 is the biomass conversion parameter.
Assumption 3 Prey X are favourite food for predator Y , the population density of predator Y will increase following the mass action principle response manner in the presence of favourite food. Hence we have where d > 0 denotes the predator's biomass conversion efficiency.
Hence the proposed model is from Eqs. (1)-(3) with positive initial condition X 0 > 0 and Y 0 > 0. Let X n = k x n and Y n = y n b . We rewrite model Eq. (4) as x n+1 = r x n e 1−x n − x n y n , where α = d k > 0.

Fixed points: existence and stability
In this section, we study the existence and stability of fixed points of the system (5)  Proof The fixed points correspond to the steady states of the model population (5) can be obtained by solving the algebraic system then all non-negative fixed points of the system (5) are E 0 (0, 0), The Jacobian matrix J (x, y) of the model (5) about the fixed point (x, y) is given by Its characteristic equation is The biological interruption for three fixed points is as follows: The fixed point E 0 represents a situation in which there is no prey and no predators. The fixed point E 1 reflects the case where there is prey but no predators. E 2 refers to the coexistence of fixed nonzero numbers of prey and predators.

Dynamic behavior of the model
Now we will study the behavior of the solutions of (5) about E 0 , E 1 and E 2 . To do this we recall the following definitions: Definition 1 Let E(x, y) be a fixed point of (5) and let μ 1 and μ 2 are the eigenvalues of (7).
The relations between roots and coefficients of the quadratic equation is given by the following Lemma.  The following Proposition confirms the stability of fixed points of system (5) under some conditions.

Proposition 2
For the fixed point E 0 of the system (5), the following statements are true: Proof The Jacobian matrix associated with (5) at E 0 is J (E 0 ) and is given by

Proposition 3
The fixed point E 0 of the system (5) is globally asymptotically stable if 0 < r < e −1 .
Proof It was shown in Proposition 2 that E 0 is locally asymptotically stable if 0 < r < e −1 . So, proving that lim n→∞ x n = 0 and lim n→∞ y n = 0 is sufficient to prove the global stability of E 0 . It follows from the first equation of system (5) that Then Then lim n→∞ x n = 0.
Note that ⇐⇒ r e e −x n ≥ y n ⇐⇒ r e ≥ y n e x n .
Since x < e x for all x ≥ 0, then x n y n < y n e x n < r e. Then it follows from second equation of system (5) that y n+1 ≤ αr e. That is y n is bounded from above. Again we see from the second equation of system (5) and Eq. (8) that y n+1 = αx n y n ≤ α(r e) n x 0 y n , n = 1, 2, ...
Thus y n ≤ (r e) n α n x n 0 y 0 .
Since y n is bounded and r e < 1 therefore lim n→∞ y n = 0.
Then the proof is completed. (5), the following statements are true:

Proposition 4 For the fixed point E 1 of the system
Proof In the same manner as in Proposition (2). The Jacobian matrix associated with (5) at E 1 is J (E 1 ) and is given by hence, the eigenvalues of J (E 1 ) are λ 1 = − ln(r ) and λ 2 = α (1 + ln(r )). Then, by using Lemma (1) we get the results of Proposition (4).
We are interested to analyze the dynamical behavior of the system around interior fixed point E 2 for its importance. In order to determine conditions that guarantee that the characteristic roots of Jacobian matrix about E 2 are inside the unit disk one can use Lemma (1). (5) is sink if and only if the following conditions hold:

Proposition 5 The positive equilibrium point E 2 of the model
The Jacobian matrix (6) at E 2 reads as , and the characteristic equation of J (E 2 ) is according to Lemma (1), the positive equilibrium point E 2 is locally asymptotically stable if The nonhyperbolic fixed point case is more involved. According to the eigenvalues of 1 or − 1, there are several scenarios. It is common to use centre manifold theory to determine the stability of the fixed point when one of the eigenvalues is on the unit circle and the other eigenvalue is inside the unit circle [38][39][40]. In order to determine conditions that guarantee that the characteristic roots in the unit circle, one can use Theorem 4.5 (Elaydi [40], p. 203). Proof By Proposition (4) the stability boundary of E 2 consists of parts of three curves, namely, Curve 3: r = e 1−α α , the points of Curve 1 on the stability boundary of E 2 meet 1 + p + q = 0, implying that they have an eigenvalue of −1 and are hence period doubling points. The points on the stability boundary of Curve 2 satisfy q = 1, i.e. they have two eigenvalues with product 1 and are hence Neimark-Sacker points. Curve 3 points on the stability boundary satisfy 1 − p + q = 0, implying that they have an eigenvalue of 1. It is simple to verify that E 2 then branches to E 1 .

Bifurcations analysis
The existence of bifurcations about equilibria is investigated in this section, which is based on theoretical studies in Sect. 3. It is well known that the classification of bifurcation types depends mainly on the natural of the eigenvalues about the equilibrium points. The following is a summary of the presence of related bifurcations around the equilibria E 0 , E 1 , and E 2 : (2), we can find out that when r = e −1 , one of the eigenvalues about the equilibrium E 0 is 1. So when the parameter varies in the small neighborhood of r = e −1 , a fold bifurcation may occur.
(ii) From Proposition (2), we can see that if r = e 1−α α holds then one of the eigenvalues about E 1 is 1 . So a fold bifurcation may occurs when the parameter varies in a small neighborhood of r = e 1−α α . And we denote the parameters satisfying r = e 1−α α as (iii) If (i) holds in Proposition (6), we can observe that one of the eigenvalues of J E 2 about E 2 is −1, while the other is neither 1 nor −1 . As a result, a perioddoubling bifurcation exists due to parameter variation in a small neighborhood of r = 3α 2−α e 1−α α . More precisely we can also represent the parameters satisfying (iv) If (ii) holds in Proposition (6), we see that the eigenvalues of J E 2 about E 2 are a pair of complex conjugate with modulus 1 . So a Neimark-Sacker bifurcation exists by the variation of parameter in a small neighborhood of r = α α−1 e 1−α α .
Precisely we represent the parameters satisfying r = α α−1 e 1−α α as The stability regions (i), (ii) and (iii) of E 2 obtained in Proposition 3 are depicted as S E2 in Fig. 1. S E0 and S E1 in Fig. 1 represent the stability regions of E 0 and E 1 , respectively. The following is the biological interpretation of Fig. 1. For all parameter values, the situation with no prey and no predators exists, but it is only stable in S E0 , that is, for r < 1 e . For every r > 1 e , the situation with a fixed number of prey but no predators exists, but only in S E2 is stable. When r > 1 e and α > 1 1+ln(r ) , coexistence of fixed nonzero numbers of prey and predators is possible, but only in S E2 is stable.
In order for (11) to undergo a Neimark-Sacker bifurcation, it is mandatory that the following discriminatory quantity, i.e., χ = 0 (see [38,40,41]), where τ 02 = 1 8 X n X n − Y n Y n + 2 X n Y n + ι X n X n − Y n Y n + 2 X n Y n (0,0) , τ 11 = 1 4 X n X n + Y n Y n + ι X n X n + Y n Y n (0,0) , τ 20 = 1 8 X n X n − Y n Y n + 2 X n Y n + ι X n X n − Y n Y n − 2 X n Y n (0,0) , τ 21 = 1 16 X n X n X n + X n Y n Y n + X n X n Y n + Y n Y n Y n +ι X n X n X n + X n Y n Y n − X n X n Y n − Y n Y n Y n | (0,0 Based on this analysis and the Neimark-Sacker bifurcation Proposition discussed in [38,[42][43][44][45], we arrive at the following Proposition.

Remark
According to bifurcation theory discussed in [43], the bifurcation is called a supercritical Neimark-Sacker bifurcation if the discriminatory quantity χ < 0. In the following section, numerical simulations guarantee that a supercritical Neimark-Sacker bifurcation occurs for the model (5).

Period-doubling bifurcation
Here we study the period-doubling bifurcation of model (5) at E 2 when parameters vary in a small neighborhood of P E 2 . Select arbitrary parameters (r * , α) from P E 2 . We consider the parameter r * as a new dependent variable, and we can get let u n = x n − x * , v n = y n − y * then the equilibrium E 2 of the discrete-time model (13) transforms into O(0, 0). By calculating we get u n+1 = 11 u n + 12 v n + 13 u 2 n + 14 u n v n + ϒ 01 u n r * + ϒ 02 u 2 n r * + O (|u n | , |r * |) 4 , v n+1 = 21 u n + 22 v n + 23 u n v n , (14) where 11 and use the translation gives where (u n , v n , r * ) = 13 (λ 2 − 11 ) u n r * = 12 X n r * + 12 Y n r * , u 2 n r * = 12 2 X 2 n r * + 2X n Y n r * + Y 2 n r * , hereafter we determine the center manifold W c (0, 0) of (15) about (0, 0) in a small neighborhood of r * [6,38,43,46]. By center manifold theorem, there exists a center manifold W c (0, 0) that can be represented as follows: where O (|X n | + |r * |) 3 is a function with order at least three in their variables (X n , r * ), and c 0 = 0, Therefore, we consider the map (15) restricted to W c (0, 0) as follows: where in order for the map (16) to undergo a period-doubling bifurcation, we require that the following discriminatory quantities are non-zero: , after calculating we obtain from the above analysis in [37] and theorem in [38,[42][43][44][45] ,we have the following Proposition Proposition 7 If 2 = 0, the map (13) undergoes a period-doubling bifurcation about the unique positive equilibrium E 2 when r * varies in a small neighborhood of O(0, 0). Moreover, if 2 > 0 ( resp . 2 < 0), then the period-2 points that bifurcate from E 2 are stable (resp. unstable).

Numerical simulations
We use numerical simulations to create bifurcation diagrams and the phase portraits of model (5) to support the previous theoretical analysis and highlight new complicated dynamical behaviors. We will use the MATCONTM package [43,47] to carry out the numerical continuation of equilibria in order to obtain bifurcation curves of system (5) for different sets of parameters and locate points of interest. The bifurcation parameters are considered in the following four cases. 0.1, 0.1), We see that model (5) has stable equilibrium point E 0 (0, 0), there are no prey or predator in the area. Figure 2 shows the correctness of discussion about equilibrium E 0 = (0, 0) in Sect. 3. From Fig. 2d we see that equilibrium E 0 = (0, 0) loses its stability and bifurcates to E 1 when r = e −1 . The MATCONTM report is label = BP, x = 0.000000 0.000000 0.367879 .

Case 2
Choosing r = 1, α = 0.4, with initial value (x 0 , y 0 ) = (0.5, 0.2), model (5) has stable equilibrium point E 1 (1, 0). It can be observed from Fig. 3a that the prey density x is fixed at 1 but the predator population y quickly reduces to zero. In stability after some calculations from Maple one gets τ 02 = 0.0761055464 + 0.5739851115ι, τ 11 = 0.1651982395 + 0.6481391100ι, τ 20 = 0.0380527732 + 0.03707699925ι, in view of (17) and (18) The numerical solution in Figs. 4 and 5 reveals that the stable fixed point E 2 becomes unstable, allowing preys and predators to coexist by a persistent positive periodic oscillation as time passes. Case 4 Varying r in range 2.3 ≤ r ≤ 3.89 and fixing α = 0.63. Figure 6 shows that prey and predator densities are fixed at the stable equilibrium point E 2 for r < 2.482007420, and E 2 loses its stability when r = 2.482007420. Further, when r > 2.482007420, the numbers of prey and predator oscillate in a stable period-2 orbit then period 4-orbit and so on, that is a period-doubling bifurcation. Moreover, a chaotic set is emerged with the increasing of r .

Conclusion
The dynamical behaviors of model (5) is discussed in this study. Model (5) contains two boundary equilibria, E 0 , E 1 , and a positive equilibria E 2 when r > max{e −1 , e 1−α α }. We explored local stability along topological types about E 0 , E 1 and E 2 by using the linearization method. From the discussion in Sect. 4.2, we know there is period- (c) stable period-4 orbit (d) stable period-8 orbit  Comparing with the discrete-time model (4) in [37] where the author used the logistic growth r X n (1 − X n ) in the first equation of the model. This logistic map has the unrealistic fact that the parabola drops below the horizontal axis as we move off to the right. A more reasonable choice would be the Ricker map r X n e 1−X n which produces very small, but still positive, values as X n becomes very large [2]. The global qualitative analysis of model (5) has not been obtained yet, we will leave it for our future work.