Abstract
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, period-22, period-24, period-30, and period-34 orbits, invariant cycles, chaotic attractors, sub-flip bifurcation, sub-(inverse) Neimark-Sacker 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 and Neimark-Sacker bifurcation, showing the nonlinear characteristics of both bifurcations.
1. 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 predator-prey 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–4]. In recent decades, mathematical models have been established to analyze various complex dynamics of the predator-prey systems in various circumstances [4–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–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 discrete-time 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–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–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 predator-prey 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 prey-predator-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 predator-prey system with Crowley-Martin functional response can be described by the following equations [25]: where and 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 predator-prey 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 stage-structured 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 predator-prey model with Crowley-Martin functional response often exhibits rich dynamics and attracts the attention of many researchers. However, the complex dynamics of the discrete-time 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 predator-prey 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 , the th iteration represents the time at .
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.
2. 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, and . If and , the corresponding fixed point is considered to be stable; if or , an unstable fixed point emerges. For the fixed point , the corresponding Jacobian matrix is the two eigenvalues of which are and . Since is larger than 1, then is unstable. For the fixed point , the Jacobian matrix is obtained as The two eigenvalues of are and . When the following conditions are established,the two eigenvalues are both less than one. Therefore, under conditions (11), is stable. The Jacobian matrix of is presented as the following matrix: Likewise, the two eigenvalues of matrix (3) arewhere We calculate the stability conditions of , that is, and , as follows:
3. Bifurcation Analysis
3.1. Flip Bifurcation
Around the stable fixed point , 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 is equal to −1. Therefore, it is explicitly said that ; that is, which directly yields At such critical point , the two eigenvalues of are and . Moreover, the value of should satisfy ; that is,
Consider the bifurcation parameter τ also as a dependent variable. We then translate the fixed point to the origin via the following translation: Consequently, via Taylor expansion near the fixed point, map (3) can be transformed as where is a polynomial with at least four orders with variables , , and . The coefficients can be calculated by where , , and are the orders of the variables , , and , respectively, in the corresponding term of . Here, we define .
Applying a reversible transformation asmap (20) becomes where
The center manifold of map (23) at the fixed point is then determined. According to the center manifold theorem, a central manifold at the fixed point can be approximated by the following: where is assumed to beUtilizing map (23) on both sides of leads to simultaneously balancing the variables in (28); then, we get v , and the expressions of and are
Considering the dynamics of map (23) restricted in the center manifold , hence we obtain a one-dimensional 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 with parameter varying in a small neighborhood of needs the satisfaction of conditions (18) and (33). Furthermore, if , a period-2 orbit bifurcates from and is stable; if , the period-2 orbit bifurcating from is unstable.
3.2. Neimark-Sacker Bifurcation
When map (3) undergoes Neimark-Sacker bifurcation, the two eigenvalues of the fixed point are a pair of conjugate complex numbers whose module is one. Hence, we directly have and , which are further determined as
Under conditions (34), we translate the fixed point to the origin by the following translation: With application of Taylor expansion, map (3) is transformed to The coefficients in the map are also calculated by (21) but . The two eigenvalues at the fixed point of map (36) are also a pair of conjugate complexes with a module of one. These two eigenvalues are written as where and . The modulus of the two eigenvalues satisfies Simultaneously, the two eigenvalues should not be real numbers or imaginary numbers, which means , . Since , we have . In the meantime, . Therefore, , which is equivalent to
Through the transformation the normal form of map (36) is obtained as in which
Calculating the two- and three-order derivatives of and at and , 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
Taking the above calculations together, the occurrence of Neimark-Sacker bifurcation at the fixed point in the discrete predator-prey system needs the satisfaction of conditions (34), (39), and (46). Moreover, if and , an attracting invariant closed curve bifurcates from the fixed point for ; if and , a repelling invariant closed curve bifurcates from the fixed point for .
4. 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 , , , , , and and assume that varies. Under such conditions, the fixed point is and the critical point for flip bifurcation is . At the critical bifurcation point, the two eigenvalues are and . Considering the dynamics of map (23) restricted to the center manifold, we then obtain the values of two quantities as and . Therefore, the predator-prey system indeed undergoes flip bifurcation at in a small neighborhood of , and the system dynamics converge to a period-2 orbit. Taking as an example, the two stable periodic points are and .
Figure 1(a) displays the flip bifurcation diagram of the discrete predator-prey system with τ varying in . 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 . 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 sub-period-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.
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
(p)
Figure 3 displays Neimark-Sacker bifurcation of the discrete predator-prey system with parametric conditions given as , , , , , and . In this case, the Neimark-Sacker bifurcation emerges at the fixed point (1.3467, 0.3849) in a small neighborhood of . We can verify that, at the critical bifurcation point, and . The values of and are determined as and . Since and , we know that an attracting invariant closed curve bifurcates from the fixed point for .
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
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 . 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.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
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 is a snap-back repeller by an iterative method and the emergence of Marotto chaos. In the case of , , , , , , and , we have and . Notice that the fixed point of with the norm of all eigenvalues of exceeding 1 in magnitude is equivalent to , , and . Let , , and represent the sets of three inequalities in the front and ; then, we obtain the neighborhood . We need to find a point , satisfying and for , where . Through three iterative calculations and numerical simulations, we get the point which meets and for . Through these calculations, we proved that is a snap-back repeller, and map (3) exhibits a chaotic behavior in the sense of Marotto.
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 , , , , , and , a new Neimark-Sacker bifurcation diagram is obtained as shown in Figure 5(a). The system undergoes Neimark-Sacker bifurcation at , and it firstly enters chaotic dynamics at about (Figures 5(a) and 6(b)). Similar nonlinear characteristics to Figures 3 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 , different dynamical behaviors are observed (Figure 5(c)). At about , the chaotic dynamics of the predator-prey system suddenly disappear and the system dynamics jump to a periodic behavior with period 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 , 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)).
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Figure 7 shows the Neimark-Sacker bifurcation diagram when the parameter values are given as , , , , , and . The critical Neimark-Sacker bifurcation point is , and the first chaotic point is at around (Figures 7(a) and 7(b)). On the route to chaos, periodic windows with period-5, period-6, period-7, period-11, and period-16 orbits and narrow chaotic band are found. When we amplify the range , new dynamical behaviors are observed (Figure 7(c)). On each branch of the period-6 orbit, which occurs after the chaotic band (Figures 8(a) and 8(b)), the predator-prey system sequently undergoes sub-Neimark-Sacker bifurcation and periodic window and inverse sub-Neimark-Sacker bifurcation with the increase of τ value (Figures 7(d) and 7(e)). This process results in a dynamic transition from period-6 orbit, to invariant curves with 6 cycles, to period-24 orbit, and then the system experiences the above dynamical behaviors again but in an inverse way (Figures 8(b)–8(f)), until the emergence of chaotic behaviors (Figure 8(g)).
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
Figure 9 displays a case in which the bifurcation and route to chaos are restricted in the range of . The parametric conditions for this case are given as , , , , , and . The fixed point is a stable attractor at . With the change of value, the fixed point loses its stability through Neimark-Sacker bifurcation at . The maximum Lyapunov exponent corresponding to Figure 9(a) is shown (Figure 9(b)); we find that the predator-prey system enters chaos at . Before the chaotic region, the predator-prey system undergoes a periodic window with period-10 orbits, and as the value of τ increases, the period orbit doubles (Figures 10(a) and 10(b)), causing the system to jump to quasiperiodic behaviors (Figure 9(d)), following which the system returns to another periodic window with period halving (Figures 9(e) and 10(d)–10(g)). And then, the system experiences intermittent chaos, which suggests the abrupt appearance of chaotic dynamics. The chaotic attractor is demonstrated by the cloud of points in Figure 10(h).
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Figure 11 shows three-dimensional bifurcation diagrams in ( space and space with varying parameter values. As shown in Figure 11(a), we find that, with the increase of parameter , the Neimark-Sacker bifurcation diagram moves toward the negative τ axis. Such movement takes the decrease of Neimark-Sacker bifurcation critical point and first chaotic point and simultaneously changes the periodic windows on the route to chaos. Figure 11(b) exhibits that, with the increase of parameter , the predator-prey system presents a dynamic transition from Neimark-Sacker bifurcation to flip bifurcation. Moreover, we find a critical case between the two bifurcations, where the critical points of Neimark-Sacker bifurcation and flip bifurcation are the same and the two bifurcations combine together. With the given parametric conditions, we have in such case. The predator-prey system firstly experiences flip bifurcation and then abruptly jumps to the latter part of the Neimark-Sacker bifurcation diagram.
(a)
(b)
5. Conclusions
In this research, the complex behaviors of discrete predator-prey systems with Crowley-Martin functional response are 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 predator-prey 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 long-term 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 phytoplankton-zooplankton system, which is one of the important predator-prey systems. With the complex dynamics found in this research, many complicated states or oscillations of predator-prey 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.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors would like to thank the National Major Science and Technology Program for Water Pollution Control and Treatment (no. 2017ZX07101-002) and the Fundamental Research Funds for the Central Universities (no. JB2017069) for their support.