# Neimark-Sacker-Turing Instability and Pattern Formation in a Spatiotemporal Discrete Predator-Prey System with Allee Effect

Authors: Huayong Zhang; Xuebing Cong; Tousheng Huang; Shengnan Ma; Ge Pan
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## Abstract

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.


---

## Body

## 1. Introduction

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 [3]. 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 [3]. 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 [8]. 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 [11]. 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 [6].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 [3] 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 [12], 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. [13], Sun et al. [14], and Rodrigues et al. [15].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 [16]. 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 [17]. In the predator-prey systems with predator and/or prey convection, the spatial patterns are often characterized by stripes [18].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 [19]. Levin and Segel originally extended the mechanism of Turing instability to ecology and theoretically investigate how diffusion leads to emergence of spatial patterns [9]. 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 [20].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. [24], Rodrigues et al. [15], and Punithan et al. [25] 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 [26]. 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 [24], 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 [29]. 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 [3], 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 [3]:(1)dNdT=4ωK-N02NN-N0K-N-ANPN+P,dPdT=μANPN+P-MP,where N(T) and P(T) are population densities of the prey and predator at time T≥0, respectively, ω is the maximum per capita growth rate of the prey, K is the prey carrying capacity, and N0 (N0<K) is the “Allee threshold”. The value of N0 is considered as a measure of the intensity of the Allee effect: the less the value of N0 is, the less prominent is the Allee effect [13] (Petrovskii et al., 2002). More specifically, the Allee effect is called “strong” if 0<N0<K (when the growth rate becomes negative for N<N0), while it is called “weak” if -K<N0≤0, and no Allee effect exists when N0≤K [31] (Morozov et al., 2004). A is the predation rate, μ is the food utilization coefficient, and M is the mortality rate of the predator.Considering dimensionless variables with the following scaling:(2)u=NK,v=PK,t=AT,and simultaneously employing the population diffusive terms on (1), then we obtain the following reaction-diffusion equations:(3)∂u∂t=auu-b1-u-uvu+v+D1Δu,∂v∂t=cuvu+v-dv+D2Δv,in which(4)a=4ωK2AK-N02,b=N0K,c=μ,d=MA.In the above equations, a is the relative intrinsic growth rate of the prey; b∈(-1,1) describes the Allee threshold; c∈(0,1] represents the energy conversion rate from prey to predator; d denotes the relative death rate of predator; Δ denotes the Laplacian operator and describes the population diffusion in space; nonnegative constants D1 and D2 are the diffusion coefficients of u and v, 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 n×n 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, u(i,j,t) and v(i,j,t) (i,j∈1,2,3,…,n and t∈N), are newly defined in the discrete space and time, representing the prey density and the predator density in the (i,j) site and at tth iteration. If we give an initial time t0, then the time at tth iteration is t0+tτ. 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 t to t+1 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:(5)ui,j,t′=ui,j,t+τδ2D1Δdui,j,t,vi,j,t′=vi,j,t+τδ2D2Δdvi,j,t,where ui,j,t′ and vi,j,t′ are the prey and predator densities after dispersal; Δd denotes the discrete form of the Laplacian operator, i.e.,(6)Δdui,j,t=ui+1,j,t+ui-1,j,t+ui,j+1,t+ui,j-1,t-4ui,j,t,Δdvi,j,t=vi+1,j,t+vi-1,j,t+vi,j+1,t+vi,j-1,t-4vi,j,t.Via discretizing the nonspatial part of (3), the predator-prey reaction stage described in the CML model is expressed by the following equations:(7)ui,j,t+1=f1ui,j,t′,vi,j,t′,vi,j,t+1=g1ui,j,t′,vi,j,t′,where f1 and g1 are the reaction functions determined by local inter- and intraspecific interactions, given by(8)f1u,v=u+τauu-b1-u-uvu+v,g1u,v=v+τcuvu+v-dv.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 [10].


### 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.,(9)Δdui,j,t=0,Δdvi,j,t=0,for all of i, j, and t. Under such conditions, the predator-prey equations can be described by the following dynamic equations:(10)ut+1=ut+τautut-b1-ut-utvtut+vt,vt+1=vt+τcutvtut+vt-dvt.For convenience of the following analysis, (10) are rewritten in the form of a map,(11)uv⟶u+τauu-b1-u-uvu+vv+τcuvu+v-dv.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),(12)u0,v0:0,0,u1,v1:1,0,u2,v2:b,0,u3,v3:u+,c-ddu+,u4,v4:u-,c-ddu-,in which(13)u±=12b+1±b+12-4b+c-dac.When c>d and b-1>2c-d/ac are established, the two fixed points (u3,v3) and (u4,v4) 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:(14)Ju,v=1+τ-3au2+2ab+1u-ab-v2u+v2-τu2u+v2τcv2u+v21+τcu2u+v2-d.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 λ1 and λ2. If |λ1|<1 and |λ2|<1, then the corresponding fixed point is either a stable node or a stable focus; if |λ1|>1 and |λ2|>1, then the corresponding fixed point is either an unstable node or an unstable focus; if |λ1|>1 and |λ2|<1 or |λ1|<1 and |λ2|>1, then the corresponding fixed point is a saddle point, which is also unstable.Substituting(u0,v0) into (14), then we obtain the Jacobian matrix as(15)Ju0,v0=1-τab001-τd.Obviously the two eigenvalues of J(u0,v0) are λ1=1-τab, λ2=1-τd. If we want the fixed point (u0,v0) to be stable, the absolute values of λ1 and λ2 must be less than one. According to this, the conditions for that (u0,v0) is stable are determined as 0<τ<min⁡{2/ab,2/d}.Similarly for the fixed points(u1,v1) and (u2,v2), we have(16)Ju1,v1=1+τab-1-τ01+τc-d,Ju2,v2=1+τab1-b-τ01+τc-d.It can be found that one of the eigenvalues of (16) is the same, described as λ2=1+τ(c-d). Because we presume that c>d, we directly have λ2>1, suggesting (u1,v1) and (u2,v2) are unstable.As for the fixed points(u3,v3) and (u4,v4), we can write their Jacobian matrix as(17)J=1+τ-ab+1u±+2ab+3c-dc-c-d2c2-τd2c2τc-d2c1+τdd-cc.For convenience, we demote the four elements of matrix (17) as a11, a12, a21, and a22. The two eigenvalues of (17) are obtained as(18)λ1,2=12-p±p2-4q,in which(19)p=-a11-a22=-2+τab+1u±-2ab-c-d2c+d-cdc2,q=a11a22-a12a21=1-τab+1u±-2ab-c-d2c+d-cdc2+τ2dc-dcab+1u±-2ab-2c-dc.According to|λ1|<1 and |λ2|<1, the conditions for that (u3,v3) or (u4,v4) is stable are determined as q<1 and |p|<1+q. With the given parametric conditions, the stability of (u3,v3) and (u4,v4) will be further determined in the numerical simulations section. In the following theoretical analysis, we denote the stable coexistent fixed point as (u∗,v∗).


## 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 [3], 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 [3]:(1)dNdT=4ωK-N02NN-N0K-N-ANPN+P,dPdT=μANPN+P-MP,where N(T) and P(T) are population densities of the prey and predator at time T≥0, respectively, ω is the maximum per capita growth rate of the prey, K is the prey carrying capacity, and N0 (N0<K) is the “Allee threshold”. The value of N0 is considered as a measure of the intensity of the Allee effect: the less the value of N0 is, the less prominent is the Allee effect [13] (Petrovskii et al., 2002). More specifically, the Allee effect is called “strong” if 0<N0<K (when the growth rate becomes negative for N<N0), while it is called “weak” if -K<N0≤0, and no Allee effect exists when N0≤K [31] (Morozov et al., 2004). A is the predation rate, μ is the food utilization coefficient, and M is the mortality rate of the predator.Considering dimensionless variables with the following scaling:(2)u=NK,v=PK,t=AT,and simultaneously employing the population diffusive terms on (1), then we obtain the following reaction-diffusion equations:(3)∂u∂t=auu-b1-u-uvu+v+D1Δu,∂v∂t=cuvu+v-dv+D2Δv,in which(4)a=4ωK2AK-N02,b=N0K,c=μ,d=MA.In the above equations, a is the relative intrinsic growth rate of the prey; b∈(-1,1) describes the Allee threshold; c∈(0,1] represents the energy conversion rate from prey to predator; d denotes the relative death rate of predator; Δ denotes the Laplacian operator and describes the population diffusion in space; nonnegative constants D1 and D2 are the diffusion coefficients of u and v, 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 n×n 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, u(i,j,t) and v(i,j,t) (i,j∈1,2,3,…,n and t∈N), are newly defined in the discrete space and time, representing the prey density and the predator density in the (i,j) site and at tth iteration. If we give an initial time t0, then the time at tth iteration is t0+tτ. 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 t to t+1 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:(5)ui,j,t′=ui,j,t+τδ2D1Δdui,j,t,vi,j,t′=vi,j,t+τδ2D2Δdvi,j,t,where ui,j,t′ and vi,j,t′ are the prey and predator densities after dispersal; Δd denotes the discrete form of the Laplacian operator, i.e.,(6)Δdui,j,t=ui+1,j,t+ui-1,j,t+ui,j+1,t+ui,j-1,t-4ui,j,t,Δdvi,j,t=vi+1,j,t+vi-1,j,t+vi,j+1,t+vi,j-1,t-4vi,j,t.Via discretizing the nonspatial part of (3), the predator-prey reaction stage described in the CML model is expressed by the following equations:(7)ui,j,t+1=f1ui,j,t′,vi,j,t′,vi,j,t+1=g1ui,j,t′,vi,j,t′,where f1 and g1 are the reaction functions determined by local inter- and intraspecific interactions, given by(8)f1u,v=u+τauu-b1-u-uvu+v,g1u,v=v+τcuvu+v-dv.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 [10].


## 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.,(9)Δdui,j,t=0,Δdvi,j,t=0,for all of i, j, and t. Under such conditions, the predator-prey equations can be described by the following dynamic equations:(10)ut+1=ut+τautut-b1-ut-utvtut+vt,vt+1=vt+τcutvtut+vt-dvt.For convenience of the following analysis, (10) are rewritten in the form of a map,(11)uv⟶u+τauu-b1-u-uvu+vv+τcuvu+v-dv.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),(12)u0,v0:0,0,u1,v1:1,0,u2,v2:b,0,u3,v3:u+,c-ddu+,u4,v4:u-,c-ddu-,in which(13)u±=12b+1±b+12-4b+c-dac.When c>d and b-1>2c-d/ac are established, the two fixed points (u3,v3) and (u4,v4) 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:(14)Ju,v=1+τ-3au2+2ab+1u-ab-v2u+v2-τu2u+v2τcv2u+v21+τcu2u+v2-d.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 λ1 and λ2. If |λ1|<1 and |λ2|<1, then the corresponding fixed point is either a stable node or a stable focus; if |λ1|>1 and |λ2|>1, then the corresponding fixed point is either an unstable node or an unstable focus; if |λ1|>1 and |λ2|<1 or |λ1|<1 and |λ2|>1, then the corresponding fixed point is a saddle point, which is also unstable.Substituting(u0,v0) into (14), then we obtain the Jacobian matrix as(15)Ju0,v0=1-τab001-τd.Obviously the two eigenvalues of J(u0,v0) are λ1=1-τab, λ2=1-τd. If we want the fixed point (u0,v0) to be stable, the absolute values of λ1 and λ2 must be less than one. According to this, the conditions for that (u0,v0) is stable are determined as 0<τ<min⁡{2/ab,2/d}.Similarly for the fixed points(u1,v1) and (u2,v2), we have(16)Ju1,v1=1+τab-1-τ01+τc-d,Ju2,v2=1+τab1-b-τ01+τc-d.It can be found that one of the eigenvalues of (16) is the same, described as λ2=1+τ(c-d). Because we presume that c>d, we directly have λ2>1, suggesting (u1,v1) and (u2,v2) are unstable.As for the fixed points(u3,v3) and (u4,v4), we can write their Jacobian matrix as(17)J=1+τ-ab+1u±+2ab+3c-dc-c-d2c2-τd2c2τc-d2c1+τdd-cc.For convenience, we demote the four elements of matrix (17) as a11, a12, a21, and a22. The two eigenvalues of (17) are obtained as(18)λ1,2=12-p±p2-4q,in which(19)p=-a11-a22=-2+τab+1u±-2ab-c-d2c+d-cdc2,q=a11a22-a12a21=1-τab+1u±-2ab-c-d2c+d-cdc2+τ2dc-dcab+1u±-2ab-2c-dc.According to|λ1|<1 and |λ2|<1, the conditions for that (u3,v3) or (u4,v4) is stable are determined as q<1 and |p|<1+q. With the given parametric conditions, the stability of (u3,v3) and (u4,v4) will be further determined in the numerical simulations section. In the following theoretical analysis, we denote the stable coexistent fixed point as (u∗,v∗).


## 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 [29]. 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 [32]. 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 [29].When map (11) undergoes Neimark-Sacker bifurcation, the stable fixed point (u∗,v∗) 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 [33]. Firstly, Neimark-Sacker bifurcation takes place if the eigenvalues in (18) are conjugate complex numbers and their moduli are one. That results in(20)p2-4q<0,q=1.Calculating on (19) yields the following conditions:(21)0<ab+1u±-2ab-c-d2c+d-cdc2<4τ,τ=τ0=cdc-d-1-cacb+1u±-2abc-2c-d.Under the satisfaction of conditions (21), the fixed point (u∗,v∗) of map (11) is then translated to the origin by the following translation:(22)w=u-u∗,z=v-v∗.Thus, the map (11) is transformed into a new map as(23)wz⟶a11w+a12z+a13w2+a14z2+a15wz+a16w3+a17w3+a18w2z+a19wz2+Oζ4a21w+a22z+a23w2+a24z2+a25wz+a26w3+a27w3+a28w2z+a29wz2+Oζ4,where ζ=w+z and Oζ4 describes a function with order at least four in the variables w,z, and the coefficients in (23) are described as follows:(24)a13=τ0-3au∗+ab+1+dc-d2c3u∗,a14=τ0d3c3u∗,a15=2τ0d2d-cc3u∗,a16=-τ0a+dc-dc2u∗2,a17=-τd4c4u∗2,a18=τ0c-d3d-cd2c4u∗2,a19=τ02c-3dd3c4u∗2,a23=-τ0dc-d2c2u∗,a24=-τ0d3c2u∗,a25=2τ0d2c-dc2u∗,a26=τ0d2c-d2c3u∗2,a27=τ0d4c3u∗2,a28=τ0d2c-dc-3dc3u∗2a29=τ0d33d-2cc3u∗2.The map (23) has a fixed point (0, 0) and the two eigenvalues at this fixed point can be written as(25)λτ0,λ-τ0=-pτ02±i^24qτ0-p2τ0=α±i^β,where i^=-1, pτ0 and qτ0 are described in (19) with τ=τ0. Based on conditions (21), the moduli of λτ0 and λ-τ0 satisfy(26)λτ0=qτ0=1.When Neimark-Sacker bifurcation occurs, the modulus of the eigenvalues (25) should be varying with the bifurcation parameter. This requires(27)d^=dλτdττ=τ0≠0.A direct calculation on d^ leads to(28)d^=dc-dcab+1u±-2ab-2c-dc.Since u±≠2abc+c-d/acb+1, then we find that the value of d^ 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, implying(29)λτ0θ≠1,θ=1,2,3,4,which is equivalent to pτ0≠-2,0,1,2. From (21), we find that p2τ0<4. It means pτ0≠-2,2, i.e.,(30)acb+1u±-2bdc-d+dc-d1-c2c2acb+1u±-2abc-2c-d-2c-2d≠0,1.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 (31)wz=a120α-a11-βw~z~on map (23), then we obtain(32)w~z~⟶α-ββαw~z~+1a12βG1w~,z~G2w~,z~,where(33)
G
1
w
~
,
z
~
=
a
12
2
a
13
β
+
a
14
β
α
-
a
11
2
+
a
12
a
15
β
α
-
a
11
w
~
2
+
a
14
β
3
z
~
2
-
2
a
14
β
2
α
-
a
11
+
a
12
a
15
β
2
w
~
z
~
+
a
12
3
a
16
β
+
a
17
β
α
-
a
11
3
+
a
12
2
a
18
β
α
-
a
11
+
a
12
a
19
β
α
-
a
11
2
w
~
3
-
a
17
β
4
z
~
3
+
3
a
17
β
3
α
-
a
11
+
a
12
a
19
β
3
w
~
z
~
2
-
3
a
17
β
2
α
-
a
11
2
+
a
12
2
a
18
β
2
+
2
a
12
a
19
β
2
α
-
a
11
w
~
2
z
~
+
O
ζ
4
,
(34)
G
2
w
~
,
z
~
=
a
12
2
a
13
α
-
a
11
-
a
12
a
23
+
a
14
α
-
a
11
-
a
12
a
24
α
-
a
11
2
+
a
12
a
15
α
-
a
11
-
a
12
2
a
25
α
-
a
11
w
~
2
+
β
2
a
14
α
-
a
11
-
a
12
a
24
z
~
2
-
2
β
a
14
α
-
a
11
-
a
12
a
24
α
-
a
11
+
β
a
12
a
15
α
-
a
11
-
a
12
2
a
25
w
~
z
~
+
a
12
3
a
16
α
-
a
11
-
a
12
4
a
26
+
a
17
α
-
a
11
-
a
12
a
27
α
-
a
11
3
+
a
12
2
a
18
α
-
a
11
-
a
12
3
a
28
α
-
a
11
+
a
12
a
19
α
-
a
11
-
a
12
2
a
29
α
-
a
11
2
w
~
3
-
β
3
a
17
α
-
a
11
-
a
12
a
27
z
~
3
+
β
2
3
α
-
a
11
a
17
α
-
a
11
-
a
12
a
27
+
a
12
a
19
α
-
a
11
-
a
12
2
a
29
w
~
z
~
2
-
β
3
α
-
a
11
2
a
17
α
-
a
11
-
a
12
a
27
+
a
12
2
a
18
α
-
a
11
-
a
12
3
a
28
+
2
α
-
a
11
a
12
a
19
α
-
a
11
-
a
12
2
a
29
w
~
2
z
~
+
O
ζ
4
.The second-order and third-order partial derivatives ofG1w~,z~ and G2w~,z~ at (0, 0) are calculated as(35)G1w~w~=2a122a13β+2a14βα-a11+2a12a15βα-a11,G1w~z~=-2a14β2α-a11-a12a15β3,G1z~z~=2a14β3,G1w~w~w~=6a123a16β+6a17βα-a113+6a122a18βα-a11+6a12a19βα-a112,G1w~w~z~=-6a17β2α-a112-2a122a18β2-4a12a19βα-a11,G1w~z~z~=6a17β3α-a11+2a12a19β3,G1z~z~z~=-6a17β4,G2w~w~=2a122a13α-a11-a12a23+2α-a112a14α-a11-a12a24+2α-a11a12a15α-a11-a122a25,G2w~z~=-2βa14α-a11-a12a24α-a11-βa12a15α-a11-a122a25,G2z~z~=2β2a14α-a11-a12a24,G2w~w~w~=6a123a16α-a11-a124a26+6a17α-a11-a12a27α-a113+6a122a18α-a11-a123a28α-a11+6a12a19α-a11-a122a29,G2w~w~z~=-6βa17α-a11-a12a27α-a112-2βa122a18α-a11-a123a28-4βa12a19α-a11-a122a29,G2w~z~z~=6β2α-a11a17α-a11-a12a27+2β2a12a19α-a11-a122a29,G2z~z~z~=-6β3a17α-a11-a12a27.According to the bifurcation theorem [33], if the map (32) undergoes Neimark-Sacker bifurcation at (u∗,v∗) around τ=τ0, the following discriminatory quantity E must not be zero, i.e.,(36)E=-Re1-2λ-λ-21-λξ11ξ20-12ξ112-ξ022+Reλ-ξ21≠0,in which(37)ξ20=18a12βG1w~w~-G1z~z~+2G2w~z~+i′G2w~w~-G2z~z~-2G1w~z~,ξ11=14a12βG1w~w~+G1z~z~+i′G2w~w~+G2z~z~,ξ02=18a12βG1w~w~-G1z~z~-2G2w~z~+i′G2w~w~-G2z~z~+2G1w~z~,ξ21=116a12βG1w~w~w~+G1w~z~z~+G2w~w~z~+G2z~z~z~+i′G2w~w~w~+G2w~z~z~-G1w~w~z~-G1z~z~z~.A calculation on the quantity E gives(38)E=-1AχMχG1w~w~+G1z~z~G1w~w~-G1z~z~+2G2w~z~-G2w~w~+G2z~z~G2w~w~-G2z~z~-2G1w~z~+NχG2w~w~+G2z~z~G1w~w~-G1z~z~+2G2w~z~-G1w~w~+G1z~z~G2w~w~-G2z~z~-2G1w~z~-132a122β2G1w~w~+G1z~z~2+G2w~w~+G2z~z~2+116a12βαG1w~w~w~+G1w~z~z~+G2w~w~z~+G2z~z~z~+βG2w~w~w~+G2w~z~z~-G1w~w~z~-G1z~z~z~-164a122β2G1w~w~-G1z~z~-2G2w~z~2+G2w~w~-G2z~z~+2G1w~z~2≠0,where(39)Aχ=32a122β21-α2+β2,Nχ=2αβ1-α1-2α-2β2-β3-4αα2-β2,Mχ=1-3α+2α2-2β2α2-β2+6α-8α2β2.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 E<0 and d^>0, then an attracting invariant closed curve bifurcates from the fixed point (u∗,v∗) for τ>τ0; otherwise, if E>0 and d^>0, a repelling invariant closed curve bifurcates from the fixed point for τ<τ0.


### 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 [10]. 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 (u∗,v∗), given by(40)u~i,j,t=ui,j,t-u∗,v~i,j,t=vi,j,t-v∗.It should be noticed that(41)Δdu~i,j,t=Δdui,j,t,Δdv~i,j,t=Δdvi,j,t.Based on the calculations in former approach [10], the eigenvalues of the operator Δd are determined as(42)λkl=4sin2⁡ϕk+sin2⁡ϕl,in which(43)ϕk=k-1πn,ϕl=l-1πn,k,l∈1,2,3,…,n.Substituting (40) and (41) into the spatiotemporal discrete predator-prey system leads to(44)u~i,j,t+1=a11u~i,j,t+τδ2D1Δdu~i,j,t+a12v~i,j,t+τδ2D2Δdv~i,j,t+Oς2,v~i,j,t+1=a21u~i,j,t+τδ2D1Δdu~i,j,t+a22v~i,j,t+τδ2D2Δdv~i,j,t+Oς2,where ς=u~i,j,t+v~i,j,t. The two-order terms Oς2 in (44) can be ignored when the perturbations are small. Using the corresponding eigenfunction Xklij of the eigenvalues λkl to multiply (44) gets(45)Xkliju~i,j,t+1=a11Xkliju~i,j,t+a12Xklijv~i,j,t+τδ2a11D1XklijΔdu~i,j,t+τδ2a12D2XklijΔdv~i,j,t,Xklijv~i,j,t+1=a21Xkliju~i,j,t+a22Xklijv~i,j,t+τδ2a21D1XklijΔdu~i,j,t+τδ2a22D2XklijΔdv~i,j,t.Summing (45) for all of i and j obtains(46)∑Xkliju~i,j,t+1=a11∑Xkliju~i,j,t+a12∑Xklijv~i,j,t+τδ2a11D1∑XklijΔdu~i,j,t+τδ2a12D2∑XklijΔdv~i,j,t,∑Xklijv~i,j,t+1=a21∑Xkliju~i,j,t+a22∑Xklijv~i,j,t+τδ2a21D1∑XklijΔdu~i,j,t+τδ2a22D2∑XklijΔdv~i,j,t.Let(47)u-t=∑Xkliju~i,j,t+1,v-t=∑Xklijv~i,j,t+1,Equations (47) can be transformed into the following form [10, 26]:(48)u-t+1=a111-τδ2D1λklu-t+a121-τδ2D2λklv-t,v-t+1=a211-τδ2D1λklu-t+a221-τδ2D2λklv-t.The Turing instability emerges when system (48) has unstable eigenvalues at the fixed point (0, 0) [10, 26]. Calculating these two eigenvalues of system (48), then we obtain(49)λ±k,l,τ=12Rk,l,τ±Rk,l,τ2-4Sk,l,τ,where(50)Rk,l,τ=-pτ-a11τD1+a22τD2τδ2λkl,Sk,l,τ=qτ1-D1τδ2λkl1-D2τδ2λkl.According to previous research works [10, 26–28], the occurrence of Turing instability demands that existing one group of kl to make the maximum value of λ±k,l larger than one. Define(51)Lτ=maxk,l=1n⁡max⁡λ+k,l,τ,λ-k,l,τ.Thus, when L(τ)>1, Turing instability occurs; when L(τ)<1, 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 L(τ)=1. Consequently, there exists a critical value τ=τ′ satisfying Lτ′=1.


## 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 [29]. 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 [32]. 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 [29].When map (11) undergoes Neimark-Sacker bifurcation, the stable fixed point (u∗,v∗) 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 [33]. Firstly, Neimark-Sacker bifurcation takes place if the eigenvalues in (18) are conjugate complex numbers and their moduli are one. That results in(20)p2-4q<0,q=1.Calculating on (19) yields the following conditions:(21)0<ab+1u±-2ab-c-d2c+d-cdc2<4τ,τ=τ0=cdc-d-1-cacb+1u±-2abc-2c-d.Under the satisfaction of conditions (21), the fixed point (u∗,v∗) of map (11) is then translated to the origin by the following translation:(22)w=u-u∗,z=v-v∗.Thus, the map (11) is transformed into a new map as(23)wz⟶a11w+a12z+a13w2+a14z2+a15wz+a16w3+a17w3+a18w2z+a19wz2+Oζ4a21w+a22z+a23w2+a24z2+a25wz+a26w3+a27w3+a28w2z+a29wz2+Oζ4,where ζ=w+z and Oζ4 describes a function with order at least four in the variables w,z, and the coefficients in (23) are described as follows:(24)a13=τ0-3au∗+ab+1+dc-d2c3u∗,a14=τ0d3c3u∗,a15=2τ0d2d-cc3u∗,a16=-τ0a+dc-dc2u∗2,a17=-τd4c4u∗2,a18=τ0c-d3d-cd2c4u∗2,a19=τ02c-3dd3c4u∗2,a23=-τ0dc-d2c2u∗,a24=-τ0d3c2u∗,a25=2τ0d2c-dc2u∗,a26=τ0d2c-d2c3u∗2,a27=τ0d4c3u∗2,a28=τ0d2c-dc-3dc3u∗2a29=τ0d33d-2cc3u∗2.The map (23) has a fixed point (0, 0) and the two eigenvalues at this fixed point can be written as(25)λτ0,λ-τ0=-pτ02±i^24qτ0-p2τ0=α±i^β,where i^=-1, pτ0 and qτ0 are described in (19) with τ=τ0. Based on conditions (21), the moduli of λτ0 and λ-τ0 satisfy(26)λτ0=qτ0=1.When Neimark-Sacker bifurcation occurs, the modulus of the eigenvalues (25) should be varying with the bifurcation parameter. This requires(27)d^=dλτdττ=τ0≠0.A direct calculation on d^ leads to(28)d^=dc-dcab+1u±-2ab-2c-dc.Since u±≠2abc+c-d/acb+1, then we find that the value of d^ 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, implying(29)λτ0θ≠1,θ=1,2,3,4,which is equivalent to pτ0≠-2,0,1,2. From (21), we find that p2τ0<4. It means pτ0≠-2,2, i.e.,(30)acb+1u±-2bdc-d+dc-d1-c2c2acb+1u±-2abc-2c-d-2c-2d≠0,1.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 (31)wz=a120α-a11-βw~z~on map (23), then we obtain(32)w~z~⟶α-ββαw~z~+1a12βG1w~,z~G2w~,z~,where(33)
G
1
w
~
,
z
~
=
a
12
2
a
13
β
+
a
14
β
α
-
a
11
2
+
a
12
a
15
β
α
-
a
11
w
~
2
+
a
14
β
3
z
~
2
-
2
a
14
β
2
α
-
a
11
+
a
12
a
15
β
2
w
~
z
~
+
a
12
3
a
16
β
+
a
17
β
α
-
a
11
3
+
a
12
2
a
18
β
α
-
a
11
+
a
12
a
19
β
α
-
a
11
2
w
~
3
-
a
17
β
4
z
~
3
+
3
a
17
β
3
α
-
a
11
+
a
12
a
19
β
3
w
~
z
~
2
-
3
a
17
β
2
α
-
a
11
2
+
a
12
2
a
18
β
2
+
2
a
12
a
19
β
2
α
-
a
11
w
~
2
z
~
+
O
ζ
4
,
(34)
G
2
w
~
,
z
~
=
a
12
2
a
13
α
-
a
11
-
a
12
a
23
+
a
14
α
-
a
11
-
a
12
a
24
α
-
a
11
2
+
a
12
a
15
α
-
a
11
-
a
12
2
a
25
α
-
a
11
w
~
2
+
β
2
a
14
α
-
a
11
-
a
12
a
24
z
~
2
-
2
β
a
14
α
-
a
11
-
a
12
a
24
α
-
a
11
+
β
a
12
a
15
α
-
a
11
-
a
12
2
a
25
w
~
z
~
+
a
12
3
a
16
α
-
a
11
-
a
12
4
a
26
+
a
17
α
-
a
11
-
a
12
a
27
α
-
a
11
3
+
a
12
2
a
18
α
-
a
11
-
a
12
3
a
28
α
-
a
11
+
a
12
a
19
α
-
a
11
-
a
12
2
a
29
α
-
a
11
2
w
~
3
-
β
3
a
17
α
-
a
11
-
a
12
a
27
z
~
3
+
β
2
3
α
-
a
11
a
17
α
-
a
11
-
a
12
a
27
+
a
12
a
19
α
-
a
11
-
a
12
2
a
29
w
~
z
~
2
-
β
3
α
-
a
11
2
a
17
α
-
a
11
-
a
12
a
27
+
a
12
2
a
18
α
-
a
11
-
a
12
3
a
28
+
2
α
-
a
11
a
12
a
19
α
-
a
11
-
a
12
2
a
29
w
~
2
z
~
+
O
ζ
4
.The second-order and third-order partial derivatives ofG1w~,z~ and G2w~,z~ at (0, 0) are calculated as(35)G1w~w~=2a122a13β+2a14βα-a11+2a12a15βα-a11,G1w~z~=-2a14β2α-a11-a12a15β3,G1z~z~=2a14β3,G1w~w~w~=6a123a16β+6a17βα-a113+6a122a18βα-a11+6a12a19βα-a112,G1w~w~z~=-6a17β2α-a112-2a122a18β2-4a12a19βα-a11,G1w~z~z~=6a17β3α-a11+2a12a19β3,G1z~z~z~=-6a17β4,G2w~w~=2a122a13α-a11-a12a23+2α-a112a14α-a11-a12a24+2α-a11a12a15α-a11-a122a25,G2w~z~=-2βa14α-a11-a12a24α-a11-βa12a15α-a11-a122a25,G2z~z~=2β2a14α-a11-a12a24,G2w~w~w~=6a123a16α-a11-a124a26+6a17α-a11-a12a27α-a113+6a122a18α-a11-a123a28α-a11+6a12a19α-a11-a122a29,G2w~w~z~=-6βa17α-a11-a12a27α-a112-2βa122a18α-a11-a123a28-4βa12a19α-a11-a122a29,G2w~z~z~=6β2α-a11a17α-a11-a12a27+2β2a12a19α-a11-a122a29,G2z~z~z~=-6β3a17α-a11-a12a27.According to the bifurcation theorem [33], if the map (32) undergoes Neimark-Sacker bifurcation at (u∗,v∗) around τ=τ0, the following discriminatory quantity E must not be zero, i.e.,(36)E=-Re1-2λ-λ-21-λξ11ξ20-12ξ112-ξ022+Reλ-ξ21≠0,in which(37)ξ20=18a12βG1w~w~-G1z~z~+2G2w~z~+i′G2w~w~-G2z~z~-2G1w~z~,ξ11=14a12βG1w~w~+G1z~z~+i′G2w~w~+G2z~z~,ξ02=18a12βG1w~w~-G1z~z~-2G2w~z~+i′G2w~w~-G2z~z~+2G1w~z~,ξ21=116a12βG1w~w~w~+G1w~z~z~+G2w~w~z~+G2z~z~z~+i′G2w~w~w~+G2w~z~z~-G1w~w~z~-G1z~z~z~.A calculation on the quantity E gives(38)E=-1AχMχG1w~w~+G1z~z~G1w~w~-G1z~z~+2G2w~z~-G2w~w~+G2z~z~G2w~w~-G2z~z~-2G1w~z~+NχG2w~w~+G2z~z~G1w~w~-G1z~z~+2G2w~z~-G1w~w~+G1z~z~G2w~w~-G2z~z~-2G1w~z~-132a122β2G1w~w~+G1z~z~2+G2w~w~+G2z~z~2+116a12βαG1w~w~w~+G1w~z~z~+G2w~w~z~+G2z~z~z~+βG2w~w~w~+G2w~z~z~-G1w~w~z~-G1z~z~z~-164a122β2G1w~w~-G1z~z~-2G2w~z~2+G2w~w~-G2z~z~+2G1w~z~2≠0,where(39)Aχ=32a122β21-α2+β2,Nχ=2αβ1-α1-2α-2β2-β3-4αα2-β2,Mχ=1-3α+2α2-2β2α2-β2+6α-8α2β2.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 E<0 and d^>0, then an attracting invariant closed curve bifurcates from the fixed point (u∗,v∗) for τ>τ0; otherwise, if E>0 and d^>0, a repelling invariant closed curve bifurcates from the fixed point for τ<τ0.


## 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 [10]. 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 (u∗,v∗), given by(40)u~i,j,t=ui,j,t-u∗,v~i,j,t=vi,j,t-v∗.It should be noticed that(41)Δdu~i,j,t=Δdui,j,t,Δdv~i,j,t=Δdvi,j,t.Based on the calculations in former approach [10], the eigenvalues of the operator Δd are determined as(42)λkl=4sin2⁡ϕk+sin2⁡ϕl,in which(43)ϕk=k-1πn,ϕl=l-1πn,k,l∈1,2,3,…,n.Substituting (40) and (41) into the spatiotemporal discrete predator-prey system leads to(44)u~i,j,t+1=a11u~i,j,t+τδ2D1Δdu~i,j,t+a12v~i,j,t+τδ2D2Δdv~i,j,t+Oς2,v~i,j,t+1=a21u~i,j,t+τδ2D1Δdu~i,j,t+a22v~i,j,t+τδ2D2Δdv~i,j,t+Oς2,where ς=u~i,j,t+v~i,j,t. The two-order terms Oς2 in (44) can be ignored when the perturbations are small. Using the corresponding eigenfunction Xklij of the eigenvalues λkl to multiply (44) gets(45)Xkliju~i,j,t+1=a11Xkliju~i,j,t+a12Xklijv~i,j,t+τδ2a11D1XklijΔdu~i,j,t+τδ2a12D2XklijΔdv~i,j,t,Xklijv~i,j,t+1=a21Xkliju~i,j,t+a22Xklijv~i,j,t+τδ2a21D1XklijΔdu~i,j,t+τδ2a22D2XklijΔdv~i,j,t.Summing (45) for all of i and j obtains(46)∑Xkliju~i,j,t+1=a11∑Xkliju~i,j,t+a12∑Xklijv~i,j,t+τδ2a11D1∑XklijΔdu~i,j,t+τδ2a12D2∑XklijΔdv~i,j,t,∑Xklijv~i,j,t+1=a21∑Xkliju~i,j,t+a22∑Xklijv~i,j,t+τδ2a21D1∑XklijΔdu~i,j,t+τδ2a22D2∑XklijΔdv~i,j,t.Let(47)u-t=∑Xkliju~i,j,t+1,v-t=∑Xklijv~i,j,t+1,Equations (47) can be transformed into the following form [10, 26]:(48)u-t+1=a111-τδ2D1λklu-t+a121-τδ2D2λklv-t,v-t+1=a211-τδ2D1λklu-t+a221-τδ2D2λklv-t.The Turing instability emerges when system (48) has unstable eigenvalues at the fixed point (0, 0) [10, 26]. Calculating these two eigenvalues of system (48), then we obtain(49)λ±k,l,τ=12Rk,l,τ±Rk,l,τ2-4Sk,l,τ,where(50)Rk,l,τ=-pτ-a11τD1+a22τD2τδ2λkl,Sk,l,τ=qτ1-D1τδ2λkl1-D2τδ2λkl.According to previous research works [10, 26–28], the occurrence of Turing instability demands that existing one group of kl to make the maximum value of λ±k,l larger than one. Define(51)Lτ=maxk,l=1n⁡max⁡λ+k,l,τ,λ-k,l,τ.Thus, when L(τ)>1, Turing instability occurs; when L(τ)<1, 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 L(τ)=1. Consequently, there exists a critical value τ=τ′ satisfying Lτ′=1.


## 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 [3], a group of ecologically feasible parameter values can be given as a = 1.63, b = 0.1, c = 0.32, d = 0.22, D1 = 0.6, D2 = 1.2, and n = 200. The values of parameters τ and δ should be set to ensure nonnegativity of the state variables u and v. In the following simulations, we let δ = 10 and τ as a varying parameter. Then according to previous research works [10, 27], the nonnegativity of u and v needs that τD1/δ2 and τD2/δ2 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τ0 = 4.9444. When τ<τ0, the fixed point (u∗,v∗) of the system is stable. Taking τ = 4.5 as an example, the two eigenvalues of the fixed point (u∗,v∗) are λ1 = 0.8307 + 0.5287i, λ2 = 0.8307 − 0.5287i, the modulus of which is less than one, suggesting (u∗,v∗) is a stable focus. When τ=τ0, we get the two eigenvalues of (u∗,v∗) as λ1 = 0.8140 + 0.5809i, λ2 = 0.8140 − 0.5809i, and |λ1,2| = 1. Then we check conditions (28), (30), and (38) for Neimark-Sacker bifurcation; we have d^ = 0.0152 ≠ 0, p(τ0) = −1.6280 ≠ 0,1 and E = −179.8092 ≠ 0. Hence, according to the Neimark-Sacker bifurcation theorem, the discrete predator-prey system indeed undergoes Neimark-Sacker bifurcation at τ=τ0. Furthermore, since E < 0 and d^ > 0, we know that an attracting invariant closed curve bifurcates from the fixed point (u∗,v∗) for τ>τ0, and the stable point (u∗,v∗) loses its stability.Figure1 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.Figure 1
(a) Neimark-Sacker bifurcation diagram and (b) maximum Lyapunov exponent corresponding to (a).a = 1.63,b = 0.1,c = 0.32, andd = 0.22.
(a) (b)Figure 2
Three typical dynamical behaviors on the route to chaos induced by the Neimark-Sacker bifurcation, (a) invariant closed curve, (b) period-31 orbit, and (c) chaotic attractor.a = 1.63, b = 0.1, c = 0.32, and d = 0.22.
(a)
τ = 5.1 (b)
τ = 5.525 (c)
τ = 5.83When 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 theL values corresponding to the three cases in Figure 2 are L(5.1) = 1.0060, L(5.525) = 1.0241, and L(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 (u∗, v∗) with small random perturbations. And deterministic initial condition is given by the following:(52)ui,j,0=1+0.01u∗,vi,j,0=1+0.01v∗,90≤i,j≤110,ui,j,0=u∗,vi,j,0=v∗,else  cases.Figures3–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.Figure 3
Pattern formation resulting from spatially symmetry breaking on invariant closed curves and starting with stochastic initial condition.a = 1.63, b = 0.1, c = 0.32, d = 0.22, D1 = 0.6, D2 = 1.2, τ = 5.1, δ = 10, and n = 200.
(a)
t = 1 (b)
t = 500 (c)
t = 1000 (d)
t = 1500 (e)
t = 2000 (f)
t = 5000Figure 4
Pattern formation resulting from spatially symmetry breaking on period-31 orbit and starting with stochastic initial condition.a = 1.63, b = 0.1, c = 0.32, d = 0.22, D1 = 0.6, D2 = 1.2, τ = 5.525, δ = 10, and n = 200.
(a)
t = 1 (b)
t = 10 (c)
t = 100 (d)
t = 500 (e)
t = 1000 (f)
t = 5000Figure 5
Pattern formation resulting from spatially symmetry breaking on chaotic attractor and starting with stochastic initial condition.a = 1.63, b = 0.1, c = 0.32, d = 0.22, D1 = 0.6, D2 = 1.2, τ = 5.83, δ = 10, and n = 200.
(a)
t = 1 (b)
t = 10 (c)
t = 100 (d)
t = 500 (e)
t = 2000 (f)
t = 5000As exhibited in Figures3–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, Figure6 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 u and v 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).Figure 6
Phase portraits and wave diagrams corresponding to the change of state variables in lattice (100, 100) of Figures3–5 with the increase of time.
(a) (b) (c) (d)
τ = 5.1 (e)
τ = 5.525 (f)
τ = 5.83To explore the nonlinear characteristics of pattern evolution process in Figures3–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, ui,j,0,vi,j,0 and u~i,j,0,v~i,j,0, the difference of which is given as ui,j,0-u~i,j,0=vi,j,0-v~i,j,0=0.0001 in merely one lattice, and then we observe the change of Pu=ui,j,t-u~i,j,t in another lattice. At the beginning, Pu stays at zero; after a certain time of pattern evolution, the amplitude of Pu suddenly bursts and Pu 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.Figure 7
Sensitive analysis of pattern evolution process to variations of initial conditions (graphs a–c) and parameterτ (graphs d–f), with the application of state variables changing process in Figure 6.
(a) (b) (c) (d)
τ = 5.1 (e)
τ = 5.525 (f)
τ = 5.83The above sensitivity analysis suggests that change of initial conditions may lead to the self-organization of different patterns. Figures8–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., t = 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).Figure 8
Pattern formation resulting from spatially symmetry breaking on invariant closed curves and starting with deterministic initial condition.a = 1.63, b = 0.1, c = 0.32, d = 0.22, D1 = 0.6, D2 = 1.2, τ = 5.1, δ = 10, and n = 200.
(a)
t = 1 (b)
t = 500 (c)
t = 2000 (d)
t = 3000 (e)
t = 5000 (f)
t = 7000 (g)
t = 10000 (h)
t = 100000 (i)
t = 200000Figure 9
Pattern formation resulting from spatially symmetry breaking on period-31 orbit and starting with deterministic initial condition.a = 1.63, b = 0.1, c = 0.32, d = 0.22, D1 = 0.6, D2 = 1.2, τ = 5.525, δ = 10, and n = 200.
(a)
t = 1 (b)
t = 500 (c)
t = 1000 (d)
t = 1500 (e)
t = 2000 (f)
t = 200000Figure 10
Pattern formation resulting from spatially symmetry breaking on chaotic attractor and starting with deterministic initial condition.a = 1.63, b = 0.1, c = 0.32, d = 0.22, D1 = 0.6, D2 = 1.2, τ = 5.83, δ = 10, and n = 200.
(a)
t = 1 (b)
t = 100 (c)
t = 500 (d)
t = 1000 (e)
t = 1500 (f)
t = 200000The decrease of symmetry in Figures8–10 as the time progresses results from the combined effects of Turing instability and Neimark-Sacker instability. A space-independent Neimark-Sacker bifurcation breaks the temporal symmetry of the system and gives rise to oscillations that are uniform in space and periodic in time. On the contrary, the diffusion-driven Turing instability breaks the spatial symmetry and manifests itself through the formation of patterns that are stationary in time and oscillatory in space. As explored in literature, the combined effects of both may lead to the emergence of spatiotemporal chaos, which presents disordered characteristics [10, 27]. The deterministic initial conditions provide order states which the predator-prey system starts and keeps at the beginning of the evolution process, but with the increase of time the spatiotemporal chaos will dominate and result in disordered states with lessening spatial symmetry.Figure11 shows the dynamical trajectories for the three cases of pattern evolution in Figures 8–10. Compared with Figure 6, we find that although the initial conditions are different, the attractors in phase portraits present similar configurations. Further sensitivity analysis demonstrates that the pattern evolution processes in these cases are also sensitive to tiny variation of parameter τ.Figure 11
Phase portraits corresponding to the change of state variables in lattice (100, 100) of Figures8–10 with the increase of time.
(a)
τ = 5.1 (b)
τ = 5.525 (c)
τ = 5.83The simulation results in Figures8–11 demonstrate that on the front part of the route to chaos (such as τ0 < τ < 5.2), patterns can still exhibit strongly symmetrical property at large evolution time (see t = 10000 in Figure 8), including horizontal, vertical, and central symmetry. Figure 12 displays a few transient strongly symmetrical patterns with the variation of parameter τ. These patterns show rich regular spatial structures. It reveals the spatiotemporal complexity of the discrete predator-prey system and the regularity in pattern self-organization.Figure 12
Various transient strongly symmetrical patterns when Neimark-Sacker-Turing instability occurs on the front part of route to chaos.a = 1.63, b = 0.1, c = 0.32, d = 0.22, D1 = 0.6, D2 = 1.2, δ = 10, n = 200, and t = 10000.
(a)
τ = 5.0 (b)
τ = 5.02 (c)
τ = 5.04 (d)
τ = 5.06 (e)
τ = 5.08 (f)
τ = 5.12 (g)
τ = 5.14 (h)
τ = 5.16 (i)
τ = 5.18


## 5. Discussion and Conclusions

This research investigates the spatiotemporal dynamics of a space- and time-discrete predator-prey system with Allee effect with the application of CML model. The development of the CML model in the present research is based on two types of discrete models in literature. The first is the type derived by directly discretizing partial differential equations with finite difference scheme, for example, the exponential discrete Lotka-Volterra model proposed by Li et al. [34]. The second type is the CML model as given by Mistro et al. [24], which abandoned the structure of partial differential equations. Taking the advantages of the two types of models, the present CML model is still developed from discretization of partial differential equations and simultaneously utilizes the structure of the CML model in Mistro et al. [24]. Therefore, the same ecological mechanisms with the corresponding continuous model as well as the advantages of CML model are both kept. The present CML model further bridges the continuous models and the discrete models [10].In the research, we find new nonlinear characteristics for the system dynamics and pattern formation under Neimark-Sacker-Turing instability. Summarizing the results obtained in theoretical analysis and numerical simulations, the concluding remarks of this research can be made as the following.(1)
The system can undergo Neimark-Sacker bifurcation inducing a route to chaos, on which three typical dynamical behaviors emerge, including invariant closed curves, periodic orbits, and chaotic attractors.(2)
Neimark-Sacker-Turing instability triggers spatial symmetry breaking on invariant closed curves, periodic orbits, and chaotic attractors and leads to the formation of heterogeneous patterns, which keep oscillating property all along.(3)
Different initial conditions result in diverse pattern evolution processes, revealing the multistability of the predator-prey system.(4)
Under the effect of Turing instability, the invariant closed curves and periodic orbits change to more complex chaotic attractors. Moreover, spatiotemporal chaos occurs when Neimark-Sacker instability and Turing instability both take place.(5)
With the application of deterministic initial conditions, transient symmetrical patterns are observed and the symmetry breaking increases as time continuously progresses. 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.(6)
Under Neimark-Sacker-Turing instability, the pattern evolution process is disordered and chaotic. However, at transient times, the predator-prey pattern may exhibit some symmetrical property. This demonstrates that ordered structures can exist in chaotic processes.(7)
Since pattern evolution processes are sensitive to variations of initial conditions and parameterτ, the system generates rich spatiotemporal complexity under the coupled effects of Neimark-Sacker instability and Turing instability. This reveals numerous possible ways for predator-prey interactions and coexistence.



---

*Source: 8713651-2018-06-26.xml*