Complex Dynamics on the Routes to Chaos in a Discrete Predator-Prey System with Crowley-Martin Type Functional Response

We present in this paper an investigation on a discrete predator-prey system with Crowley-Martin type functional response to know its complex dynamics on the routes to chaos which are induced by bifurcations. Via application of the center manifold theorem and bifurcation theorems, occurrence conditions for flip bifurcation and Neimark-Sacker bifurcation are determined, respectively. Numerical simulations are performed, on the one hand, verifying the theoretical results and, on the other hand, revealing new interesting dynamical behaviors of the discrete predator-prey system, including period-doubling cascades, period-2, period-3, period-4, period-5, period-6, period-7, period-8, period-9, period-11, period-13, period-15, period-16, period-20, period22, period-24, period-30, and period-34 orbits, invariant cycles, chaotic attractors, sub-flip bifurcation, sub-(inverse) NeimarkSacker bifurcation, chaotic interior crisis, chaotic band, sudden disappearance of chaotic dynamics and abrupt emergence of chaos, and intermittent periodic behaviors. Moreover, three-dimensional bifurcation diagrams are utilized to study the transition between flip bifurcation and Neimark-Sacker bifurcation, and a critical case between the two bifurcations is found. This critical bifurcation case is a combination of flip bifurcation andNeimark-Sacker bifurcation, showing the nonlinear characteristics of both bifurcations.


Introduction
Predator-prey interaction shows widespread existence in nature and can take many forms, such as resource-consumer, plant-herbivore, and phytoplankton-zooplankton forms [1].Due to the ubiquity and importance of the predator-prey interaction between populations, the dynamics of predatorprey systems have attracted the attention of many scholars, and the research on predator-prey systems has become one of the dominant themes in ecology [2][3][4].In recent decades, mathematical models have been established to analyze various complex dynamics of the predator-prey systems in various circumstances [4][5][6].
The Lotka-Volterra predator-prey model, proposed in the pioneering works of Lotka and Volterra [7,8], is probably the simplest and most basic predator-prey model in the field.The development of subsequent predator-prey models is mostly based on the Lotka-Volterra model [3,9,10].The development may be a change of functional response (e.g., Holling type I, II, and III classifications) or numerical response to describe different ecological processes in the predator-prey interaction [11,12].Several researchers argued that the nonautonomous model is closer to the realistic predator-prey system than the autonomous model [3].Therefore, predator-prey models with a time delay were also studied because the predator-prey interaction may have time lag [13].In order to incorporate the influence of seasonal variation into the predator-prey interaction, predator-prey models with periodically varying parameters were also proposed and investigated [14][15][16].
In earlier approaches, most predator-prey models were time-continuous.However, the continuous models hardly described the dynamics when the populations have nonoverlapping generations or the number of populations is small [17,18].For these cases, applying discrete-time models should be more reasonable and appropriate [17].In recent years, more and more scholars have payed attention to discretetime predator-prey models.Compared with the continuous models, the discrete dynamic models often have advantages in describing richer nonlinear characteristics and complexity, such as various periodic orbits, quasiperiodic behaviors, and chaotic attractors [19][20][21].Moreover, a few research works found that the discrete model can lead to more accurate results than the corresponding continuous model in describing predator-prey dynamics [6,22,23].
Through the research of He and Lai [2], Liu and Xiao [19], Jing and Yang [20], Agiza et al. [5], Hu et al. [6], Huang and Zhang [21], and others, they found out that flip bifurcation and Neimark-Sacker bifurcation are two basic bifurcations occurring in discrete predator-prey systems.The two bifurcations often start the routes to chaos, revealing the dynamic transition from the nonchaotic state to the chaotic state.Moreover, complex dynamics on the routes to chaos can be observed, including cascades of period doubling, periodic windows, intermittent chaos, and chaotic crisis [19][20][21].For example, via calculating Lyapunov exponents and fractal dimension, Agiza et al. determined diverse strange attractors of a predator-prey system [5].
When we study the nonlinear characteristics of predatorprey systems, the functional response of the predator to the prey is an essential factor which affects the dynamical properties.From an ecological point of view, functional responses may be determined by the prey escape capability, prey habitat property, and predator hunting capability [24].Skalski and Gilliam [12] suggested that three classical preypredator-dependent functional responses-Beddington-DeAngelis, Crowley-Martin, and Hassell-Varley functional responses-can provide a better description of predator feeding over a range of predator-prey abundance presence.
Particularly, Crowley-Martin functional response is much more suitable for the case where the predator feeding rate is decreased by a higher predator density even when the prey density is high [12].The Crowley-Martin functional response depends on both the predator and the prey and takes into account the interaction between predators, regardless of whether a predator is looking for prey.It is similar to the classical Beddington-DeAngelis functional response; in particular, it shows the interspecies interference with each other.Usually, it is employed to demonstrate the predator-prey stability.And, therefore, it shows ecological significance for researching the predator-prey models with Crowley-Martin functional response.A classical predatorprey system with Crowley-Martin functional response can be described by the following equations [25]: where () and V() represent the densities of the prey and the predator, respectively; parameters , , , , , and  are positive constants; the parameter  is the intrinsic growth rate of the prey and the parameter  is the mortality rate of the predator;  and / stand for the effect of capture rate and conversion factor denoting the newly born predators for each captured prey; the parameters  and  are the saturating parameters of Crowley-Martin functional response,  measures the magnitude of interference among prey, and  expresses the interference among the predators.
In recent years, the research on the predator-prey systems with Crowley-Martin functional response has attracted a lot of attention.For example, Yang studied the stability and instability of positive equilibrium in a predator-prey model with Crowley-Martin functional response and time delay [24].Asymptotic properties of a stochastic predatorprey model with Crowley-Martin functional response were studied by Liu et al., such as global existence and stochastic boundedness of the positive solution [26].A stagestructured predator-prey system with Crowley-Martin functional response was investigated, and the persistence of the system and global asymptotic stability of the positive equilibrium were assessed [27].Sivakumar et.al. analyzed a diffusive density-dependent predator-prey model with Crowley-Martin functional response and revealed Neimark-Sacker bifurcation, Turing bifurcation, and spatiotemporal patterns [28].Dong et al. investigated the positive solutions of a spatiotemporal predator-prey system with Crowley-Martin functional response and obtained a complete understanding of the uniqueness and nonuniqueness of positive solutions [29,30].Li and Wu discussed the extinction and persistence results of time-dependent positive solutions to the system [31].From previous studies, it can be seen that the predatorprey model with Crowley-Martin functional response often exhibits rich dynamics and attracts the attention of many researchers.However, the complex dynamics of the discretetime predator-prey model with Crowley-Martin functional response are still in lack of research.
The discretized form of system (1a) and (1b) was never investigated in previous studies.Since the growth, death, feeding, and migration of the predator and prey individuals always occur periodically, we can observe the dynamics of the predator-prey system by a particular time scale.This time scale can be defined by the generation span of the predator and prey populations and it measures the regeneration time of both populations.Its value is mainly determined by the population types; on the other hand, it is influenced by the change of population size or environmental conditions.In this research, we denote the time scale on which predatorprey dynamics are described and observed as parameter .Applying the forward Euler scheme with time interval , system (1a) and (1b) is now transformed to a discrete predator-prey system: where  represents the sequence number of iterations.If we give the initial time as  0 , the th iteration represents the time at  0 + .
In recent decades, the research on discrete predator-prey systems has become an important topic.Via the research of Li and Wu [31], the dynamical properties of system (1a) and (1b) have been revealed a lot.However, nonlinear characteristics of discrete system (2a) and (2b) are still scarcely known.In the literature, Ren et al. have actually investigated a Leslie-Gower type discrete predator-prey model with Crowley-Martin functional response and found the existence of Marotto chaos [32].They also controlled chaotic orbits to be a fixed point by a feedback control method.Different from the study of Ren et al. [32], this research focuses on a detailed exploration of the complex dynamics, such as periodic windows, subbifurcations, and chaos crisis, on the routes to chaos induced by the bifurcations of system (2a) and (2b).
In order to facilitate the study of system (2a) and (2b), we rewrite the equations of system (2a) and (2b) as a map; that is, In this research, the bifurcations and complex dynamics of the discrete predator-prey system with Crowley-Martin functional response are explored through map (3).The exploration is arranged as follows.Section 2 will give the stability analysis on the fixed points of map (3).Section 3 will provide bifurcation analysis and determine occurrence conditions for flip bifurcation and Neimark-Sacker bifurcation.Section 4 will present numerical simulations for complex dynamic behaviors of the discrete predator-prey system.Section 5 will describe the conclusions.

Stability Analysis
The fixed points of map (3) can be calculated through the following equations: which directly yield where  * is a positive root of the following cubic equations: According to previous results in the literature [21], the sufficient conditions for map (3) possessing a unique positive fixed point are determined as follows: To determine the stability of the fixed points of map (3), the Jacobian matrix associated with the map is applied and can be described as follows: For the Jacobian matrix of ( 8), we substitute the values of the three fixed points into it and calculate the two eigenvalues, an unstable fixed point emerges.For the fixed point ( 0 , V 0 ), the corresponding Jacobian matrix is the two eigenvalues of which are 1+ and 1−.Since 1+ is larger than 1, then ( 0 , V 0 ) is unstable.For the fixed point , the Jacobian matrix is obtained as The two eigenvalues of ( 1 , V 1 ) are 1− and 1−+/(1+ ).When the following conditions are established, the two eigenvalues are both less than one.Therefore, under conditions ( 11), ( 1 , V 1 ) is stable.The Jacobian matrix of ( 2 , V 2 ) is presented as the following matrix: Likewise, the two eigenvalues of matrix (3) are where We calculate the stability conditions of ( 2 , V 2 ), that is, | 1 | < 1 and | 2 | < 1, as follows:

Bifurcation Analysis
3.1.Flip Bifurcation.Around the stable fixed point ( 2 , V 2 ), flip bifurcation of map ( 3) is studied via the flip bifurcation theorem and center manifold theorem [27].According to the flip bifurcation theorem, the predator-prey system undergoes flip bifurcation when one of the two eigenvalues of which directly yields At such critical point  * , the two eigenvalues of ( 2 , V 2 ) are Consider the bifurcation parameter  also as a dependent variable.We then translate the fixed point ( 2 , V 2 ,  * ) to the origin via the following translation: Consequently, via Taylor expansion near the fixed point, map (3) can be transformed as where ((|| + || + |τ|) 4 ) is a polynomial with at least four orders with variables ||, ||, and |τ|.The coefficients   can be calculated by where  1 ,  2 , and  3 are the orders of the variables , , and τ, respectively, in the corresponding term of   .Here, we define  0 / 0 = 1.Applying a reversible transformation as map (20) becomes where  4 ) . ( The center manifold of map (23) at the fixed point (0, 0, 0) is then determined.According to the center manifold theorem, a central manifold at the fixed point (0, 0, 0) can be approximated by the following: where ℎ * (w, τ) is assumed to be Utilizing map (23) on both sides of z = ℎ * (w, τ) leads to simultaneously balancing the variables in (28); then, we get v  0 =  3 = 0, and the expressions of  1 and  2 are Considering the dynamics of map (23) restricted in the center manifold   (0, 0, 0), hence we obtain a onedimensional mapping as The coefficients in map (30) are described as follows: [ On the basis of map (30), the predator-prey system experiencing flip bifurcation needs the requirement of the following two conditions: which are equivalent to Hence, the occurrence of flip bifurcation at the fixed point ( 2 , V 2 ) with parameter  varying in a small neighborhood of  * needs the satisfaction of conditions ( 18) and (33).Furthermore, if  2 > 0, a period-2 orbit bifurcates from ( 2 , V 2 ) and is stable; if  2 < 0, the period-2 orbit bifurcating from ( 2 , V 2 ) is unstable.

(42b)
Calculating the two-and three-order derivatives of  1 (w, z) and  2 (w, z) at w = 0 and τ = 0, the following results are obtained: ( On the basis of map (41), the predator-prey system experiencing Neimark-Sacker bifurcation needs the requirement of the following condition: where Using the above formula, we can get in which (47) Taking above calculations together, the occurrence of Neimark-Sacker bifurcation at the fixed point ( 2 , 2 ) in the discrete predator-prey system needs the satisfaction of conditions ( 34), (39), and (46).Moreover, if  0 and 0, invariant from the fixed point for > ; > 0  > 0, a repelling closed curve bifurcates from the fixed point for  <  0 .

Numerical Results
To show the bifurcations and complex dynamics of the discrete predator-prey system, numerical simulations are carried out.We fix the parameters as  = 3,  = 2,  = 0.2,  = 1,  = 1, and  = 2 and assume that  varies.Under such conditions, the fixed point is (2.8086, 1.3436) and the critical point for flip bifurcation is  * = 0.7501.At the critical bifurcation point, the two eigenvalues are  1 = −1 and  2 = 0.8899.Considering the dynamics of map (23) restricted to the center manifold, we then obtain the values of two quantities as  1 = −2.6665and  2 = −0.0036.Therefore, the predator-prey system indeed undergoes flip bifurcation at ( 2 , V 2 ) in a small neighborhood of  =  * , and the system dynamics converge to a period-2 orbit.Taking  = 0.8 as an example, the two stable periodic points are (3.1929,1.1228) and (2.2583, 1.094).
Figure 1(a) displays the flip bifurcation diagram of the discrete predator-prey system with  varying in [0.6, 1.08].It shows a period-doubling cascade in orbits of periods 2, 4, 8, and 16 (Figures 2(a)-2(d)).As determined by the maximum Lyapunov exponent in Figure 1(b), we find that the system firstly experiences chaos at  = 0.9698.Figures 1(c)-1(e) are local amplifications of the bifurcation diagram, displaying three periodic windows occurring in the chaotic region.In these periodic windows, period-6, period-5, and period-3 orbits appear (Figures 2(g), 2(j), and 2(m)).Moreover, a subperiod-doubling cascade emerges in each periodic window, leading to interior chaotic crisis and dynamic transition from periodic behaviors to chaos (Figures 2(m)-2(p)).Figure 2 shows various attractors in the flip bifurcation diagram with the increase of  value.
The bifurcation diagram in Figure 3(a) and the maximum Lyapunov exponent diagram in Figure 3(b) suggest that the Neimark-Sacker bifurcation induces a route to chaos, leading to a dynamic transition from a fixed point, to invariant curves, with periodic windows occurring in between, to chaotic dynamics.The maximum Lyapunov exponent  diagram says that the predator-prey system firstly shows chaotic dynamics at  = 1.4072.Local amplifications in Figures 3(c)-3(h) illustrate that the invariant curves may show many isolated cycles and periodic windows can appear in the zones of both invariant curves and chaos.In the invariant curve zone, the periodic windows often emerge intermittently, whereas in the chaotic zone, the periodic windows often exhibit sub-period-doubling cascade.In the periodic windows of Figures 3(c)-3(h), we find period-9, period-13, period-20, period-22 orbits.Figure 4 explicitly displays the alternation between periodic behaviors and invariant curves or chaotic dynamics with the increase of  value.
Phase portraits for various values of  corresponding to Figure 3(a) are shown in Figure 4.It is worth noting that Figure 4(l) is a chaotic attractor in the sense of Marotto.According to the definition of snap-back repeller and Marotto chaos [32], we can prove that the fixed point ( 2 , V 2 ) is a snapback repeller by an iterative method and the emergence of Marotto chaos.In the case of  = 2,  = 2,  = 2,  = 1.75,  = 0.1,  = 0.1, and  = 1.48, we have  1,2 = 0.006 ± 1.2833 ĩ and | 1,2 | = 1.2833 > 1.Notice that the fixed point of  with the norm of all eigenvalues of (, V) exceeding 1 in magnitude is equivalent to 1++ > 0, 1-+ > 0, and  > 1.
With the variation of other parameter values (e.g., parameter ), the predator-prey system may exhibit richer dynamical behaviors in the Neimark-Sacker bifurcation diagram.When we set the parameter values as  = 2,  = 2,  = 2,  = 1.85,  = 0.1, and  = 0.1, a new Neimark-Sacker  and 4 are found in this case, such as route to chaos, invariant curves, chaotic attractors, and periodic windows, which here exhibit period-5, period-9, period-24, and period-34 orbits.However, in the periodic window of [1.3134, 1.3189], different dynamical behaviors are observed (Figure 5(c)).At about  = 1.3134, the chaotic dynamics of the predator-prey system suddenly and system with 7 (Figures 6(a) and 6(b)).On each branch of the period-7 orbit, sub-Neimark-Sacker bifurcation occurs (Figure 5(d)).In the phase portrait, this reflects as each periodic point changing to a cycle (Figure 6(c)).The Neimark-Sacker bifurcation in each branch of period-7 orbit also triggers a route to chaos like Figure 3(a), leading to 7 small chaotic subsets.Chaotic interior crisis then takes place and leads to an abrupt transition from chaotic subsets to a narrow chaotic band (Figures 5(c) and 6(d)).The chaotic band is in the range [1.317, 1.3171], and after that, the system dynamics turn to be periodic again (see period-15 orbit as shown in Figure 6(e)).Similarly, on each branch of the periodic orbit, sub-flip bifurcation occurs and induces a period-doubling cascade (Figures 5(e) and 6(e)-6(g)).Chaotic subsets and interior crisis appear sequently, taking the system entering chaotic dynamics again (Figure 6(h)).studied.Like many studies on the discrete predator-prey systems [5,6,20,21,32], we found that the discrete system can undergo flip bifurcation and Neimark-Sacker bifurcation at the stably coexistent fixed point.The two bifurcations can both trigger the route to chaos; that is, chaotic dynamics appear along with the emergence of bifurcations.On the routes to chaos, we also find a period-doubling cascade in orbits of periods 2, 4, and 8, periodic windows, invariant cycles, chaotic attractors, and Marotto chaos.Compared with the study of Ren et al. [32], who also investigated a discrete predator-prey system with Crowley-Martin functional response, many new interesting dynamic behaviors are revealed, including period-7, period-9, period-11, period-13, period-15, period-22, period-30, and period-34 orbits, sub-(inverse) Neimark-Sacker bifurcation, chaotic interior crisis, chaotic band, and intermittent periodic behaviors.Moreover, in the three-dimensional bifurcation diagram, we find the transition between flip bifurcation and Neimark-Sacker bifurcation and a critical bifurcation case between the two bifurcations.These results indicate that the discrete predatorprey system has more complex and richer dynamics than the corresponding continuous system.
The generation of complex dynamics in the discrete predator-prey system is related to two important nonlinear mechanisms, flip bifurcation and Neimark-Sacker bifurcation, which cause the system to jump from steady state to oscillatory states, such as periodic and quasiperiodic states, and trigger routes to chaos.The occurrence of these behaviors induced by the two bifurcations can ecologically explain the complex dynamic states of populations.For example, the Australian entomologist Nicholson initiated a series of longterm experiments in which he kept populations of blowflies (Lucilia cuprina) in cages under a variety of resource renewal regimes.In one experiment, he found that the numbers of adult blowflies fluctuated as perturbed periodic state for the first 400 days [33].Other experiments have investigated the long-term dynamics of freshwater organisms such as the cladoceran Daphnia in the laboratory.In one of the first such experiments, Slobodkin observed irregular fluctuations in time series that lasted between 200 and 400 days [34].Application of the approaches for analyzing time series by computing Lyapunov exponents indicates that these time series may be quasiperiodic.Moreover, chaotic dynamics were also found in the natural world or laboratory experiments.Many authors have emphasized the presence of chaotic dynamics in a wide variety of communities with interacting species [34], such as outbreaks of planktons in the phytoplanktonzooplankton system, which is one of the important predatorprey systems.With the complex dynamics found in this research, many complicated states or oscillations of predatorprey populations in the real world may be possibly explained.Also, this research could be useful for researchers who will further explore the discrete-time predator-prey systems.

Figure 2 :
Figure 2: Phase portraits for various values of  corresponding to Figure 1(a).

Figure 4 :
Figure 4: Phase portraits for various values of  corresponding to Figure 3(a).