Bifurcation and complex dynamics of a discrete-time predator-prey system with simplified Monod-Haldane functional response

*Correspondence: srana.mthdu@gmail.com Department of Mathematics, University of Dhaka, Dhaka, 1000, Bangladesh Abstract In this paper, we investigate the dynamics of a discrete-time predator-prey system with simplified Monod-Haldane functional response. The existence and local stability of positive fixed point of the discrete dynamical system is analyzed algebraically. It is shown that the system undergoes a flip bifurcation and a Neimark-Sacker bifurcation in the interior of R+ by using bifurcation theory. Numerical simulation results not only show the consistence with the theoretical analysis but also display new and interesting dynamical behaviors, including phase portraits, period-11 orbits, attracting invariant circle, cascade of period-doubling bifurcation from period-11 leading to chaos, quasi-periodic orbits, and the sudden disappearance of the chaotic dynamics and attracting chaotic set. The Lyapunov exponents are numerically computed to characterize the complexity of the dynamical behaviors.


Introduction
The dynamics of predator-prey interaction is the starting point for many variations that yield more realistic biological and mathematical problems in population ecology. Predation is a direct interaction which occurs when individuals from one population derive their nourishment by capturing and ingesting individuals from another population. There are many articles devoted to the study of predator-prey interaction both from the experimental and the modeling point of view. It is well known the Lotka-Voltera predator-prey model is one of the fundamental population models; a predator-prey interaction has been described firstly by two pioneers Lotka (1924) and Voltera (1926) in two independent works. After them, more realistic prey-predator model IAEES www.iaees.org were introduced by Holling suggesting three types of functional responses for different species to model the phenomena of predation (Holling, 1965).
Qualitative analyses of prey-predator models describe by set of differential equations were studied by many authors (Brauer and Castillo, 2001;Hastings and Powell, 1991;Klebanoff and Hastings, 1994;May, 1974;Murray, 1998;Zhu et al., 2002). Another possible way to understand a prey-predator interaction is by using discrete-time models. These models are more reasonable than the continuous time models when populations have non-overlapping generations (Brauer and Castillo, 2001;Murray, 1998) and lead to unpredictable dynamic behaviors from a biological point of view. This suggests the possibility that the governing laws of ecological systems may be relatively simple and therefore discoverable. The author (May, 1975(May, , 1976 had clearly documented the rich array of dynamic behavior possible in simple discrete-time models. Recently, there is a growing evidence showing that the dynamics of the discrete-time prey-predator models can present a much richer set of patterns than those observed in continuous-time models (Agiza et al., 2009;Danca et al., 1997;Elsadany et al., 2012;Hasan et al., 2012;He and Lai, 2011;Jing and Yang, 2006;Li, 1975;Liu, 2007;Hu et al., 2011;He and Li, 2014). However, there are few articles discussing the dynamical behaviors of predator-prey models, which include bifurcations and chaos phenomena for the discrete-time models. The authors (He and Lai, 2011;Jing, 2006;Liu, 2007;Hu et al., 2011) obtained the flip bifurcation by using the center manifold theorem and bifurcation theory. But in (Agiza et al., 2009;Danca et al., 1997;Elsadany et al., 2012), the authors only showed the flip bifurcation and Hopf bifurcation by using numerical simulations. In this work, we confine our interest to present, by using both analytic and numerical methods, the domains of the values of the parameters under which the system predicts that the populations will be able to persist at a steady state, the conditions for flip and/or Neimark-Sacker bifurcations by using the normal form theory of the discrete system (see section 4, Kuznetsov, 1998) and the domain for the presence of chaos in the system by measuring the maximum Lyapunov exponents.
In ecology, many species have no overlap between successive generations, and thus their population evolves in discrete-time steps (Murray, 1998). Such a population dynamics is described by difference equation.
Let n x denotes the number of prey population and n y the number of predator population in the n th generation. Our model is described by the following system of nonlinear difference equations in non- (1) In the system (1), the prey grows logistically with intrinsic growth rate r and carrying capacity one in the absence of predation. The predator consumes the prey with functional response Holling type I. All have positive values that stand for prey intrinsic growth rate, per capita searching efficiency of the predator, conversion rate, and the death rate of the predator, respectively. From mathematical and biological point of view, we will pay attention on the dynamical behaviors of (1) in the closed first quadrant 2  R . Starting with initial population size   0 0 , y x , the iteration of system (1) is uniquely determined a trajectory of the states of population output in the following form Our results in this paper are extension to those in (Danca et al., 1997;Elsadany et al., 2012). This paper is organized as follows. In Section 2, we discuss the existence and local stability of positive fixed point for system (1) in 2  R . In Section 3, we show that there exist some values of the parameters such that (1) undergoes the flip bifurcation and the Neimark-Sacker bifurcation in the interior of 2  R . In section 4, we present the numerical simulations which not only illustrate our results with theoretical analysis but also exhibit complex dynamical behaviors such as the cascade periodic-doubling bifurcation in periods 2, 4, 8, 9, 10, 20orbits, quasi-periodic orbits and chaotic sets. Finally a short discussion is given in Section 5.

Existence and Local Stability of Fixed Points
In this section, we shall first discuss the existence of fixed points for (1), then study the stability of the fixed point by the eigenvalues for the Jacobian matrix of (1) at the fixed point.It is clear that the system (1) has the following fixed points in the ) , ( y x -plane: To discuss the existence of fixed points, we say that fixed points will not exist if any one of its components is negative. The fixed point 0 E always exists. The existence condition for 1 E is 1  r . Finally, the feasibility condition for the positive fixed point 2 ). The Jacobian matrix due to the linearization of (1) evaluated at 2 E is given by r d y x J and the characteristic equation of the Jacobian matrix J can be written as Using Jury's criterion (Elaydi, 1996), we have necessary and sufficient condition for local stability of the fixed point 2 E which are given in the following proposition.
, then system (1) has a positive fixed point 2 E and (i) it is a sink if one of the following conditions holds: (ii) it is a source if one of the following conditions holds: (iii) it is non-hyperbolic if one of the following conditions holds: (iv) it is a saddle for the other values of parameters except those values in (i)-(iii).
Following Jury's criterion, we can see that one of the eigenvalues of   2 E J is 1  and the others are neither 1 nor 1  if the term (iii.1) of Proposition 1 holds. Therefore, there may be flip bifurcation of the fixed point 2 Also when the term (iii.2) of Proposition 1 holds, we can obtain that the eigenvalues of   2 E J are a pair of conjugate complex numbers with module one. The conditions in the term (iii.2) of Proposition 1 can be written as the following set: and if the parameter r varies in the small neighborhood of 2 E NS ; then the Neimark-Sacker bifurcation will appear.

Flip Bifurcation and Neimark-Sacker Bifurcation
In this section, we choose the parameter r as a bifurcation parameter to study the flip bifurcation and the Neimark-Sacker bifurcation of 2 E by using bifurcation theory in (see Section 4 in Kuznetsov, 1998; see also Guckenheimer and Holmes, 1983;Robinson, 1999;Wiggins, 2003).
We first discuss the flip bifurcation of (1) at 2 (1) into the origin, then system (1) becomes In order to normalize p with respect to q , we denote means the standard scalar product in R 2 : Following the algorithms given in (Kuznetsov, 1998), the sign of the critical normal form coefficient   1 1 r  , which determines the direction of the flip bifurcation, is given by the following formula: From the above analysis and the theorem in (Kuznetsov, 1998;Guckenheimer and Holmes, 1983;Robinson, 1999;Wiggins, 2003), we have the following result. In Section 4, we will give some values of the parameters such that   0 1 1  r  , thus the flip bifurcation occurs as r varies (see Figure 1). We next discuss the existence of a Neimark-Sacker bifurcation by using the Neimark-Sacker theorem in (Kuznetsov, 1998;Guckenheimer and Holmes, 1983;Robinson, 1999;Wiggins, 2003 For 2 r r  , the eigenvalues of the matrix associated with the linearization of the map (6) For the above argument and the theorem in (Kuznetsov, 1998;Guckenheimer and Holmes, 1983;Robinson, 1999;Wiggins, 2003), we have the following result. In Section 4 we will choose some values of the parameters so as to show the process of a Neimark-Sacker IAEES www.iaees.org bifurcation for system (1) in Figure 2 by numerical simulation.

Numerical Simulations
In this section, our aim is to present numerical simulations to explain the above theoretical analysis, especially the bifurcation diagrams, phase portraits and Maximum Lyapunov exponents for system (1)   For case (i). The bifurcation diagrams of system (1) in ( y x r   ) space and in ( x r  ) pane are given in Fig. 1(a-b) , which determines the direction of the flip bifurcation and shows the correctness of Theorem1.From Fig. 1(b), we see that the fixed point 2 E is stable for 3 . 3  r and loses its stability at the flip bifurcation parameter value 3 . 3  r , we also observe that there is a cascade of period doubling bifurcations for 3 . 3  r . The maximum Lyapunov exponents corresponding to Fig. 1(b) are computed and plotted in Fig. 1(c), confirming the existence of the chaotic regions and period orbits in the parametric space.
For case (ii). The bifurcation diagrams of system (1) in the ( y x r   ) space, the ( x r  ) plane and the ( y r  ) plane are given in Fig. 2(a-b-c) Fig. 2(b-c), we observe that the fixed point 2 E of map (1) is stable for 2  r and loses its stability at 2  r and an invariant circle appears when the parameter r exceeds 2 , we also observe that there are period-doubling phenomenons. The maximum Lyapunov exponents corresponding to Fig. 2(b-c) are computed and plotted in Fig. 2(d), confirming the existence of the chaotic regions and period orbits in the parametric space. From Fig. 2(d), we observe that some Lyapunov exponents are bigger than 0, some are smaller than 0, so there exist stable fixed points or stable period windows in the chaotic region. In general the positive Lyapunov exponent is considered to be one of the characteristics implying the existence of chaos. The bifurcation diagrams for x and y together with maximum Lyapunov exponents is presented in Fig. 2(e). Fig. 2

Fig. 2(f) is the local amplification corresponding to
The phase portraits which are associated with Fig. 2(a) are disposed in Fig. 3  . In Fig. 4