Research Article | Open Access
Neimark-Sacker-Turing Instability and Pattern Formation in a Spatiotemporal Discrete Predator-Prey System with Allee Effect
A spatiotemporal discrete predator-prey system with Allee effect is investigated to learn its Neimark-Sacker-Turing instability and pattern formation. Based on the occurrence of stable homogeneous stationary states, conditions for Neimark-Sacker bifurcation and Turing instability are determined. Numerical simulations reveal that Neimark-Sacker bifurcation triggers a route to chaos, with the emergence of invariant closed curves, periodic orbits, and chaotic attractors. The occurrence of Turing instability on these three typical dynamical behaviors leads to the formation of heterogeneous patterns. Under the effects of Neimark-Sacker-Turing instability, pattern evolution process is sensitive to tiny changes of initial conditions, suggesting the occurrence of spatiotemporal chaos. With application of deterministic initial conditions, transient symmetrical patterns are observed, demonstrating that ordered structures can exist in chaotic processes. Moreover, when local kinetics of the system goes further on the route to chaos, the speed of symmetry breaking becomes faster, leading to more fragmented and more disordered patterns at the same evolution time. The rich spatiotemporal complexity provides new comprehension on predator-prey coexistence in the ways of spatiotemporal chaos.
Since the pioneering work of Lotka and Volterra on population dynamics [1, 2], predator-prey relationship has become a central topic for biological studies over the last one hundred years because of its significance in many problems . To study the interactions between predator and prey in space, it is important to determine the specific forms of predator-prey functional response and spatial movement of the predator and prey populations . Theoretically, there are a fair amount of works by using nonlinear reaction-diffusion equations to study the spatiotemporal dynamics of predator-prey systems [4, 5]. Many former studies demonstrated that spatial pattern formation is an important nonlinear phenomenon in the spatially extended predator-prey systems [6, 7]. Spatial pattern is a ubiquitous scenario in nature and often takes the role of modifying temporal dynamics and regulating population stability . After the investigations in about four decades, pattern formation of predator-prey systems has been widely concerned and become one of timely central topics in biological science [9, 10]. It has been a great interest to understand spatial patterns and related mechanisms of interacting species in predator-prey systems.
The researchers have developed various reaction-diffusion predator-prey models in order to study the pattern formation of predator-prey systems in varying circumstances, such as different functional responses, growth forms, and population motions. In literature, it is found that many types of functional responses can contribute greatly to the spatiotemporal complexity of predator-prey systems, such as Beddington-DeAngelis type, Leslie-Gower type, Holling type IV, and ratio-dependent type functional response. Upadhyay et al. discovered that zooplankton-phytoplankton systems with Holling type IV functional response have a common mechanism as wave of chaos, which can result in the self-organization of complex patterns . Banerjee and Petrovskii demonstrated that a ratio-dependent predator-prey system can develop patterns both inside and outside of the Turing parameter domain and Hopf bifurcation is essential for the onset of spatiotemporal chaos .
Different growth forms of prey can also exert influences on pattern formation. Generally, the most common growth forms include logistic form, exponential form, and the form with Allee effect. Many works have considered continuous-time models for the predator-prey systems with Allee effect. Numerical results in Rao and Kang  indicated that the spatial Michaelis-Menten type predation model with Allee effects exhibits rich dynamical transitions from wave growth not only to striped-like/spots patterns, but also to regular patterns, which implied that synergistic effects of diffusion and Allee effects can contribute to the patterns of the solutions in different ways. The dynamics of a reaction-diffusion predator-prey system with strong Allee effect in the prey population has been considered in Wang et al.’s research , which revealed that the impact of Allee effect essentially increases the system spatiotemporal complexity. More abundant investigations and results relating to the Allee effect influencing predator-prey pattern formation can be found in Petrovskii et al. , Sun et al. , and Rodrigues et al. .
Moreover, many researchers further found that the population motions play a key role in spatiotemporal pattern formation of predator-prey systems. As important natural phenomena, cross-diffusion and convection of the predator and prey are often considered and studied in spatially extended predator-prey systems. Ghorai and Poria found that the cross-diffusion supports the formation of a wide variety of spatial and spatiotemporal patterns . Ling et al. explored the impacts of cross-diffusion on the formation of spatial patterns in a ratio-dependent predator-prey system with zero-flux boundary conditions and determined that the cross-diffusion can trigger the emergence of spatial patterns which is however impossible under the same conditions when cross-diffusion is absent . In the predator-prey systems with predator and/or prey convection, the spatial patterns are often characterized by stripes .
The nonlinear mechanism of pattern formation, Turing instability, can date back to the pioneering work of Turing who studied a reaction-diffusion model which included two chemical species, an activator and an inhibitor, in 1952 . Levin and Segel originally extended the mechanism of Turing instability to ecology and theoretically investigate how diffusion leads to emergence of spatial patterns . Their method becomes a routine framework in the studying of the spatial pattern formation of predator-prey systems. Much research effort has been devoted to the mechanisms behind the emergence of spatial patterns and a variety of pattern formation mechanisms have been found. These mechanisms can be basically classified into two categories. One category is that spatial patterns emerge from locally stable kinetic dynamics due to diffusion and/or convection, often called pure Turing instability, whereas the other category is that spatial patterns emerge from locally unstable kinetic dynamics (such as limit cycle and chaos) through spatially heterogeneous perturbations, and one of the mechanisms is known as Hopf-Turing instability .
Most of the reaction-diffusion predator-prey models developed in literature are time-and space-continuous. Recently, a new type of mathematical model was introduced into biology to study the pattern formation, i.e., coupled map lattices (CMLs). In the 1980s, the CML was used by Kaneko to explore the spatial-temporal structure of coupled logistic map, he found very complex spatiotemporal dynamics, including frozen chaos, defect turbulence, spatiotemporal intermittency, and fully developed spatiotemporal chaos [21–23]. Recently in 2012, Mistro et al. , Rodrigues et al. , and Punithan et al.  demonstrated that the CML model shows importance in understanding the predator-prey dynamics on a fragmented habitat, which can exhibit spatiotemporal chaos and spatiotemporal multistability. Huang and Zhang further revealed that considering the combining effect of pattern formation and bifurcations in the CML can result in rich and complicated dynamical behaviors and nonlinear characteristics . Diverse predator-prey pattern self-organizations and complex pattern transitions can be induced by the coupled effects of Turing instability and Neimark-Sacker instability [10, 27, 28].
On the basis of former research works, we aim at studying the complex dynamics of a space- and time-discrete predator-prey system under the influence of Allee effect with the application of CML. For the discrete predator-prey system considered, four properties should be specified. First, the habitat where predator and prey populations dwell is distinctly patchy or even fragmented , and the individuals migrate among different patches. Second, the births and deaths of populations all occur at discrete times, or within certain intervals of time . Third, the population numbers are small and in this case the application of CML model may be the most appropriate. Fourth, the Allee effect occurring in the predator-prey system represents a positive correlation between the population density and the per capita growth rate in population dynamics. It deserves to explore the spatiotemporal structures of such predator-prey system under the combining effect of Neimark-Sacker instability and Turing instability, which can help to further understand the spatiotemporal complexity of the discrete predator-prey system and are still few investigated in literature. The structure of the present research work is arranged as follows. Section 2 describes the CML model of the discrete predator-prey system with Allee effect and presents the stability analysis. In Section 3, we perform Neimark-Sacker bifurcation analysis and Turing instability analysis. Section 4 provides the numerical results of the complex dynamical behaviors, Neimark-Sacker bifurcation and pattern formation. In Section 5, we provide the discussion and conclusions of the research.
2. CML Model and Stability Analysis
2.1. CML Model of the Spatiotemporal Discrete Predator-Prey System
In this research, we investigate a predator-prey system with Allee effect and Michaelis-Menten type functional response. Allee effect, which was introduced by as well as named after Allee, refers to a positive correlation between population density and the per capita growth rate [3, 30]. Recently, Allee effect has been extensively explored to investigate its impact on dynamical features of predator-prey systems [12, 13, 15]. According to the research of Rao and Kang , the Michaelis-Menten type predator-prey system can exhibit great spatiotemporal complexity under the influence of Allee effect. Therefore, it holds significance and importance to further investigate the spatiotemporal dynamics of the discrete form of such predator-prey system, which may show more complicated nonlinear characteristics.
The Michaelis-Menten type reaction-diffusion predator-prey system with Allee effect in prey can be described as follows :where and are population densities of the prey and predator at time , respectively, is the maximum per capita growth rate of the prey, is the prey carrying capacity, and () is the “Allee threshold”. The value of is considered as a measure of the intensity of the Allee effect: the less the value of is, the less prominent is the Allee effect  (Petrovskii et al., 2002). More specifically, the Allee effect is called “strong” if (when the growth rate becomes negative for ), while it is called “weak” if , and no Allee effect exists when  (Morozov et al., 2004). is the predation rate, is the food utilization coefficient, and is the mortality rate of the predator.
Considering dimensionless variables with the following scaling:and simultaneously employing the population diffusive terms on (1), then we obtain the following reaction-diffusion equations:in whichIn the above equations, is the relative intrinsic growth rate of the prey; describes the Allee threshold; represents the energy conversion rate from prey to predator; denotes the relative death rate of predator; denotes the Laplacian operator and describes the population diffusion in space; nonnegative constants and are the diffusion coefficients of and , respectively, measuring the speed of individual movements.
The CML model of spatiotemporal discrete predator-prey system is developed based on discretizing (3). For making the discretization, we give a time interval , a space interval , and a two-dimensional rectangular domain which includes lattices (the length of each lattice is ). Each lattice represents one site and is ascribed to two numbers, i.e., the prey density and the predator density. Correspondingly, two state variables, and ( and ), are newly defined in the discrete space and time, representing the prey density and the predator density in the site and at th iteration. If we give an initial time , then the time at th iteration is . According to the framework of coupled map lattice described in previous approach [10, 15, 24–28], the dynamics of the spatiotemporal discrete predator-prey at each discrete step from to iteration consists of two distinctly different stages, (a) the dispersal stage and (b) the “reaction” stage. The dispersal stage can be obtained by discretizing the spatial terms of (3), described by the following equations:where and are the prey and predator densities after dispersal; denotes the discrete form of the Laplacian operator, i.e.,
Via discretizing the nonspatial part of (3), the predator-prey reaction stage described in the CML model is expressed by the following equations:where and are the reaction functions determined by local inter- and intraspecific interactions, given by
Equations (5)–(8) describe the governing equations of the CML model of the spatiotemporal discrete predator-prey system. For applying the CML model, periodic boundary conditions are utilized in this research .
2.2. Stability Analysis on Homogeneous Stationary States
A homogeneous stationary state of the spatiotemporal discrete predator-prey system is a state where the prey and predator densities are homogeneous in space and simultaneously keep steady in time. Its stability determines whether the predator-prey system can stably stay at the corresponding homogeneous stationary state under some external disturbances. Generally, the stationary states of the discrete predator-prey system satisfy the following conditions: i.e.,for all of , , and . Under such conditions, the predator-prey equations can be described by the following dynamic equations:
For convenience of the following analysis, (10) are rewritten in the form of a map,
Mathematically, the homogeneous stationary states of the discrete predator-prey system are exactly equivalent to the fixed points of the map (11). Therefore, the results for solving the fixed points of (11) and the corresponding stability analysis are consistent with that on the homogeneous stationary states. According to the definition of the fixed point of a map, direct calculation yields five nonnegative fixed points for map (11),in whichWhen and are established, the two fixed points and are positive. Since we intend to study the predator-prey coexistence in this research, hereafter we always presume these two conditions of positive fixed points are satisfied.
In order to determine the stability of the fixed points, the method of Jacobian matrix is often applied. The Jacobian matrix associated with map (11) at any point is described as follows:
Substitute the value of each fixed point into (14) and then calculate the two eigenvalues of corresponding Jacobian matrix. We denote these two eigenvalues as and . If and , then the corresponding fixed point is either a stable node or a stable focus; if and , then the corresponding fixed point is either an unstable node or an unstable focus; if and or and , then the corresponding fixed point is a saddle point, which is also unstable.
Substituting into (14), then we obtain the Jacobian matrix asObviously the two eigenvalues of are , . If we want the fixed point to be stable, the absolute values of and must be less than one. According to this, the conditions for that is stable are determined as .
Similarly for the fixed points and , we haveIt can be found that one of the eigenvalues of (16) is the same, described as . Because we presume that , we directly have , suggesting and are unstable.
According to and , the conditions for that or is stable are determined as and . With the given parametric conditions, the stability of and will be further determined in the numerical simulations section. In the following theoretical analysis, we denote the stable coexistent fixed point as .
3. Analysis on Neimark-Sacker Bifurcation and Turing Instability
3.1. Neimark-Sacker Bifurcation Analysis
Besides the homogeneous stationary states, the discrete predator-prey can also generate many other complex spatially homogeneous states, which always keep oscillating with time. To determine the transition from homogeneous stationary states to homogeneous oscillating states, bifurcation analysis is one of the most reliable and efficient ways. It should be noticed that such transition in the spatiotemporal discrete predator-prey system is mathematically equivalent to the transition from the fixed point to other attractors in the map (11). Generally, discrete system may exhibit various types of bifurcations, including fold bifurcation, flip bifurcation, and Neimark-Sacker bifurcation . The Neimark-Sacker bifurcation is the discrete analogue of the Hopf bifurcation that occurs in continuous systems and is highly relevant to the discrete dual congestion control algorithm . Often in literature, the explicit Neimark-Sacker bifurcation is also called Hopf bifurcation for map, the calculations of which can be seen in Cheng and Cao .
When map (11) undergoes Neimark-Sacker bifurcation, the stable fixed point will change to an invariant closed curve. In the following bifurcation analysis, parameter is chosen as the bifurcation parameter. According to the Neimark-Sacker bifurcation theorem, the occurrence of such bifurcation needs the satisfaction of many conditions . Firstly, Neimark-Sacker bifurcation takes place if the eigenvalues in (18) are conjugate complex numbers and their moduli are one. That results inCalculating on (19) yields the following conditions:
Under the satisfaction of conditions (21), the fixed point of map (11) is then translated to the origin by the following translation:Thus, the map (11) is transformed into a new map aswhere and describes a function with order at least four in the variables , and the coefficients in (23) are described as follows:
When Neimark-Sacker bifurcation occurs, the modulus of the eigenvalues (25) should be varying with the bifurcation parameter. This requiresA direct calculation on leads toSince , then we find that the value of is indeed not equal to zero. Another condition for the occurrence of Neimark-Sacker bifurcation is that the eigenvalues (25) are exactly neither real numbers nor imaginary numbers, implyingwhich is equivalent to . From (21), we find that . It means , i.e.,
On the basis of map (23), the canonical form is studied to obtain the last condition for the occurrence of Neimark-Sacker bifurcation. Applying the following transformation on map (23), then we obtainwhere
The second-order and third-order partial derivatives of and at (0, 0) are calculated as
According to the bifurcation theorem , if the map (32) undergoes Neimark-Sacker bifurcation at around , the following discriminatory quantity must not be zero, i.e.,in whichA calculation on the quantity giveswhere
Based on the above calculations, the occurrence of Neimark-Sacker bifurcation in the discrete predator-prey system needs the satisfaction of conditions (21), (30), and (38). Moreover, if and , then an attracting invariant closed curve bifurcates from the fixed point for ; otherwise, if and , a repelling invariant closed curve bifurcates from the fixed point for .
3.2. Turing Instability Analysis
When spatially heterogeneous perturbations occur at the homogeneous states, the spatiotemporal discrete predator-prey system may converge to heterogeneous states . To determine the occurrence conditions for the stable heterogeneous states, Turing instability analysis is performed. Firstly, spatially heterogeneous perturbations are introduced to perturb the stable homogeneous stationary state , given byIt should be noticed that
Based on the calculations in former approach , the eigenvalues of the operator are determined asin which
Substituting (40) and (41) into the spatiotemporal discrete predator-prey system leads towhere . The two-order terms in (44) can be ignored when the perturbations are small. Using the corresponding eigenfunction of the eigenvalues to multiply (44) getsSumming (45) for all of and obtains
According to previous research works [10, 26–28], the occurrence of Turing instability demands that existing one group of to make the maximum value of larger than one. DefineThus, when , Turing instability occurs; when , the spatiotemporal discrete predator-prey system returns to the homogeneous states under spatially heterogeneous perturbations. Moreover, the threshold condition for the occurrence of Turing bifurcation requires . Consequently, there exists a critical value satisfying .
4. Numerical Results
The purpose of numerical simulations is to demonstrate the Neimark-Sacker bifurcation as well as pattern formation under the influence of Neimark-Sacker-Turing instability. For such demonstration, the parametric conditions should be provided firstly. Based on former research work of Rao and Kang , a group of ecologically feasible parameter values can be given as = 1.63, = 0.1, = 0.32, = 0.22, = 0.6, = 1.2, and = 200. The values of parameters and should be set to ensure nonnegativity of the state variables and . In the following simulations, we let = 10 and as a varying parameter. Then according to previous research works [10, 27], the nonnegativity of and needs that and are both less than 0.5.
Under the above provided parametric conditions, the fixed point of the discrete predator-prey system is determined as (0.6538, 02972), and the critical point for Neimark-Sacker bifurcation is calculated as = 4.9444. When , the fixed point of the system is stable. Taking = 4.5 as an example, the two eigenvalues of the fixed point are = 0.8307 + 0.5287i, = 0.8307 − 0.5287i, the modulus of which is less than one, suggesting is a stable focus. When , we get the two eigenvalues of as = 0.8140 + 0.5809i, = 0.8140 − 0.5809i, and = 1. Then we check conditions (28), (30), and (38) for Neimark-Sacker bifurcation; we have = 0.0152 ≠ 0, = −1.6280 ≠ 0,1 and = −179.8092 ≠ 0. Hence, according to the Neimark-Sacker bifurcation theorem, the discrete predator-prey system indeed undergoes Neimark-Sacker bifurcation at . Furthermore, since < 0 and > 0, we know that an attracting invariant closed curve bifurcates from the fixed point for , and the stable point loses its stability.
Figure 1 shows the Neimark-Sacker bifurcation diagram of the discrete predator-prey system. From the bifurcation diagram and corresponding maximum Lyapunov exponent, it is known that the Neimark-Sacker bifurcation induces a route to chaos, leading to a dynamic transition from fixed point, to invariant curves, with periodic windows occurring in-between, to chaotic dynamics. Invariant curves suggest that the discrete predator-prey system follows dynamical behaviors which are homogeneous in space and quasiperiodically oscillating in time. In the periodic windows, we find period 11, 20, 31, 41, 51, 62, 72, and 82 orbits, revealing spatially homogeneous and periodically oscillating behaviors of the system. The maximum Lyapunov exponent diagram in Figure 1(b) says that the predator-prey system begins entering spatially homogeneous chaotic dynamics with the increase of the value of and finally enters chaotic dynamics at around = 5.778. Figure 2 explicitly displays the three typical dynamical behaviors in the phase portraits, i.e., invariant closed curve, periodic orbit (period 31), and chaotic attractor.
(a) τ = 5.1
(b) τ = 5.525
(c) τ = 5.83
When heterogeneous perturbations occur on the spatially homogeneous states, the discrete predator-prey system may experience Turing instability and the system dynamics converges to stable spatially heterogeneous states. Under the previous given parametric conditions, the critical point for Turing bifurcation is overlapping with the Neimark-Sacker bifurcation critical point. Via calculation, we find the values corresponding to the three cases in Figure 2 are (5.1) = 1.0060, (5.525) = 1.0241, and (5.83) = 1.0385, respectively. Therefore, it is verified that Turing instability can indeed take place on the route to chaos induced by Neimark-Sacker bifurcation. For the three cases in Figure 2, Turing instability is spatial symmetry breaking on spatially homogeneous periodically, quasiperiodically, and chaotically oscillating states, resulting in different pattern formation. For simulating the pattern formation, two types of initial conditions are applied, i.e., stochastic initial condition and deterministic initial condition. Stochastic initial condition is given by perturbing the homogeneous stationary state (, ) with small random perturbations. And deterministic initial condition is given by the following:
Figures 3–5 show the evolution of spatiotemporal patterns of prey population (the predator patterns exhibit similar configuration) corresponding to spatial symmetry breaking occurring on the three dynamical behaviors in Figure 2, starting from the stochastic initial condition. At the beginning, the prey patterns evolve to a spatial configuration which is composed of large patches of higher and lower prey densities (such as Figures 3(b), 4(c), and 5(c)). As the time progresses, the large prey patches gradually fragment and change to smaller prey patches, the shape of which may present as stripes, arcs, spots, rings, etc. These smaller patches display disordered distribution in space and may combine together to form complex structure. After long-term evolution of patterns, we find that the prey patterns are always oscillating similarly in time and never reach equilibrium configuration.
(a) = 1
(b) = 500
(c) = 1000
(d) = 1500
(e) = 2000
(f) = 5000
(a) = 1
(b) = 10
(c) = 100
(d) = 500
(e) = 1000
(f) = 5000
(a) = 1
(b) = 10
(c) = 100
(d) = 500
(e) = 2000
(f) = 5000
As exhibited in Figures 3–5, pattern formation in the three different cases presents similar evolution trend and process. However, two main differences are found in these pattern evolutions. When the local kinetics of the discrete predator-prey system goes further on the route to chaos, first, the pattern evolution speed becomes faster, and second, the prey patterns evolving to long time shows more fragmented, more disordered characteristics with smaller patches.
In order to exhibit the dynamical trajectory of pattern evolution, Figure 6 is plotted using the change of state variables at one lattice in Figures 3–5, which here is chosen as lattice (100, 100), with the increase of time. As shown in Figures 6(a)–6(c), three new attractors emerge in the phase portraits. Comparing with Figure 2, it can be found that, under the influence of Neimark-Sacker-Turing instability, the invariant cycle and periodic orbit break and transform to more complex attractors. The clouds of random points in Figures 6(b)-6(c) suggest that and may experience almost all states in the black area. Moreover, the evolution process of prey patterns in the three cases as time progresses presents disordered waves, as displayed by Figures 6(d)–6(f).
(d) = 5.1
(e) = 5.525
(f) = 5.83
To explore the nonlinear characteristics of pattern evolution process in Figures 3–5, sensitivity analysis is carried out with the variations of initial conditions and parameter and the results are shown in Figure 7. In Figures 7(a)–7(c), we apply two close initial conditions, and , the difference of which is given as in merely one lattice, and then we observe the change of in another lattice. At the beginning, stays at zero; after a certain time of pattern evolution, the amplitude of suddenly bursts and follows chaotic oscillations. For the pattern formation in the three cases of Figures 3–5, this phenomenon suggests the occurrence of spatiotemporal chaos, which can result in the self-organization of numerous different spatiotemporal patterns with the variation of initial conditions. Simultaneously, it also reveals that the discrete predator-prey system holds the property of multistability when Neimark-Sacker instability and Turing instability both occur. Figures 7(d)–7(f) show the sensitivity of pattern evolution to the variation of parameter , which is controlled as 0.00001. Even with such tiny variation of parameter , the pattern evolution process in Figures 3–5 will follow totally different trajectories.
(d) = 5.1
(e) = 5.525
(f) = 5.83
The above sensitivity analysis suggests that change of initial conditions may lead to the self-organization of different patterns. Figures 8–10 demonstrate the evolution of prey patterns with the application of deterministic initial condition as described in (52). In these three cases, the pattern evolution starts from a square in the background. Then, a target wave pattern emerges replacing the square. As the rings of target wave pattern gradually expand, a new different structure may emerge at the core of the target and the target wave breaks from the corners, leading to the symmetry breaking increasing with time. In the pattern evolution process, symmetrical irregular transient patterns dominate the pattern evolution process all along, such as in Figure 8(g), from which one can see that the degree of symmetry is still strong as vertical, horizontal, and central symmetry. However, at large iteration time (e.g., = 200000), the patterns for the three cases all become disordered and irregular and the spatial symmetry of the patterns reduces to be weak, like the horizontal symmetry observed in Figures 8(i), 9(f), and 10(f).
(a) = 1
(b) = 500
(c) = 2000
(d) = 3000
(e) = 5000
(f) = 7000
(g) = 10000
(h) = 100000
(i) = 200000
(a) = 1
(b) = 500
(c) = 1000
(d) = 1500
(e) = 2000
(f) = 200000