Stability and Bifurcation Analysis of a Discrete Predator-Prey Model with Mixed Holling Interaction

: In this paper, a discrete Lotka-Volterra predator-prey model is proposed that considers mixed functional responses of Holling types I and III. The equilibrium points of the model are obtained, and their stability is tested. The dynamical behavior of this model is studied according to the change of the control parameters. We ﬁnd that the complex dynamical behavior extends from a stable state to chaotic attractors. Finally, the analytical results are clariﬁed by some numerical simulations. by conducting some numerical simulations to conﬁrm our theoretical results. The model displays fascinating dynamical behavior, including stable behavior, bifurcation, phase portraits, and chaos. It shows that the equilibrium points of the model


Introduction
Mathematical biology tries to model, study, analyze, and interpret biological phenomena, such as the interactions among individuals of different species, which can be predatory, competitive, or mutual. Some simple mathematical models have been suggested to understand the prey-predator interaction. The first attempt in this field was accomplished independently by Lotka [Lotka (1925)] and Volterra [Volterra (1926)], and it is known as the Lotka-Volterra model. This consists of two ordinary differential equations. Many researchers became interested in these models due to their complex behavior [Elettreby and El-Metwally (2007); Elettreby (2009) ;Danca, Codreanu and Bako (1997)]. Many natural factors, such as viz, delay, daily, seasonal effects, and regularity effects, cause complex behaviors in these models [Albert (2012); Gakkhar, Singh and Singh (2012); Hu, Li and Yan (2009)]. Delays due to maturation time of prey, gestation time of the predator, and hunting can induce periodic solutions as well as chaos in the prey-predator model. Many research papers have studied fractional analysis and its applications in prey-predator phenomena  ;].
Holling introduced three kinds of functional responses for different species to model the phenomena of predation [Holling (1965)].
He et al. [He and Li (2014)] investigated the dynamics of a discrete prey-predator model with the functional response of Holling type III in the closed first quadrant R 2 + . Their model consists of two difference equations. The prey-predator interaction term is presented by Holling type III in both equations. They showed that the system undergoes a flip bifurcation and a Neimark-Sacker bifurcation. Saratchandran et al. [Saratchandran, Ajithprasad and Harikrishnan (2015)] proposed a discrete Lotka-Volterra model for prey-predator interactions. The interaction term in the two equations is of Holling type I. They showed how the condition for the survival of the predator depends on the natural death rate of the predator and the efficiency of predation. Their model supported different dynamical regimes asymptotically, including predator extinction, stable fixed p oint, l imit cycle a ttractors for coexistence of predator and prey, and more complex dynamics involving chaotic attractors.
In this paper, we formulate a discrete Lotka-Volterra prey-predator model with a mixed functional response of Holling types I and III in R 2 + . First, we discuss the existence and stability of the equilibrium points of the system. Second, we study the effects on the long-term behavior of the model of changing each control parameter. Finally, we show some numerical simulations to illustrate and support the analytical results.
The remainder of this paper is organized as follows. In Section 2, we introduce a discrete prey-predator model with mixed functional response of Holling types I and III. We study the stability of its equilibrium solutions in Section 3. In Section 4, we study the effects of changing the parameter values and their corresponding regimes. Also, we give different dynamical behaviors, including bifurcation diagrams and different types of attractors. Section 5 discusses our conclusions.

The model
Let y n , x n represent the population of predator species (tiger) density and prey species (zebra) density, respectively, at a discrete time step n. Suppose that the evolution of the prey population follows the logistic map x n+1 = ax n (1 − x n ), where x n ∈ (0, 1) and the positive parameter a represents the constant intrinsic growth rate which is restricted to the interval (0, 4) [May (1976)].
The discrete form of the Lotka-Volterra model can be rewritten as where c (0 < c < 1) is the natural death rate of the predator population. The term p(x n ) y n represents the prey-predator interaction that decreases the prey population, p(x n ) measures the efficiency of predation, and q(x n ) y n represents the prey-predator interaction that increases the predator population. The function q(x n ) is the measure of how effectively the predator can exploit the advantage of predation to increase its population [Shonkwiler and Herod (2009)]. Here, we suppose that q (x n ) = (1 + λ b x n ), where the positive parameter λ is the evolution rate of the survival predator. The parameter λ measures how effectively the predators develop their predation skills to the new generation for their survival. The positive parameter b represents the efficiency of predation (the fraction of the predator population y n that can prey is denoted by b < 1). Thus the predator population increases by the factor b λ x n in each generation.
We suppose that p(x n ) is a Holling type III functional saturation response of a predator or sigmoidal function of the form [Huang, Ruan and Song (2014) where d and e are positive constants. The parameter d represents the predation of the predator coefficient, which reflects the predator's ability, and e is the half saturated. We get the following discrete predator-prey model with the mixed functional response of Holling types I and III: x n+1 = a x n (1 − x n ) − d x 2 n y n e+x 2 n y n+1 = (1 + b λ x n ) y n − c y n The main idea in this model is the use of a mixed functional response of Holling types I and III.

Equilibrium points and their stability analysis
It is challenging to solve system (2) analytically, so we will use qualitative analysis to study it. We will try to find its equilibrium points and study their stabilities. We consider the following Lotka-Volterra prey-predator system with Holling mixed functional responses of types I and III: Proposition 3.1. The equilibrium points of model (3) are Proof. To calculate the equilibrium points (x, y) of system (3), we let x n+1 = x n =x, and y n+1 = y n =ȳ. Using Eq.
Then the growth rate a should be greater than 1 (a > 1). Also, the existence of the third equilibrium point (x 3 ,ȳ 3 ) implies thatȳ 3 should be positive. Then the existence condition of the third equilibrium point is a > 1 and λ > a c b (a−1) . The topoligical nature of the equilibrium point can be a sink, source, or saddle, or it can be non-hyperbolic, according to the eigenvalues for the Jacobian matrix evaluated at an equilibrium point (x,ȳ) [Albert (2012)]. To study the stability of the above equilibrium points, we must compute the coefficients a i j of the Jacobian matrix J [Muller and Kuttler (2015)]. Using system (3), we get The stability of the equilibrium points can be established by calculating the eigenvalues m of the Jacobian matrix J corresponding to each equilibrium point using the characteristic equation |J − m I| = 0, where I is the identity matrix and m is the eigenvalue [Muller and Kuttler (2015); Chou and Friedman (2016)].
Proof. The characteristic equation of the equilibrium point E 1 is Thus the eigenvalues are m 1 = a and m 2 = 1 − c. Then the trivial equilibrium point E 1 = (0, 0) is locally asymptotically stable if each eigenvalue satisfies |m i | < 1, i = 1, 2.
|m 1 | = |a| < 1 then − 1 < a < 1, but we assume that a is positive, hence, the first condition becomes 0 < a < 1. Also, Since c is the death rate of the predator, this leads to 0 < c < 1. Then the trivial steady state E 1 = (0,0) is stable if 0 < a < 1 and 0 < c < 1. The two conditions 0 < a < 1, 0 < c < 1 make the trivial equilibrium point E 1 = (0,0) stable, which means that both prey and predator vanish. This case is biologically unimportant.
Proposition 3.3. The second equilibrium point E 2 = (x 2 ,ȳ 2 ) = a−1 a , 0 of system (3) is stable if the following two inequalities are satisfied: Proof. For the equilibrium point E 2 = a−1 a , 0 , (a > 1), the characteristic equation of the second equilibrium point has the form which can be written as Thus the eigenvalues of the characteristic equation are Then the second equilibrium point |m 1 | = |2 − a| < 1 then 1 < a < 3.
which gives the following two conditions: , since the death rate parameter c must be less than 1 (c < 1), so the first condition a (c−2) b (a−1) < λ is always satisfied. Therefore, the equilibrium point (x 2 ,ȳ 2 ) is stable if the following two conditions are satisfied: Note that the stability of the second equilibrium point (x 2 ,ȳ 2 ) implies the instability of the first trivial equilibrium point (x 1 ,ȳ 1 ). Also, it means the presence of the prey and extinction of the predator. Proposition 3.4. The third equilibrium point E 3 of system (3) is stable if the following condition is satisfied: Proof. The existence of this equilibrium point implies that λ > a c b (a−1) . The characteristic equation of the third equilibrium point which can be written as Then the equilibrium point (x 3 ,ȳ 3 ) is stable, using the Jury test [Keshet (2005) which means that both eigenvalues |m 1,2 | < 1. This condition can be split into two conditions: Let us discuss each condition separately. With some calculations, the first condition, γ < 1, leads to The left side of the inequality of the second condition, −1 − γ < β , leads to . Also, the right side of the inequality of the second condition, β < 1 + γ, leads to which is always satisfied.
From the above discusion, we conclude that the third steady state ( 3.1 Bifurcation analysis of E 3 Since the interior equilibrium point implies the coexistence of prey and predator, this situation is more real and vital than the other two equilibrium points. Thus, we discuss the bifurcation analysis of the equilibrium point E 3 of system (3).
Definition 3.1. The Neimark-Sacker bifurcation (NSB) is the principal tool to determine the existence of quasi-periodic orbits for a map [Asheghi (2018)]. NSB is analogous to the Hopf bifurcation that occurs in a continuous system. In the circumstance that a stable focus loses its stability as a specific parameter is varied with the consequent birth of quasi-periodicity is called a supercritical NSB. In this case, a stable focus enclosed by an unstable closed curve loses its stability, with the consequent disappearance of the closed invariant curve as a parameter is varied. Definition 3.2. The flip bifurcation (FB) occurs when a new limit cycle emerges from a current limit cycle, and the period of the new limit cycle is twice that of the old one. Sometimes flip bifurcation is called period-doubling bifurcation. Definition 3.3. The fold bifurcation is a collision or disappearance of two equilibria in the system. Fold bifurcation is also called saddle-node bifurcation. Lemma 3.1. The interior equilibrium point E 3 transforms from the stable state to: Proof.
(II) Flip bifurcation occurs when a single eigenvalue equals −1. In the context of the characteristic equation of the associated Jacobian matrix J 3 , the condition of the flip bifurcation is expressed as 1 + β + γ = 0. This gives Therefore, the equilibrium point E 3 loses its stability through flip bifurcation when Note that a fold bifurcation or transcritical bifurcation occurs when a real eigenvalue passes through 1, or it is defined as 1 − β + γ = 0. Substituting by β and γ, we get which is not possible because it contradicts the condition of the third equilibrium point E 3 . Then there is no fold bifurcation for this system.

Some numerical results
In Fig. 1, we plot the generation n versus the populations x n of the prey and y n of the predator to check their qualitative behavior. In Figs. 1(a) and 1(b), all of the curves of the populations of the prey x n and predator y n decrease to zero as n increases (which means that both prey and predator vanish) when the conditions 0 < a < 1 and 0 < c < 1 are satisfied. This means that the first steady state (0, 0) is stable. In Figs. 1(a) and 1(b), we change the parameters a, c and fix the other parameters at b = 0.6, d = 0.4, e = 3.1, and λ = 5. This means that the steady state (0, 0) will be stable if the conditions 0 < a < 1 and 0 < c < 1 are satisfied, despite the values of the other parameters. In Fig. 1(c), where a = 1.1 > 1, c = 0.4 < 1, we find that the curve of the prey x n converges to the value a−1 a = 0.09, while the curve of the predator goes to zero. In Fig. 1(d), where a = 1.7 > 1, c = 0.1 < 1, the predator's y n curve goes to infinity, which is not acceptable. This is because the parameters a, c do not satisfy the stability conditions.  Fig. 2, we plot the generation n versus the populations x n of the prey and y n of the predator to investigate their qualitative behavior. In Figs. 2(a) and 2(b), all of the curves of the population x n of the prey converge to the value a−1 a , while the curves of the population y n of the predator tend to zero. This is due to satisfying the stability conditions of the second steady state, 1 < a < 3, λ < a c b (a−1) . In Figs. 2(c) and 2(d), the parameter values lead to instability of the second steady state (x 1 ,ȳ 1 ). The population curves of the prey and predator move away from this steady state. In Fig. 2(c), the population curves of the prey and predator go to zero. This is because the stability conditions of the second steady state are not satisfied (a < 1). Also, in Fig. 2(d), the predator population curve moves away from zero and goes to infinity (λ = 3.1 > a c b (a−1) = 0.254). Note that when the second steady state (x 1 ,ȳ 1 ) is stable, the trivial steady state will be unstable, and when the trivial steady state is stable, then the second steady state (x 1 ,ȳ 1 ) is unstable.
In Fig. 3, we plot the generation n versus the populations x n , y n to test their qualitative behavior. In Figs. 3 Figure 3: Qualitative behavior of the populations x n of the prey and y n of the predator that show the stability of the third steady state to determine the zone of stability located between the upper (bifurcation) and lower (extinction) curve. We note that increasing the growth parameter a leads to a smaller stability area. In Fig. 4(a), we fix c = 0 .5 a nd e = 3 .1 a nd p lot r egions f or t hree values of b(= 0.1, 0.2, 0.3). From the figure, w e n ote t hat i ncreasing t he value o f b l eads t o a decrease in the region of stability and shifts the regions downward. In Fig. 4(b), we fix the parameters b = 0.6, e = 3.1 and plot regions for three values of parameter c(= 0.5, 0.6, 0.7). We observe that increasing the value of c leads to an increase in the region of stability and shifts the regions upward. Also,In Figs. 4(c) and 4(d), we plot parameter b, the fraction of predators that can prey, vs. the parameter λ to determine the region of stability. We fix t he p arameters a = 2 .1, e = 3 .1 a nd p lot r egions f or t hree v alues o f parameter c(= 0.01, 0.03, 0.05). We observe that increasing the value of c leads to an increase in the region of stability and shifts the regions upward. Also, we note that increasing the value of b leads to a decrease in the region of stability and shrinks the regions until they vanish. In Fig. 4(d), we increase the value of a from 2.1 to 2.4 and plot regions corresponding to three values of parameter e(= 1.1, 2.1, 3.1). We observe that increasing the value of c leads to a decrease in the region of stability. Also, increasing the value of b shrinks the stability regions until they vanish.
It is clear from these figures t hat c hanging t he p arameters a, b, c ,λ l eads t o a r ange of dynamical behaviors. We will numerically study the complex dynamical behaviors of our system. We will use two sets of parameter values of the system. For each case, we will study the bifurcations and the nature of the attractors by independently changing the In this case, we use a = 2.4, d = 0.4, and e = 3.1. In Fig. 5, we plot the respective populations x n , y n of the prey and predator versus the parameter λ . We start the graph at λ = 1 to focus on its important parts. Note that at λ = 1, we get (x,ȳ) = (0.583, 0) = (x 2 ,ȳ 2 ). Increasing λ leads to a decrease of the prey population and an increase of the predator population. This is consistent with the above analytical results. Then the second steady state (x 2 ,ȳ 2 ) becomes unstable and the third steady state appears. The third steady state (x 3 ,ȳ 3 ) is stable until λ reaches the threshold value λ th . After that, the bifurcation structure for the prey and predator is as shown in the figure. Also, there are many periodic windows within the chaotic regime. As λ increases, we see a transition phase from stability to bifurcation within a limit cycle, to a periodic window, and finally to chaos.
In Fig. 6, we choose the parameter values d = 0.4, e = 3.1, and λ = 1.8. We plot the populations x n , y n of the prey and predator, respectively, versus the parameter a in the interval (1, 4). A bifurcation structure for the populations of the prey and predator appears when the parameter a reaches a critical value. As a increases, we note a transition from stability to bifurcation within a simple limit cycle, periodic window, and finally to chaos. Case II: b = 0.8, c = 0.6 In this case, we use the values a = 2.8, d = 0.4, and e = 2.1. In Fig. 7, we plot the populations x n , y n for the prey and predator, respectively, vs. λ to investigate the dynamical behavior of our system. We start the graph at λ = 1 to focus on its important parts. Note that at λ = 1, we have (x,ȳ) = (0.6428, 0) = (x 2 ,ȳ 2 ). Increasing λ leads to a decrease of the prey population and an increase of the predator population. This is consistent with the above analytical results. Then the second steady state (x 2 ,ȳ 2 ) becomes unstable and the third steady state appears. The third steady state (x 3 ,ȳ 3 ) becomes stable until λ reaches a threshold value λ th . After that, the bifurcation structure for the prey and predator is as shown in the figure. Also, there are many periodic windows within the chaotic regime. As λ increases, we see a transition phase from stability to bifurcation within a limit cycle, to a periodic window, and finally to chaos.

Discussion and conclusion
In this paper, we have proposed and investigated the complex behavior of a discrete prey-predator model with the mixed functional response of Holling types I and III. The main idea is to study the effects of changing model parameters on the dynamics of the model. The equilibrium points of the model are obtained. There are two equilibrium boundary points and a unique interior equilibrium point. The conditions of the local stability of these points have been proved using the Jury criterion. Also, we study the bifurcation analysis of the interior equilibrium point. We have followed the classical work of Lotka-Volterra by conducting some numerical simulations to confirm our theoretical results. The model displays fascinating dynamical behavior, including stable behavior, bifurcation, phase portraits, and chaos. It shows that the equilibrium points of the model lose stability via bifurcation. These results reveal rich dynamics of the discrete-time models.