# Bifurcations, Complex Behaviors, and Dynamic Transition in a Coupled Network of Discrete Predator-Prey System

Authors: Tousheng Huang; Huayong Zhang; Shengnan Ma; Ge Pan; Zhaodeng Wang; Hai Huang
Journal: Discrete Dynamics in Nature and Society
 (2019)

---

## Abstract

The nonlinear dynamics of predator-prey systems coupled into network is an important issue in recent biological advances. In this research, we consider each node of the coupled network represents a discrete predator-prey system, and the network dynamics is investigated. By applying Jacobian matrix, center manifold theorem and bifurcation theorems, stability of fixed points, flip bifurcation and Neimark-Sacker bifurcation of the discrete predator-prey system are analyzed. Via the method of Lyapunov exponents, the nonchaos-chaos transition of the coupled network along the routes to chaos induced by bifurcations is determined. Numerical simulations are performed to demonstrate the bifurcations, various attractors and dynamic transitions of the coupled network. Via comparison, we find that the coupled network exhibits far richer and more complex behaviors than single predator-prey system, including period-doubling cascades in orbits of period-2, period-4, period-8, invariant closed curves, dynamic windows for periodic orbits and invariant curves, quasiperiodic orbits, tori, and chaotic sets. Moreover, the attractors of the coupled network show more diverse and complicated structures. These results may provide a new perspective on the predator-prey dynamics in complex networks.


---

## Body

## 1. Introduction

Complex networks have been widely studied in the field of interdisciplinary, involving nonlinear science, complexity science, ecology, biology, and many other disciplines. The complex networks are topological abstractions of a large number of real complex systems, as well as topological foundations for the analysis of complexity. Via analyzing the characteristics of complex networks, many research works have played an important role in understanding origination and degree of complexity [1, 2]. Moreover, lots of complex networks investigations have contributed to solving practical problems, attracting attentions in the fields of science and engineering [3].Complex networks were originally within the scope of graph theory and can be traced back to the problem of the famous Seven Bridges of Konigsberg raised by Euler. The simplest network is regular network, which is an important basis of the complex networks. Three common regular networks studied in literature are globally coupled network, star network, and nearest-neighbor coupled network [4]. In 1950s, mathematicians Erdos and Renyi came up with an approach to construct random networks [5]. Regular networks and random networks are two types of extreme complex networks. Later, the researchers further constructed two important types of networks, small world network, and scale-free network [6, 7]. The small world networks have two statistical characteristics, namely, large family coefficients and small average distance [6]. In the scale-free network, node degree obeys power-law distribution [7]. The above networks reflect important characteristics of real complex networks and show great significance for exploring complex networks [8, 9].When each node in a complex network represents a nonlinear dynamic system, the complex network becomes a very complex dynamic system which demonstrates complicated nonlinear characteristics. The research on such dynamic complex networks mainly focused upon their diverse dynamic behaviors. For example, synchronization dynamics is one of the most interesting topics in research of dynamic complex networks [10, 11]. Wang et al. investigated the synchronization of coupled heterogeneous complex networks and derived sufficient criteria for quasi network synchronization [10]. The complex network control problems are attracted as well [12, 13]. Wang and Chen studied the control of a scale-free dynamical network by applying local feedback injections to a fraction of network nodes [13]. The researchers also concerned many other interesting dynamics of complex networks, such as the emergence of avalanche dynamics and cascading failures.Among the rich dynamics of complex networks, the transition of dynamics in coupled dynamical networks is a commonly concerned issue [4, 14]. Li et al. studied the transition to chaos in complex dynamical networks and verified that, during the transition from nonchaotic to chaotic states, if the network topology is more heterogeneous, the coupling strength required to achieve the transition can be decreased [14]. Zhang et al. researched the emergence of chaos in complex dynamical networks and found that the chaotic state will emerge without changing each node’s parameter if these nodes are connected through a certain type of network [4]. Moreover, a possible bifurcation route from chaos to hyperchaos was discovered in coupled dynamical networks [15]. How does network topology influence the dynamical transition of the complex networks has also become one focus of the research [14]. These approaches greatly improved the understanding of nonchaos-chaos transition in coupled dynamical networks.In recent decades, mathematical biology has developed to be a hot research field, and mathematical models have been established to analyze various biological phenomena. Due to ubiquity and importance of the predator-prey relationship, the dynamics of predator-prey systems has attracted the attentions of many scholars [16–20]. However, the dynamics of coupled networks of predator-prey systems is still seldom investigated and documented in literature.Early studies of predator-prey systems mainly utilized continuous models, via which the stability of equilibrium and existence and attractivity of limit cycles are often studied [21, 22]. When the populations have short life time and generations are nonoverlapping, it is more suitable applying discrete models to illustrate predator-prey systems [23]. Generally, discrete predator-prey models can often demonstrate more abundant dynamical behaviors than the continuous ones and therefore are particularly important for revealing nonlinear characteristics of discrete predator-prey systems [24, 25]. A great deal of research works manifested that the discrete predator-prey systems can exhibit various bifurcations and complex dynamic behaviors [16–20, 24, 25]. Through the research of He and Lai [16], Liu and Xiao [17], Jing and Yang [18], Agiza et al. [19], Hu et al. [20], Huang and Zhang [25], and etc., they found flip bifurcation and Neimark-Sacker bifurcation are basic bifurcations occurring in discrete predator-prey systems. The two bifurcations often start a route to chaos, revealing the dynamic transition from nonchaotic states to chaotic states. Moreover, complex dynamic behaviors can emerge on the routes to chaos, including cascades of period-doubling, periodic windows, quasiperiodic orbits and chaotic sets [17, 18, 25]. Via calculating the Lyapunov exponents and fractal dimension, Agiza et al. even determined diverse strange attractors of the predator-prey systems [19].The former research works on discrete predator-prey systems normally studied the dynamics of single system, but rarely took into account of the coupled networks of predator-prey system. Since natural predator-prey systems are usually not isolated and closely connected with each other, the coupled networks of predator-prey systems can be more accurate to reflect the characteristics of real ecological systems and deserve further investigations. On the other hand, the coupled networks of predator-prey systems may exhibit more complex dynamic behaviors. An investigation on coupled networks may provide new phenomena of predator-prey dynamics.In this research, we will study the bifurcations and dynamic transition between nonchaos and chaos in a coupled network of discrete predator-prey system. The exploration is arranged as follows. Section2 will give the model for the coupled dynamical network of discrete predator-prey system. Section 3 will analyze the stability of fixed points, period-doubling bifurcation, and Neimark-Sacker bifurcation when the discrete predator-prey system is not coupled. In Section 4, we will determine the sufficient conditions for emergence of chaotic state in the coupled dynamical network. In Section 5, numerical simulations will be carried out to demonstrate the complex dynamic behaviors. And finally, Section 6 will provide the conclusions.


## 2. Coupled Dynamical Network of Discrete Predator-Prey System

In this research, the complex dynamical network is considered to be modeled on the basis of a discrete Leslie-Gower predator-prey model. Leslie-Gower model is an important extension of Lotka-Volterra model. The Leslie-Gower predator-prey model assumes that interacting species grow according to the logistic law and that the environmental carrying capacity for the predator is not a constant but proportional to the population size of the prey [26]. However, due to the rarity of the prey, the predator can switch over to other food, but its growth is still limited by the fact that its favorite prey is not available in abundance [26].The predation relationship studied in the predator-prey system is considered as Crowley-Martin type functional response, which can accommodate interference among predators [27]. The Crowley-Martin type functional response is used for data sets that indicate feeding rate that is affected by predator density [28]. It is assumed that predator-feeding rate decreases by higher predator density even when prey density is high, and therefore the effects of predator interference in feeding rate remain important all the time whether an individual predator is handling or searching for a prey at a given instant of time [28].With Crowley-Martin type functional response considered, the continuous form of the Leslie-Gower predator-prey model can be described as follows:(1a)dudt=R1u1-uK-muv1+Au+Bv+ABuv,(1b)dvdt=R2v1-Dvu+E,in which u(t) and v(t) represent population densities of prey and predator at the time oft respectively;R1,R2,A,B,m,D,E, andK are positive constants;R1 and K are the intrinsic growth rate and the carrying capacity of the prey; the constants m, A, and B are the saturating Crowley-Martin type functional response parameters, where A measures the magnitude of interference among prey, B expresses the interference among the predator, and m represents the maximum consumption rate;R2 describes the growth rate of predator, D is the maximum value which per capita reduction rate of v can attain, and E measures the extent to which environment provides protection to the predator.In order to reduce the parameter numbers in the predator-prey model, we change the variables and parameters as follows:u→Ku, v→R1v/m,t→t/R1. Then (1a) and (1b) can be transformed into(2a)dudt=u1-u-v1+au1+bv,(2b)dvdt=cv1-dvu+e,where a=AK, b=BR1/m, c=R2/R1, d=DR1/mk, and e=E/K.Applying the forward Euler scheme to (2a) and (2b), the discrete-time Leslie-Gower predator-prey system can be expressed by the following equations:(3a)uk+1=fuk,vk=uk+τuk1-uk-vk1+auk1+bvk,(3b)vk+1=guk,vk=vk+cτvk1-dvkuk+e,in which k represents the sequence number of iterations and τ describes the time step size. It should be noticed that t(k+1)=t(k)+τ.Consider a coupled dynamical network ofN linearly and diffusively coupled identical nodes, with a full and diagonal coupling, where each node is a two-dimensional dynamical system described by (3a) and (3b). Therefore, the state equations of the coupled dynamical network of discrete predator-prey system are specified by [4, 14](4a)uik+1=fuik,vik-δ∑j=1Naijfujk,vjk,(4b)vik+1=guik,vik-δ∑j=1Naijgujk,vjk,where i and j denote the sequence number of the nodes in the coupled dynamical network, δ represents the coupling strength of the network, and Λ=(aij)N×N∈RN×N denotes the coupling configuration of the network. If there is a connection between node i and node j(i≠j), then aij=aji=1; otherwise, aij=aji=0(i≠j). Simultaneously, let aii=-di, i=1,2,…,N, where di is the degree of node i and can be defined to be the number of its outreaching connections, (5)aii=-∑j=1,j≠iNaij=-∑j=1,j≠iNaji,i=1,2,…,N,which means that the network is fully connected in the sense of having no isolated clusters.To investigate the dynamics of coupled dynamical network (4a) and (4b), we need to study the dynamics of one single node firstly, i.e., system (3a) and (3b). Via Jacobian determinant, bifurcation theorem, and center manifold theorem, stability analysis and bifurcation analysis are performed on system (3a) and (3b) to learn how the system dynamics responds to parameter variation. On the basis of investigations on system (3a) and (3b), then the dynamics of the coupled dynamical network is analyzed via applying the method of maximum Lyapunov exponent.


## 3. Stability Analysis and Bifurcation Analysis on the Discrete Predator-Prey System

### 3.1. Stability Analysis

For the convenience of analysis, we rewrite the expression of system (3a) and (3b) as the form of a map, i.e., (6)uv→u+τu1-u-v1+au1+bvv+cτv1-dvu+e.According to the definition of a fixed point for the map, the fixed points of map (6) can be calculated by solving the following equations:(7a)u=u+τu1-u-v1+au1+bv,(7b)v=v+cτv1-dvu+e.Direct calculation obtains four nonnegative fixed points, i.e.,(8)u0,v0:0,0,u1,v1:0,ed,u2,v2:1,0,u3,v3:u∗,v∗,where u∗ is a positive real root of the following cubic equation: (9)-abu3+ab-b-ad-abeu2+b+a-1d+a-1be-1u+d+be-e=0.If the inequalities are satisfied(10)e<eb+d,b+ad+abe>minab,d+eb+1,then the map (6) has an unique positive equilibrium (u∗,v∗), (11)u∗=-q′2+q′22+p′333+-q′2-q′22+p′333-c13,v∗=u∗+ed,p′=c2-c123,q′=c1272c12-9c2+c3,c1=b-ab+ad+abeab,c2=d-b-ad+be-abe+1ab,c3=e-d-beab.The stability of the fixed points is determined by the method of Jacobian matrix. The Jacobian matrix associated with any point of map (6) is (12)Ju,v=1+τ1-2u-v1+au21+bv-τu1+au1+bv2τcdv2u+e21+τc1-2dvu+e.The four fixed pointsu0,v0, u1,v1, u2,v2, and u3,v3 are substituted into the Jacobian matrix (12), respectively, and the eigenvalues λ1 and λ2 are calculated. If λ1<1 and λ2<1 are found, the corresponding fixed point is stable; if λ1>1 or λ2>1, an unstable fixed point emerges. The stability of the four fixed points can be determined as follows.(1)
The two eigenvalues ofJu0,v0 are λ1=1+τ and λ2=1+bτ. Obviously, we have λ1>1 and λ2>1, then we know u0,v0 is unstable.(2)
The two eigenvalues ofJu1,v1 are λ1=1+τ1-e/(d+be) and λ2=1-cτ. Hence, the conditions for that u1,v1 is stable can be determined as(13)e-d-be>0,0<τ<min2c,2d+bee-d-be.(3)
The two eigenvalues ofJu2,v2 are λ1=1-τ and λ2=1+cτ. Since λ2>1, the fixed point u2,v2 is also unstable.(4)
Foru3,v3, the corresponding Jacobian matrix is(14)Ju3,v3=a11a12a21a22,where(15)a11=1+τ1-2u∗-u∗+e1+au∗2d+bu∗+be,a12=-d2τu∗1+au∗d+bu∗+be2,a21=cτd,a22=1-cτ.Calculating the two eigenvalues ofJu3,v3 as (16)λ1,2=12-p±p2-4q,where(17a)p=-2+τc+2u∗+u∗+e1+au∗2d+bu∗+be-1,(17b)q=1+τ1-c-2u∗-u∗+e1+au∗2d+bu∗+be+cτ2du∗+e+bu∗+e2+du∗1+au∗1+au∗2d+bu∗+be2+2u∗-1.By calculatingλ1<1 and λ2<1, the stable conditions for (u3,v3) are determined as (18)q<1,-1+q<p<1+q.


### 3.2. Flip Bifurcation Analysis

Flip bifurcation is performed around the fixed pointu3,v3 when parameter variation occurs in the discrete predator-prey system. Parameter τ is chosen as the bifurcation parameter. According to the flip bifurcation theorem [29], the flip bifurcation occurs when one of the two eigenvalues of u3,v3 is equal to −1. In such case, we have q=p-1, and direct calculation yields(19)τ∗=min⁡4-y-y2-4x,4-y+y2-4x,where(20)x=cdu∗+e+bu∗+e2+du∗1+au∗1+au∗2d+bu∗+be2+2u∗-1,y=1-c-2u∗-u∗+e1+au∗2d+bu∗+be.When condition (19) establishes, the other eigenvalue of u3,v3 is obtained as λ2=1-p. For the occurrence of flip bifurcation, we need λ2≠1, namely, (21)τ∗c+2u∗+u∗+e1+au∗2d+bu∗+be-1≠2,4.Let(22)w=u-u∗,z=v-v∗,τ~=τ-τ∗,and we can translate the fixed point u3,v3 to the origin. Simultaneously, via doing Taylor expansion near the fixed point, the map (6) can be rewritten as (23)wzτ~→a11w+a12z+12a13w2+12a14z2+a15wz+a16wτ~+a17zτ~+16a18w3+12a19wz2+12a110w2z+16a111z3+12a112w2τ~+12a113z2τ~+a114wzτ~+ow+z+τ~4a21w+a22z+12a23w2+12a24z2+a25wz+a26wτ~+a27zτ~+16a28w3+12a29wz2+12a210w2z+12a211w2τ~+12a212z2τ~+a213wzτ~+ow+z+τ~4τ~,where a11, a12, a21, and a22 are described in (15) and(24)a13=τ-2+2au∗+e1+au∗3d+bu∗+be,a14=2τbd3u∗1+au∗d+bu∗+be3,a15=-d2τ1+au∗2d+bu∗+be2,a16=1-2u∗-u∗+e1+au∗2d+bu∗+be,a17=-d2u∗1+au∗d+bu∗+be2,a18=-6τa2u∗+e1+au∗4d+bu∗+be,a19=2bd3τ1+au∗2d+bu∗+be3,a110=2ad2τ1+au∗3d+bu∗+be2,a111=-6b2d4u∗τ1+au∗d+bu∗+be4,a112=-2+2au∗+e1+au∗3d+bu∗+be,a113=2bd3u∗1+au∗d+bu∗+be3,a114=-d21+au∗2d+bu∗+be2,a23=-2τcdu∗+e,a24=-2τcdu∗+e,a25=2τcu∗+e,a26=cd,a27=-c,a28=6τcdu∗+e2,a29=2τcdu∗+e2,a210=-4τcu∗+e2,a211=-2cdu∗+e,a212=-2cdu∗+e,a213=2cu∗+e.In (24), all theτ values satisfy τ=τ∗.On the basis of map (23), we make a reversible transformation as follows: (25)wzτ~→a12a120-1-a11λ2-a110001w~z~τ~.Then the map (23) is transformed to a norm form (26)w~z~τ~→-1000λ20001w~z~τ~+1a121+λ2F1w~,z~,τ~F2w~,z~,τ~0,in which (27a)
F
1
w
~
,
z
~
,
τ
~
=
a
12
2
2
a
13
λ
2
-
a
11
-
a
12
a
23
w
~
+
z
~
2
+
1
2
a
14
λ
2
-
a
11
-
a
12
a
24
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
+
a
12
a
15
λ
2
-
a
11
-
a
12
a
25
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
+
a
12
a
16
λ
2
-
a
11
-
a
12
a
26
w
~
+
z
~
τ
~
+
a
17
λ
2
-
a
11
-
a
12
a
27
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
τ
~
+
a
12
3
6
a
18
λ
2
-
a
11
-
a
12
a
28
w
~
+
z
~
3
+
a
12
2
a
19
λ
2
-
a
11
-
a
12
a
29
w
~
+
z
~
λ
2
-
a
11
z
~
+
1
+
a
11
w
~
2
+
a
12
2
2
a
110
λ
2
-
a
11
-
a
12
a
210
w
~
+
z
~
2
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
+
1
6
a
111
λ
2
-
a
11
λ
2
-
a
11
z
~
+
1
+
a
11
w
~
3
+
a
12
2
2
a
112
λ
2
-
a
11
-
a
12
a
211
w
~
+
z
~
2
τ
~
+
1
2
a
113
λ
2
-
a
11
-
a
12
a
212
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
τ
~
+
a
12
a
114
λ
2
-
a
11
-
a
12
a
213
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
τ
~
+
o
w
~
+
z
~
+
τ
~
4
,
(27b)
F
2
w
~
,
z
~
,
τ
~
=
a
12
2
2
a
13
1
+
a
11
+
a
12
a
23
w
~
+
z
~
2
+
1
2
a
14
1
+
a
11
+
a
12
a
24
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
+
a
12
a
15
1
+
a
11
+
a
12
a
25
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
+
a
12
a
16
1
+
a
11
+
a
12
a
26
w
~
+
z
~
τ
~
+
a
17
1
+
a
11
+
a
12
a
27
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
τ
~
+
a
12
3
6
a
18
1
+
a
11
+
a
12
a
28
w
~
+
z
~
3
+
a
12
2
a
19
1
+
a
11
+
a
12
a
29
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
+
a
12
2
2
a
110
1
+
a
11
+
a
12
a
210
w
~
+
z
~
2
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
+
1
6
a
111
1
+
a
11
λ
2
-
a
11
z
~
+
1
+
a
11
w
~
3
+
a
12
2
2
a
112
1
+
a
11
+
a
12
a
211
w
~
+
z
~
2
τ
~
+
a
113
1
+
a
11
2
+
a
12
a
212
2
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
τ
~
+
a
12
a
114
1
+
a
11
+
a
12
a
213
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
τ
~
+
o
w
~
+
z
~
+
τ
~
4
.Then we determine the central manifold of map (26) at the fixed point (0, 0, 0). By using the central manifold theorem [29], there is a central manifold WC(0,0,0) existing at the fixed point (0,0,0), which can be approximated by(28)WC0,0,0=w~,z~,τ~∈R3∣z~=h∗w~,τ~,h∗0,0=0,Dh∗0,0=0.In (28), h∗w~,τ~ is assumed to be(29)h∗w~,τ~=e0τ~+e1w~2+e2w~τ~+e3τ~2+ow~+τ~3.According to z~=h∗w~,τ~, we can get (30)e0τ~+e1-w~+F1w~,h∗w~,τ~,τ~a121+λ22+e2τ~-w~+F1w~,h∗w~,τ~,τ~a121+λ2+e3τ~2-λ2e0τ~+e1w~2+e2w~τ~+e3τ~2-F2w~,h∗w~,τ~,τ~a121+λ2=ow~+τ~3.Using the balance of (30) to determine the values of e0, e1, e2, and e3, we get the following results:(31)e0=0,e1=12a121-λ22a122a131+a11+a12a23+a141+a11+a12a241+a112-2a12a151+a11+a12a251+a11),e2=1a121+λ22a171+a112+a12a27-a161+a11-a122a26,e3=0.Consider the dynamics of map (26) restricted in the central manifold WC(0,0,0). Hence, a one-dimensional map is obtained, i.e., (32)F:w~→-w~+μ1w~2+μ2w~τ~+μ3w~2τ~+μ4w~τ~2+μ5w~3+ow~+τ~4.The coefficients in map (32) are expressed as follows:(33)μ1=1a121+λ2a1222a13λ2-a11-a12a23+12a14λ2-a11-a12a241+a112-a12a15λ2-a11-a12a251+a11,μ2=11+λ2a16λ2-a11-a12a26-1+a11a121+λ2a17λ2-a11-a12a27,μ3=121+λ2a12a112+2a12a13e2+2a16e1λ2-a11-a122a16e1+2a12a23e2+a12a211+2λ2-a111+a11a114+2e2λ2-2a11-1a15λ2-a11-a12a25+2a12a2131+a11+e2λ2-a111+a11a121+λ2a14λ2-a11-a12a24+e1λ2-a11a121+λ2a17λ2-a11-a12a27+1+a112a121+λ2a113λ2-a11-a12a212,μ4=e21+λ2a16λ2-a11-a12a26+e2λ2-a11a121+λ2a17λ2-a11-a12a27,μ5=161+λ2a12a12a18+6a13e1λ2-a11-a1226a23e1+a12a28+3a122a2101+a11-3a12a291+a112-3a12a110λ2-a111+a11+3a19λ2-a111+a112+6a14λ2-a11-a12a25λ2-2a11-1-e1λ2-a111+a11a121+λ2a14λ2-a11-a12a24+a111λ2-a111+a1136a121+λ2.If the map (32) experiences flip bifurcation, two determinant values, η1 and η2, should be nonzero, i.e., (34a)η1=∂2F∂w~∂τ~+12∂F∂τ~∂2F∂w~2≠0,w~=0,τ~=0,(34b)η2=16∂3F∂w~3+12∂2F∂w~22≠0,w~=0,τ~=0.It can be obtained by direct calculation, (35)η1=μ2≠0,η2=μ5+μ12≠0.According to the flip bifurcation theorem [29], the conditions of the flip bifurcation of the discrete predator-prey system are described as follows: if conditions (19) and (21) as well as (35) establish, the discrete predator-prey system undergoes flip bifurcation at the fixed point u3,v3 when the parameter τ varies in a small neighborhood of τ∗. Moreover, if η2>0, the period-2 orbit bifurcates from u3,v3 which is stable; if η2<0, the period-2 orbit bifurcates from u3,v3 which is unstable.


### 3.3. Neimark-Sacker Bifurcation Analysis

According to the Neimark-Sacker bifurcation theorem [29], the occurrence of Neimark-Sacker bifurcation first requires that the two eigenvalues in (16) are a pair of conjugate complex numbers whose module is one. This means p2-4q<0 and q=1, which lead to(36a)y2-4x<0,(36b)τ0=-yx.When parametric conditions satisfy (36a) and (36b), we translate the fixed point u3,v3 to the origin by the translation (37)w=u-u∗,z=v-v∗.Applying the Taylor expansion to map (6), we obtain a new map(38)wz→a11w+a12z+12a13w2+12a14z2+a15wz+16a18w3+12a19wz2+12a110w2z+16a111z3+ow+z4,a21w+a22z+12a23w2+12a24z2+a25wz+16a28w3+12a29wz2+12a210w2z+ow+z4.where the coefficients in map (38) are described in (24), but with τ=τ0. Under the conditions of (36a) and (36b), the two eigenvalues of the Jacobian matrix of map (38) at the fixed point (0, 0) are also a pair of conjugate complex with module as one. We here write these two eigenvalues as(39)λτ0,λ-τ0=-pτ02±i~24qτ0-p2τ0=α±i~β,where i~=-1, and explicitly, λ=qτ0=1.At the same time, the occurrence of Neimark-Sacker bifurcation needs that the derivation ofλ to parameter τ is not zero at τ=τ0. With a calculation, we have(40)γ=dλτdττ=τ0=121xτ2+yτ+12xτ+yτ=τ0=-y/x=-y2≠0,On the other hand, to avoid that λτ0 is a real number or a pure imaginary number, we also need that (41)λτ0n≠1,n=1,2,3,4.Condition (41) is equivalent to pτ0≠-2,0,1,2. Since pτ0<4qτ0=4, we know pτ0≠-2,2. Hence, condition (41) becomes -τ0y≠2,3, i.e., (42)τ0c+2u∗+u∗+e1+au∗2d+bu∗+be-1≠2,3.Then we study the norm form of map (38). With the invertible transformation(43)wz→a110α-a11-βw~z~,the map (38) can be transformed into the following new map:(44)w~z~→α-ββαw~z~+1a12βG1w~,z~G2w~,z~,where(45a)
G
1
w
~
,
z
~
=
a
12
2
a
13
2
β
+
a
14
2
β
α
-
a
11
2
+
a
12
a
15
α
-
a
11
β
w
~
2
+
a
14
2
β
3
z
2
-
a
14
α
-
a
11
β
2
+
a
12
a
15
β
2
w
~
z
~
+
a
12
3
a
18
6
+
a
12
a
19
2
β
α
-
a
11
2
+
a
12
2
a
110
2
β
α
-
a
11
+
1
6
a
111
β
α
-
a
11
3
w
~
3
+
a
12
a
19
2
β
3
+
a
111
2
β
3
α
-
a
11
w
~
z
~
2
-
a
12
a
19
α
-
a
11
β
2
+
a
12
2
a
110
2
β
2
+
a
111
2
α
-
a
11
2
β
2
w
~
2
z
~
-
a
111
6
β
4
z
~
3
+
o
w
~
+
z
~
4
,
(45b)
G
2
w
~
,
z
~
=
a
12
2
2
a
13
α
-
a
11
-
a
12
a
23
+
1
2
a
14
α
-
a
11
-
a
12
a
24
α
-
a
11
2
+
a
12
α
-
a
11
a
15
α
-
a
11
-
a
12
a
25
)
w
~
2
+
a
14
α
-
a
11
2
-
a
12
a
24
2
β
2
z
~
2
-
a
14
α
-
a
11
-
a
12
a
24
α
-
a
11
β
+
a
12
a
15
α
-
a
11
-
a
12
a
25
β
w
~
z
~
+
a
12
3
6
a
18
α
-
a
11
-
a
12
a
28
+
a
12
2
a
19
α
-
a
11
-
a
12
a
29
α
-
a
11
2
+
a
12
2
a
110
α
-
a
11
2
-
a
12
a
210
2
α
-
a
11
+
a
111
6
α
-
a
11
4
w
~
3
+
a
12
2
a
19
α
-
a
11
-
a
12
a
29
β
2
+
a
111
2
α
-
a
11
2
β
2
w
~
z
~
2
-
(
a
12
a
19
α
-
a
11
-
a
12
a
29
α
-
a
11
β
+
a
12
2
a
110
α
-
a
11
2
-
a
12
a
210
2
β
+
a
111
2
α
-
a
11
3
β
w
~
2
z
~
-
a
111
6
β
3
α
-
a
11
z
~
3
+
o
w
~
+
z
~
4
.Calculate the two- and three-order derivatives ofG1(w~,z~) and G2(w~,z~) at w~=0 and τ~=0. The following results are obtained: (46)G1w~w~=a122a132β+a142βα-a112+a12a15α-a11β,G1w~z~=-a14α-a11β2+a12a15β2,G1z~z~=a142β3,G1w~w~w~=a123a186+a12a192βα-a112+a122a1102βα-a11+16a111βα-a113,G1w~w~z~=-a12a19α-a11β2+a122a1102β2+a1112α-a112β2,G1w~z~z~=a12a192β3+a1112β3α-a11,G1z~z~z~=-a1116β4,G2w~w~=a1222a13α-a11-a12a23+12a14α-a11-a12a24α-a112+a12α-a11a15α-a11-a12a25,G2w~z~=-a14α-a11-a12a24α-a11β+a12a15α-a11-a12a25β,G2z~z~=a14α-a112-a12a242β2,G2w~w~w~=a123a18α-a116-a12a286+a12a19α-a112-a12a292α-a112+a122a110α-a112-a12a2102α-a11+a1116α-a114,G2w~w~z~=-a12a19α-a11-a12a29α-a11β+a122a110α-a112-a12a2102β+a1112α-a113β,G2w~z~z~=a122a19α-a11-a12a29β2+a1112α-a112β2,G2z~z~z~=-a1116β3α-a11.In order for map (44) experiencing Neimark-Sacker bifurcation at τ=τ0, the following determinant value must not be zero:(47)a0=-Re1-2λ-λ-21-λξ11ξ20-12ξ112-ξ022+Reλ-ξ21≠0,where (48)ξ20=18G1w~w~-G1z~z~+2G2w~z~+i~G2w~w~-G2z~z~-2G1w~z~,ξ11=14G1w~w~+G1z~z~+i~G2w~w~+G2z~z~,ξ02=18G1w~w~-G1z~z~-2G2w~z~+i~G2w~w~-G2z~z~+2G1w~z~,ξ21=116G1w~w~w~+G1w~z~z~+G2w~w~z~+G2z~z~z~+i~G2w~w~w~+G2w~z~z~-G1w~w~z~-G1z~z~z~.By using the above formulas, condition (47) becomes(49)a0=-1A0M0G1w~w~+G1z~z~G1w~w~-G1z~z~+2G2w~z~-G2w~w~+G2z~z~G2w~w~-G2z~z~-2G1w~z~+N0G2w~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-164a122β2G1w~w~-G1z~z~-2G2w~z~2+G2w~w~-G2z~z~+2G1w~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~≠0,in which (50)A0=32a122β21-α2+β2,N0=2αβ1-α1-2α-2β2-β3-4αα2-β2,M0=1-3α+2α2-2β2α2-β2+6α-8α2β2.Taking the above calculations together, the occurrence of Neimark-Sacker bifurcation at the fixed point(u3,v3) in the discrete predator-prey system needs the satisfaction of conditions (36a) and (36b), (42), and (49). Moreover, if a<0 and γ>0, then an attracting invariant closed curve bifurcates from the fixed point for τ>τ0; otherwise, if a>0 and γ>0, an repelling invariant closed curve bifurcates from the fixed point for τ<τ0.


## 3.1. Stability Analysis

For the convenience of analysis, we rewrite the expression of system (3a) and (3b) as the form of a map, i.e., (6)uv→u+τu1-u-v1+au1+bvv+cτv1-dvu+e.According to the definition of a fixed point for the map, the fixed points of map (6) can be calculated by solving the following equations:(7a)u=u+τu1-u-v1+au1+bv,(7b)v=v+cτv1-dvu+e.Direct calculation obtains four nonnegative fixed points, i.e.,(8)u0,v0:0,0,u1,v1:0,ed,u2,v2:1,0,u3,v3:u∗,v∗,where u∗ is a positive real root of the following cubic equation: (9)-abu3+ab-b-ad-abeu2+b+a-1d+a-1be-1u+d+be-e=0.If the inequalities are satisfied(10)e<eb+d,b+ad+abe>minab,d+eb+1,then the map (6) has an unique positive equilibrium (u∗,v∗), (11)u∗=-q′2+q′22+p′333+-q′2-q′22+p′333-c13,v∗=u∗+ed,p′=c2-c123,q′=c1272c12-9c2+c3,c1=b-ab+ad+abeab,c2=d-b-ad+be-abe+1ab,c3=e-d-beab.The stability of the fixed points is determined by the method of Jacobian matrix. The Jacobian matrix associated with any point of map (6) is (12)Ju,v=1+τ1-2u-v1+au21+bv-τu1+au1+bv2τcdv2u+e21+τc1-2dvu+e.The four fixed pointsu0,v0, u1,v1, u2,v2, and u3,v3 are substituted into the Jacobian matrix (12), respectively, and the eigenvalues λ1 and λ2 are calculated. If λ1<1 and λ2<1 are found, the corresponding fixed point is stable; if λ1>1 or λ2>1, an unstable fixed point emerges. The stability of the four fixed points can be determined as follows.(1)
The two eigenvalues ofJu0,v0 are λ1=1+τ and λ2=1+bτ. Obviously, we have λ1>1 and λ2>1, then we know u0,v0 is unstable.(2)
The two eigenvalues ofJu1,v1 are λ1=1+τ1-e/(d+be) and λ2=1-cτ. Hence, the conditions for that u1,v1 is stable can be determined as(13)e-d-be>0,0<τ<min2c,2d+bee-d-be.(3)
The two eigenvalues ofJu2,v2 are λ1=1-τ and λ2=1+cτ. Since λ2>1, the fixed point u2,v2 is also unstable.(4)
Foru3,v3, the corresponding Jacobian matrix is(14)Ju3,v3=a11a12a21a22,where(15)a11=1+τ1-2u∗-u∗+e1+au∗2d+bu∗+be,a12=-d2τu∗1+au∗d+bu∗+be2,a21=cτd,a22=1-cτ.Calculating the two eigenvalues ofJu3,v3 as (16)λ1,2=12-p±p2-4q,where(17a)p=-2+τc+2u∗+u∗+e1+au∗2d+bu∗+be-1,(17b)q=1+τ1-c-2u∗-u∗+e1+au∗2d+bu∗+be+cτ2du∗+e+bu∗+e2+du∗1+au∗1+au∗2d+bu∗+be2+2u∗-1.By calculatingλ1<1 and λ2<1, the stable conditions for (u3,v3) are determined as (18)q<1,-1+q<p<1+q.


## 3.2. Flip Bifurcation Analysis

Flip bifurcation is performed around the fixed pointu3,v3 when parameter variation occurs in the discrete predator-prey system. Parameter τ is chosen as the bifurcation parameter. According to the flip bifurcation theorem [29], the flip bifurcation occurs when one of the two eigenvalues of u3,v3 is equal to −1. In such case, we have q=p-1, and direct calculation yields(19)τ∗=min⁡4-y-y2-4x,4-y+y2-4x,where(20)x=cdu∗+e+bu∗+e2+du∗1+au∗1+au∗2d+bu∗+be2+2u∗-1,y=1-c-2u∗-u∗+e1+au∗2d+bu∗+be.When condition (19) establishes, the other eigenvalue of u3,v3 is obtained as λ2=1-p. For the occurrence of flip bifurcation, we need λ2≠1, namely, (21)τ∗c+2u∗+u∗+e1+au∗2d+bu∗+be-1≠2,4.Let(22)w=u-u∗,z=v-v∗,τ~=τ-τ∗,and we can translate the fixed point u3,v3 to the origin. Simultaneously, via doing Taylor expansion near the fixed point, the map (6) can be rewritten as (23)wzτ~→a11w+a12z+12a13w2+12a14z2+a15wz+a16wτ~+a17zτ~+16a18w3+12a19wz2+12a110w2z+16a111z3+12a112w2τ~+12a113z2τ~+a114wzτ~+ow+z+τ~4a21w+a22z+12a23w2+12a24z2+a25wz+a26wτ~+a27zτ~+16a28w3+12a29wz2+12a210w2z+12a211w2τ~+12a212z2τ~+a213wzτ~+ow+z+τ~4τ~,where a11, a12, a21, and a22 are described in (15) and(24)a13=τ-2+2au∗+e1+au∗3d+bu∗+be,a14=2τbd3u∗1+au∗d+bu∗+be3,a15=-d2τ1+au∗2d+bu∗+be2,a16=1-2u∗-u∗+e1+au∗2d+bu∗+be,a17=-d2u∗1+au∗d+bu∗+be2,a18=-6τa2u∗+e1+au∗4d+bu∗+be,a19=2bd3τ1+au∗2d+bu∗+be3,a110=2ad2τ1+au∗3d+bu∗+be2,a111=-6b2d4u∗τ1+au∗d+bu∗+be4,a112=-2+2au∗+e1+au∗3d+bu∗+be,a113=2bd3u∗1+au∗d+bu∗+be3,a114=-d21+au∗2d+bu∗+be2,a23=-2τcdu∗+e,a24=-2τcdu∗+e,a25=2τcu∗+e,a26=cd,a27=-c,a28=6τcdu∗+e2,a29=2τcdu∗+e2,a210=-4τcu∗+e2,a211=-2cdu∗+e,a212=-2cdu∗+e,a213=2cu∗+e.In (24), all theτ values satisfy τ=τ∗.On the basis of map (23), we make a reversible transformation as follows: (25)wzτ~→a12a120-1-a11λ2-a110001w~z~τ~.Then the map (23) is transformed to a norm form (26)w~z~τ~→-1000λ20001w~z~τ~+1a121+λ2F1w~,z~,τ~F2w~,z~,τ~0,in which (27a)
F
1
w
~
,
z
~
,
τ
~
=
a
12
2
2
a
13
λ
2
-
a
11
-
a
12
a
23
w
~
+
z
~
2
+
1
2
a
14
λ
2
-
a
11
-
a
12
a
24
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
+
a
12
a
15
λ
2
-
a
11
-
a
12
a
25
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
+
a
12
a
16
λ
2
-
a
11
-
a
12
a
26
w
~
+
z
~
τ
~
+
a
17
λ
2
-
a
11
-
a
12
a
27
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
τ
~
+
a
12
3
6
a
18
λ
2
-
a
11
-
a
12
a
28
w
~
+
z
~
3
+
a
12
2
a
19
λ
2
-
a
11
-
a
12
a
29
w
~
+
z
~
λ
2
-
a
11
z
~
+
1
+
a
11
w
~
2
+
a
12
2
2
a
110
λ
2
-
a
11
-
a
12
a
210
w
~
+
z
~
2
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
+
1
6
a
111
λ
2
-
a
11
λ
2
-
a
11
z
~
+
1
+
a
11
w
~
3
+
a
12
2
2
a
112
λ
2
-
a
11
-
a
12
a
211
w
~
+
z
~
2
τ
~
+
1
2
a
113
λ
2
-
a
11
-
a
12
a
212
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
τ
~
+
a
12
a
114
λ
2
-
a
11
-
a
12
a
213
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
τ
~
+
o
w
~
+
z
~
+
τ
~
4
,
(27b)
F
2
w
~
,
z
~
,
τ
~
=
a
12
2
2
a
13
1
+
a
11
+
a
12
a
23
w
~
+
z
~
2
+
1
2
a
14
1
+
a
11
+
a
12
a
24
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
+
a
12
a
15
1
+
a
11
+
a
12
a
25
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
+
a
12
a
16
1
+
a
11
+
a
12
a
26
w
~
+
z
~
τ
~
+
a
17
1
+
a
11
+
a
12
a
27
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
τ
~
+
a
12
3
6
a
18
1
+
a
11
+
a
12
a
28
w
~
+
z
~
3
+
a
12
2
a
19
1
+
a
11
+
a
12
a
29
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
+
a
12
2
2
a
110
1
+
a
11
+
a
12
a
210
w
~
+
z
~
2
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
+
1
6
a
111
1
+
a
11
λ
2
-
a
11
z
~
+
1
+
a
11
w
~
3
+
a
12
2
2
a
112
1
+
a
11
+
a
12
a
211
w
~
+
z
~
2
τ
~
+
a
113
1
+
a
11
2
+
a
12
a
212
2
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
τ
~
+
a
12
a
114
1
+
a
11
+
a
12
a
213
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
τ
~
+
o
w
~
+
z
~
+
τ
~
4
.Then we determine the central manifold of map (26) at the fixed point (0, 0, 0). By using the central manifold theorem [29], there is a central manifold WC(0,0,0) existing at the fixed point (0,0,0), which can be approximated by(28)WC0,0,0=w~,z~,τ~∈R3∣z~=h∗w~,τ~,h∗0,0=0,Dh∗0,0=0.In (28), h∗w~,τ~ is assumed to be(29)h∗w~,τ~=e0τ~+e1w~2+e2w~τ~+e3τ~2+ow~+τ~3.According to z~=h∗w~,τ~, we can get (30)e0τ~+e1-w~+F1w~,h∗w~,τ~,τ~a121+λ22+e2τ~-w~+F1w~,h∗w~,τ~,τ~a121+λ2+e3τ~2-λ2e0τ~+e1w~2+e2w~τ~+e3τ~2-F2w~,h∗w~,τ~,τ~a121+λ2=ow~+τ~3.Using the balance of (30) to determine the values of e0, e1, e2, and e3, we get the following results:(31)e0=0,e1=12a121-λ22a122a131+a11+a12a23+a141+a11+a12a241+a112-2a12a151+a11+a12a251+a11),e2=1a121+λ22a171+a112+a12a27-a161+a11-a122a26,e3=0.Consider the dynamics of map (26) restricted in the central manifold WC(0,0,0). Hence, a one-dimensional map is obtained, i.e., (32)F:w~→-w~+μ1w~2+μ2w~τ~+μ3w~2τ~+μ4w~τ~2+μ5w~3+ow~+τ~4.The coefficients in map (32) are expressed as follows:(33)μ1=1a121+λ2a1222a13λ2-a11-a12a23+12a14λ2-a11-a12a241+a112-a12a15λ2-a11-a12a251+a11,μ2=11+λ2a16λ2-a11-a12a26-1+a11a121+λ2a17λ2-a11-a12a27,μ3=121+λ2a12a112+2a12a13e2+2a16e1λ2-a11-a122a16e1+2a12a23e2+a12a211+2λ2-a111+a11a114+2e2λ2-2a11-1a15λ2-a11-a12a25+2a12a2131+a11+e2λ2-a111+a11a121+λ2a14λ2-a11-a12a24+e1λ2-a11a121+λ2a17λ2-a11-a12a27+1+a112a121+λ2a113λ2-a11-a12a212,μ4=e21+λ2a16λ2-a11-a12a26+e2λ2-a11a121+λ2a17λ2-a11-a12a27,μ5=161+λ2a12a12a18+6a13e1λ2-a11-a1226a23e1+a12a28+3a122a2101+a11-3a12a291+a112-3a12a110λ2-a111+a11+3a19λ2-a111+a112+6a14λ2-a11-a12a25λ2-2a11-1-e1λ2-a111+a11a121+λ2a14λ2-a11-a12a24+a111λ2-a111+a1136a121+λ2.If the map (32) experiences flip bifurcation, two determinant values, η1 and η2, should be nonzero, i.e., (34a)η1=∂2F∂w~∂τ~+12∂F∂τ~∂2F∂w~2≠0,w~=0,τ~=0,(34b)η2=16∂3F∂w~3+12∂2F∂w~22≠0,w~=0,τ~=0.It can be obtained by direct calculation, (35)η1=μ2≠0,η2=μ5+μ12≠0.According to the flip bifurcation theorem [29], the conditions of the flip bifurcation of the discrete predator-prey system are described as follows: if conditions (19) and (21) as well as (35) establish, the discrete predator-prey system undergoes flip bifurcation at the fixed point u3,v3 when the parameter τ varies in a small neighborhood of τ∗. Moreover, if η2>0, the period-2 orbit bifurcates from u3,v3 which is stable; if η2<0, the period-2 orbit bifurcates from u3,v3 which is unstable.


## 3.3. Neimark-Sacker Bifurcation Analysis

According to the Neimark-Sacker bifurcation theorem [29], the occurrence of Neimark-Sacker bifurcation first requires that the two eigenvalues in (16) are a pair of conjugate complex numbers whose module is one. This means p2-4q<0 and q=1, which lead to(36a)y2-4x<0,(36b)τ0=-yx.When parametric conditions satisfy (36a) and (36b), we translate the fixed point u3,v3 to the origin by the translation (37)w=u-u∗,z=v-v∗.Applying the Taylor expansion to map (6), we obtain a new map(38)wz→a11w+a12z+12a13w2+12a14z2+a15wz+16a18w3+12a19wz2+12a110w2z+16a111z3+ow+z4,a21w+a22z+12a23w2+12a24z2+a25wz+16a28w3+12a29wz2+12a210w2z+ow+z4.where the coefficients in map (38) are described in (24), but with τ=τ0. Under the conditions of (36a) and (36b), the two eigenvalues of the Jacobian matrix of map (38) at the fixed point (0, 0) are also a pair of conjugate complex with module as one. We here write these two eigenvalues as(39)λτ0,λ-τ0=-pτ02±i~24qτ0-p2τ0=α±i~β,where i~=-1, and explicitly, λ=qτ0=1.At the same time, the occurrence of Neimark-Sacker bifurcation needs that the derivation ofλ to parameter τ is not zero at τ=τ0. With a calculation, we have(40)γ=dλτdττ=τ0=121xτ2+yτ+12xτ+yτ=τ0=-y/x=-y2≠0,On the other hand, to avoid that λτ0 is a real number or a pure imaginary number, we also need that (41)λτ0n≠1,n=1,2,3,4.Condition (41) is equivalent to pτ0≠-2,0,1,2. Since pτ0<4qτ0=4, we know pτ0≠-2,2. Hence, condition (41) becomes -τ0y≠2,3, i.e., (42)τ0c+2u∗+u∗+e1+au∗2d+bu∗+be-1≠2,3.Then we study the norm form of map (38). With the invertible transformation(43)wz→a110α-a11-βw~z~,the map (38) can be transformed into the following new map:(44)w~z~→α-ββαw~z~+1a12βG1w~,z~G2w~,z~,where(45a)
G
1
w
~
,
z
~
=
a
12
2
a
13
2
β
+
a
14
2
β
α
-
a
11
2
+
a
12
a
15
α
-
a
11
β
w
~
2
+
a
14
2
β
3
z
2
-
a
14
α
-
a
11
β
2
+
a
12
a
15
β
2
w
~
z
~
+
a
12
3
a
18
6
+
a
12
a
19
2
β
α
-
a
11
2
+
a
12
2
a
110
2
β
α
-
a
11
+
1
6
a
111
β
α
-
a
11
3
w
~
3
+
a
12
a
19
2
β
3
+
a
111
2
β
3
α
-
a
11
w
~
z
~
2
-
a
12
a
19
α
-
a
11
β
2
+
a
12
2
a
110
2
β
2
+
a
111
2
α
-
a
11
2
β
2
w
~
2
z
~
-
a
111
6
β
4
z
~
3
+
o
w
~
+
z
~
4
,
(45b)
G
2
w
~
,
z
~
=
a
12
2
2
a
13
α
-
a
11
-
a
12
a
23
+
1
2
a
14
α
-
a
11
-
a
12
a
24
α
-
a
11
2
+
a
12
α
-
a
11
a
15
α
-
a
11
-
a
12
a
25
)
w
~
2
+
a
14
α
-
a
11
2
-
a
12
a
24
2
β
2
z
~
2
-
a
14
α
-
a
11
-
a
12
a
24
α
-
a
11
β
+
a
12
a
15
α
-
a
11
-
a
12
a
25
β
w
~
z
~
+
a
12
3
6
a
18
α
-
a
11
-
a
12
a
28
+
a
12
2
a
19
α
-
a
11
-
a
12
a
29
α
-
a
11
2
+
a
12
2
a
110
α
-
a
11
2
-
a
12
a
210
2
α
-
a
11
+
a
111
6
α
-
a
11
4
w
~
3
+
a
12
2
a
19
α
-
a
11
-
a
12
a
29
β
2
+
a
111
2
α
-
a
11
2
β
2
w
~
z
~
2
-
(
a
12
a
19
α
-
a
11
-
a
12
a
29
α
-
a
11
β
+
a
12
2
a
110
α
-
a
11
2
-
a
12
a
210
2
β
+
a
111
2
α
-
a
11
3
β
w
~
2
z
~
-
a
111
6
β
3
α
-
a
11
z
~
3
+
o
w
~
+
z
~
4
.Calculate the two- and three-order derivatives ofG1(w~,z~) and G2(w~,z~) at w~=0 and τ~=0. The following results are obtained: (46)G1w~w~=a122a132β+a142βα-a112+a12a15α-a11β,G1w~z~=-a14α-a11β2+a12a15β2,G1z~z~=a142β3,G1w~w~w~=a123a186+a12a192βα-a112+a122a1102βα-a11+16a111βα-a113,G1w~w~z~=-a12a19α-a11β2+a122a1102β2+a1112α-a112β2,G1w~z~z~=a12a192β3+a1112β3α-a11,G1z~z~z~=-a1116β4,G2w~w~=a1222a13α-a11-a12a23+12a14α-a11-a12a24α-a112+a12α-a11a15α-a11-a12a25,G2w~z~=-a14α-a11-a12a24α-a11β+a12a15α-a11-a12a25β,G2z~z~=a14α-a112-a12a242β2,G2w~w~w~=a123a18α-a116-a12a286+a12a19α-a112-a12a292α-a112+a122a110α-a112-a12a2102α-a11+a1116α-a114,G2w~w~z~=-a12a19α-a11-a12a29α-a11β+a122a110α-a112-a12a2102β+a1112α-a113β,G2w~z~z~=a122a19α-a11-a12a29β2+a1112α-a112β2,G2z~z~z~=-a1116β3α-a11.In order for map (44) experiencing Neimark-Sacker bifurcation at τ=τ0, the following determinant value must not be zero:(47)a0=-Re1-2λ-λ-21-λξ11ξ20-12ξ112-ξ022+Reλ-ξ21≠0,where (48)ξ20=18G1w~w~-G1z~z~+2G2w~z~+i~G2w~w~-G2z~z~-2G1w~z~,ξ11=14G1w~w~+G1z~z~+i~G2w~w~+G2z~z~,ξ02=18G1w~w~-G1z~z~-2G2w~z~+i~G2w~w~-G2z~z~+2G1w~z~,ξ21=116G1w~w~w~+G1w~z~z~+G2w~w~z~+G2z~z~z~+i~G2w~w~w~+G2w~z~z~-G1w~w~z~-G1z~z~z~.By using the above formulas, condition (47) becomes(49)a0=-1A0M0G1w~w~+G1z~z~G1w~w~-G1z~z~+2G2w~z~-G2w~w~+G2z~z~G2w~w~-G2z~z~-2G1w~z~+N0G2w~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-164a122β2G1w~w~-G1z~z~-2G2w~z~2+G2w~w~-G2z~z~+2G1w~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~≠0,in which (50)A0=32a122β21-α2+β2,N0=2αβ1-α1-2α-2β2-β3-4αα2-β2,M0=1-3α+2α2-2β2α2-β2+6α-8α2β2.Taking the above calculations together, the occurrence of Neimark-Sacker bifurcation at the fixed point(u3,v3) in the discrete predator-prey system needs the satisfaction of conditions (36a) and (36b), (42), and (49). Moreover, if a<0 and γ>0, then an attracting invariant closed curve bifurcates from the fixed point for τ>τ0; otherwise, if a>0 and γ>0, an repelling invariant closed curve bifurcates from the fixed point for τ<τ0.


## 4. Analysis on Nonchaos-Chaos Transition in the Coupled Dynamical Network

For studying the dynamic transition between nonchaos and chaos in the coupled dynamical network, Lyapunov exponent is an important tool. As well known, the Lyapunov exponent can determine whether chaos exists or not. If the maximum Lyapunov exponent is positive, chaotic dynamics exists in nonlinear systems. In this research, the analyzing method as described in Zhang et al. [4] and Li et al. [14] is applied. Via studying the eigenvalues of coupling matrix of the coupled dynamics network, the positivity and negativity of the maximum Lyapunov exponent can be determined. According to the relationship between maximum Lyapunov exponent and chaos, the conditions for nonchaotic and chaotic network dynamics can be found.Before analyzing the dynamics of the coupled dynamical network, we need firstly describe the properties of the coupling matrixΛ=(aij)N×N. According to the research and description in Wu and Chua [30] and Zhang et al. [4], if the coupling matrix is an irreducible matrix and satisfies the conditions defined in Section 2, then Λ=(aij)N×N has the following properties:(1) 0 is the eigenvalue of the matrix Λ, and the corresponding eigenvector is(1,1,…,1)T, where T means the transpose of a matrix.(2) All other eigenvalues of the matrix Λ are real and strictly negative.On the basis of the above properties, we define and rank theN eigenvalues of the coupling matrix Λ as the following: (51)0=λ1>λ2≥⋯≥λN.For convenience of analysis, we rewrite the equations of the coupled dynamical network (see (4a) and (4b)) in a vector form, (52)sk+1=Fsk,δ,in which s(k)=[s1(k),s2(k),…,sN(k)]T and si(k)=[ui(k),vi(k)](i=1,2,…,N). Generally, the Lyapunov exponents of the coupled dynamical network can be calculated as(53)γi=limr→+∞1rln⁡DFrs0·ξi,i=1,2,…,N,where DFr(s0) is the Jacobian matrix of the r iterated map starting from a random initial state s0 and ξi(i=1,2,…,N) is a set of orthonormal vectors in the tangent space of the map (52).We denote the maximum Lyapunov exponent of the discrete predator-prey system (3a) and (3b) as h0, which can be described as [31](54)h0=limr→+∞1r∑k=0r-1ln⁡DF0sk,where F0=(f,g), i.e., the functions defined in system (3a) and (3b). Choose parameter τ as varying parameter. Then from (54), we know that there exists a relationship between h0 and τ when the τ value keeps in variation. We here express this relationship as the following: (55)h0=φτ.Generally, φ can be determined by numerical method.Furthermore, according to Zhang et al. [4] and Li et al. [14], there is also a relation between h0 and γi(i=1,2,…,N) which can be described as (56)γiλi=h0+ln⁡1-δλi,i=1,2,…,N,where λi is the ith eigenvalue of the coupling matrix Λ as described in (51). Because the coupling strength δ is positive, then the Lyapunov exponents of the coupled dynamical network can be ranked as follows:(57)γNλN=h0+ln⁡1-δλN≥γN-1λN-1≥…≥γ2λ2≥γ1λ1=h0.The purpose here is to investigate the aggregative dynamics of the coupled dynamical network and determine the nonchaotic and chaotic regions with the variation parameterτ. In this research, we consider the case that the parameters of all coupling nodes are the same; i.e., the dynamics of the discrete predator-prey system at any node is attracted to the same attractor if the initial conditions are provided the same. Then, if the coupled dynamical network exhibits chaotic dynamics, there is at least one positive Lyapunov exponent in (57), implying γN>0. On this basis, we can classify the Lyapunov exponents of the coupled dynamical network into two types, and each type of Lyapunov exponents satisfies the following: (58a)γN≥γN-1≥⋯≥γT+1>0,(58b)h0=γ1<γ2≤⋯≤γT<0.Here, T(1≤T≤N-1) is a positive integer and γT is the largest negative Lyapunov exponent.Combining (53)-(58a) and (58b), we can obtain the condition where at least one positive Lyapunov exponent emerges in the coupled dynamical network, i.e., (59)φτ>ln⁡11-δλN.Via solving inequality (59), the range of parameterτ in which the coupled dynamical network exhibits chaotic dynamics is determined. Therefore, with the increase ofτvalue, the distribution of nonchaotic and chaotic regions is known. Combining with the bifurcation analysis of the discrete predator-prey system, we can investigate the dynamic transition between nonchaos and chaos in the coupled network when flip bifurcation or Neimark-Sacker bifurcation take place. In the following numerical simulations, we will also demonstrate how the coupled dynamical network influences the bifurcations of the discrete predator-prey system.


## 5. Numerical Simulations

In this section, numerical simulations are performed to verify the previous theoretical results and exhibit complex and rich dynamical behaviors of the discrete predator-prey system and the coupled dynamical network. Bifurcation diagrams, phase portraits, and maximum Lyapunov exponents are plotted. Via comparing the dynamics of single predator-prey system and the coupled network, how the coupled network influences the predator-prey dynamics is investigated.


### 5.1. Bifurcations and Complex Dynamics of the Discrete Predator-Prey System

To exhibit the flip bifurcation of the discrete predator-prey system, we give a group of parameter values asa=0.1, b=0.001, c=10, d=0.05, and e=0.037 and choose parameter τ as bifurcation parameter. At the beginning, the occurrence of the flip bifurcation under such parametric conditions is demonstrated. The fixed point of the discrete predator-prey system is determined as (0.0125, 0.9897), and the critical point for flip bifurcation is calculated as τ∗=0.2052. When τ<τ∗, the fixed point (u3,v3) of the system is stable. Taking τ=0.2 for example, the two eigenvalues of the fixed point (u3,v3) are λ1=-0.9489, λ2=0.9466, the modulus of which are both less than one, suggesting (u3,v3) is a stable node. When τ=τ∗, we get the two eigenvalues of (u3,v3) as λ1=-1, λ2=0.9452. Then map (26) restricted to the center manifold is got as F:w~→-w~+4.0288w~2-9.7433w~τ~+57.8564w~2τ~-15.1681μ5w~3+o((w~+τ~)4). On the basis of the one-dimensional map F, the two discriminatory quantities η1 and η2 are got as η1=-9.7443, η2=1.0635. Hence, according to the flip bifurcation theorem, the discrete predator-prey system indeed undergoes flip bifurcation at τ=τ∗. Furthermore, because η2>0, we know that a period-2 orbit bifurcating from (u3,v3) is stable and the fixed point (u3,v3) becomes unstable. Taking τ=0.21, for example, the two eigenvalues of (u3,v3) are λ1=-1.0463, λ2=0.9439, respectively, therefore (u3,v3) is a saddle point. Simultaneously, the discrete predator-prey system is attracted to the two stable periodic points (0.0128, 0.8841) and (0.0131, 1.0923).Figure1(a) shows the flip bifurcation diagram of the discrete predator-prey system with parameter τ ranging in [0.2,0.3] under the above given parametric conditions. Explicitly, a period-doubling cascade is observed and the dynamic transition from periodic orbits to chaotic behaviors is found. Figures 1(b) and 1(c) are local amplifications of Figure 1(a), demonstrating the detail of periodic windows with period-6 and period-5 orbits. In the periodic windows, period-doubling cascade also occurs, resulting to period-10, period-20, period-40 orbits and period-12, period-24, period-48 orbits, etc. Moreover, other periodic windows can also be found in the chaotic region. Figure 2 shows the maximum Lyapunov exponents corresponding to the three bifurcation diagrams in Figure 1. It can be seen in Figure 2 that the discrete predator-prey system first goes into chaotic region at about τ=0.2642.Figure 1
(a) Flip bifurcation diagram of the discrete predator-prey system corresponding to variation ofτ ranging in [0.2,0.3]. a=0.1, b=0.001, c=10, d=0.05, and e=0.037. (b) and (c) Local amplifications of (a) in [0.27,0.2712] and [0.2812,0.2822] showing the detail of periodic windows.
(a) (b) (c)Figure 2
Diagrams of maximum Lyapunov exponent corresponding to the three graphs of Figure1.
(a) (b) (c)Figure3 shows different attractors in the phase portraits corresponding to the flip bifurcation diagram in Figure 1(a). Figures 3(a)–3(d) display the period-1, period-2, period-4, period-8 orbits. The doubling of periodic points in these orbits coincides with the period-doubling process in Figure 1(a). Figures 3(e), 3(g), and 3(i) show a cloud of random points, indicating chaotic behaviors. The maximum Lyapunov exponents of these dynamical behaviors are 0.0756, 0.1334, and 0.3265, respectively, confirming that these behaviors are indeed chaotic. Figures 3(f) and 3(h) display the period-6 and period-5 orbits; in fact, they correspond to the periodic windows as shown in Figures 1(b) and 1(c). The dynamic transition from period to chaos demonstrated in Figure 3 suggests a disordering process occurring in the discrete predator-prey system.Figure 3
Phase diagrams for different values of parameterτ corresponding to Figure 1(a).
(a)
τ = 0.2 (b)
τ = 0.24 (c)
τ = 0.26 (d)
τ = 0.263 (e)
τ = 0.267 (f)
τ = 0.2703 (g)
τ =0.2715 (h)
τ = 0.2814 (i)
τ = 0.295Figures4–6 demonstrate the dynamics of Neimark-Sacker bifurcation of the discrete predator-prey system. Before describing these figures, an example is provided to illustrate the occurrence of Neimark-Sacker bifurcation. Given parameter values as a=0.1, b=1.5, c=0.5, d=3.5, and e=0.05, the fixed point (u3,v3) is obtained as (0.8313, 0.2518), and the critical point for Neimark-Sacker bifurcation is determined as τ0=2.8234. When τ<τ0, taking τ=2.6, for instance, we get the two eigenvalues of (u3,v3) are λ1,2=-0.7138±0.4682i~ and λ1,2=0.8537; therefore (u3,v3) is a stable focus. When we take τ=τ0, the two eigenvalues are λ1=-0.8611±0.5085i~ and λ1,2=1. Besides, through calculation, we have γ=0.6592 and a=-16.8134. It is verified that the discrete predator-prey system undergoes Neimark-Sacker bifurcation at τ=τ0. Moreover, on account of γ>0 and a<0, we have an attracting invariant closed curve bifurcating from the fixed point for τ>τ0.Figure 4
(a) Neimark-Sacker bifurcation diagram of the discrete predator-prey system corresponding to variation ofτ ranging in [2.8,3.2]. a=0.1, b=1.5, c=0.5, d=3.5, and e=0.05. (b) and (c) Local amplifications of (a) in [3.045,3.07] and [3.065,3.11] showing the detail of periodic windows.
(a) (b) (c)Figure 5
Diagrams of maximum Lyapunov exponent corresponding to the three graphs of Figure4.
(a) (b) (c)Figure 6
Phase diagrams for different values of parameterτ corresponding to Figure 4(a).
(a)
τ= 2.815 (b)
τ= 2.9 (c)
τ= 3.015 (d)
τ= 3.037 (e)
τ= 3.06 (f)
τ= 3.07 (g)
τ= 3.08 (h)
τ= 3.15 (i)
τ= 3.19Figure4 shows the Neimark-Sacker bifurcation diagram with parameter τ ranging in [2.8, 3.2]. Under the above given parametric conditions, when τ<τ0, a locally asymptotically stable fixed point dominates in the phase space, whereas when τ>τ0 and the value of τ is at the nearby of τ0, the system dynamics converges to invariant circles. Via amplifying the bifurcation diagram in Figure 4(a), periodic windows with period-26 and period-7 orbits are found, as shown in Figures 4(b) and 4(c). Moreover, as the value of parameter τ continuously increases, the discrete predator-prey system will present chaotic behaviors. As demonstrated by the maximum Lyapunov exponent in Figure 5, the occurrence of chaos begins at about τ=3.0648. And in the chaotic region, periodic windows can repeatedly emerge (see Figures 5(b) and 5(c)).Figure6 shows a series of representative dynamical behaviors corresponding to the Neimark-Sacker bifurcation diagram. Figure 6(a) shows a locally asymptotically stable focus, Figures 6(b) and 6(d) display two different invariant circles, Figures 6(f), 6(h), and 6(i) all exhibit a cloud of random points, which suggest the emergence of chaos, and Figures 6(c), 6(e), and 6(g) demonstrate a few periodic windows occurring in-between the regions of invariant circles and chaos. The variation of these dynamical behaviors actually reveals the dynamic transition on the route to chaos induced by Neimark-Sacker bifurcation.


### 5.2. Bifurcations and Complex Dynamics of the Coupled Dynamical Network

To illustrate the nonchaos-chaos transition in the coupled dynamical network, here we present one example via applying the globally coupled network. Considering a globally coupled network withN nodes, the coupling matrix of which can be described as (60)Λ=-N+111⋯11-N+11⋯1⋮⋱⋱⋯1111⋯1111⋯-N+1.Through solving the solutions ofλE-Λ=0, the eigenvalues of the coupling matrix are determined as λ1=0 and λ2=λ3=⋯=λN=-N, respectively. According to (57), the maximum Lyapunov exponent of the coupled dynamical network is μN(λN)=h0+ln⁡1-cλN. When μN(λN)=h0+ln⁡1+cN>0, the coupled dynamical network exhibits chaotic behaviors. This condition can be further converted via (59) to determine the range of parameterτ where chaotic dynamics emerges in the coupled dynamical network (4a) and (4b).Figure7 exhibits the flip bifurcation diagram of the coupled dynamical network when N=5. Comparing with Figure 1, the bifurcation diagram here shows a few variations under the coupling effects of the discrete predator-prey system. First, the chaotic region extends. Following the calculation in previous section, the coupled network first experiences chaotic dynamics at about τ=0.2517, which comes earlier than that in the single discrete predator-prey system. Likewise, the periodic windows in Figure 1 also shrink or even disappear in the bifurcation diagrams of the coupled dynamical network. Second, the dynamics at the bifurcation critical point is greatly changed. As shown in Figures 7 and 8(b), invariant circles emerge when the value of parameter τ stays around τ∗. Third, new periodic window may appear in the chaotic region. As demonstrated in Figures 7(b) and 8(f), the new periodic window shows the existence of period-8 orbit, the occurrence of which results from that maximum Lyapunov exponent of the discrete predator-prey system is very small, and therefore the network dynamics is not chaotic and stays at the periodic orbit of the discrete predator-prey system.Figure 7
Flip bifurcation diagram of the coupled dynamical network with different coupling strengths. Globally coupled network is applied,N=5. Other parametric conditions are the same with Figure 1.
(a)
δ = 0.01 (b)
δ = 0.015 (c)
δ = 0.02Figure 8
Phase diagrams for different values of parameterτ corresponding to Figure 7(b), demonstrating the predator-prey dynamics at one node of the coupled network.
(a)
τ = 0.2 (b)
τ = 0.2052 (c)
τ = 0.252 (d)
τ = 0.256 (e)
τ = 0.262 (f)
τ = 0.263 (g)
τ = 0.265 (h)
τ = 0.27Figure8 displays the predator-prey dynamics in the phase diagrams at one node of the coupled network. As demonstrated in the figures, it can be found that the coupled dynamical network shows more complex dynamical behaviors than single discrete predator-prey system. At the bifurcation critical point and around, complex invariant circles occur and replace periodic behaviors (Figure 8(b)). Like the discrete predator-prey system, the coupled network also experiences period-doubling process, generating period-2, period-4, period-8 orbits. When the network dynamics goes into period-8 orbit zone, chaotic dynamics comes sequently, exhibiting various chaotic attractors (Figures 8(c)-8(d) and 8(f)–8(h)). In the chaotic region, the network dynamics once shortly falls back to periodic behaviors (Figure 8(e)).When the value ofN is given as 20, the flip bifurcation diagram of the coupled dynamical network remains similar with thatN = 5 (see Figure 9). However, with the increase of nodes, the phase diagrams of the coupled dynamical network may present more complex and interesting configurations (Figure 10). Around the bifurcation critical point τ∗, the phase diagrams of one node just show two simple invariant circles (Figures 10(a)-10(b)).Figure 9
Flip bifurcation diagram of the coupled dynamical network with different coupling strengths. Globally coupled network is applied,N=20. Other parametric conditions are the same with Figure 1.
(a)
δ = 0.001 (b)
δ = 0.002 (c)
δ = 0.003Figure 10
Phase diagrams corresponding to Figure9(c), demonstrating the predator-prey dynamics of the coupled network. (a)-(c) are plotted with the (u,v) states at one node.
(a)
τ = 0.206 (b)
τ = 0.207 (c)
τ = 0.25 (d)
τ = 0.254When the discrete predator-prey system undergoes Neimark-Sacker bifurcation, the bifurcation dynamics of the coupled network will show another situation. Comparing with Figure4, the change of the Neimark-Sacker bifurcation diagram under network coupling effect mainly reflects in the following aspects. First, the bifurcation critical point where dynamics transits from fixed point to attracted invariant circles occurs at a smaller value of parameter τ. In Figure 11, for δ=0.005, τ0=2.803; for δ=0.02, τ0=2.7465; for δ=0.05, τ0=2.6372. It can be found that the stronger the coupling effect, the smaller the bifurcation critical point. Second, theτvalue where the coupled network first goes into chaotic region comes earlier. Via the calculation in (59), the chaotic dynamics first takes place at about τ=2.8082, 2.7918, and 2.7847 for the three cases in Figure 11. Third, even though the dynamics of coupled network goes into the chaotic region, periodic windows can also emerge. Moreover, in Figure 11, there exists another dynamic window, where the network dynamics returns to invariant circles.Figure 11
Neimark-Sacker bifurcation diagram of the coupled dynamical network with different coupling strengths. Globally coupled network is applied,N=5. Other parametric conditions are the same with Figure 4.
(a)
δ = 0.005 (b)
δ = 0.02 (c)
δ = 0.02The Neimark-Sacker bifurcation of the coupled dynamical network can induce rich and complex dynamical behaviors as demonstrated in Figures12. When the value of parameter τ comes across the bifurcation critical point, the predator-prey dynamics at each node of the coupled network follows the behavior of invariant circle (Figure 12(a)). As the τ value grows, the network dynamics will suddenly become chaotic (Figure 12(b)). Following this chaotic region, the network dynamics returns to be nonchaotic, since the maximum Lyapunov exponent is smaller than zero. As shown in Figures 12(c)–12(g), different tori dominate in the phase diagrams. With further increasing τ value, dynamic windows of invariant circle and periodic orbit sequently take place (Figures 12(h)-12(i)). Finally, after the periodic window, the network dynamics again transits to chaos. Figures 12(j)–12(l) exhibit several chaotic attractors of the coupled dynamical network.Figure 12
Phase diagrams for different values of parameterτ corresponding to Figure 11(a), demonstrating the predator-prey dynamics at one node of the coupled network.
(a)
τ= 2.805 (b)
τ= 2.81 (c)
τ= 2.815 (d)
τ= 2.819 (e)
τ= 2.82 (f)
τ= 2.825 (g)
τ= 2.84 (h)
τ= 2.92 (i)
τ= 3.015 (j)
τ = 3.02 (k)
τ= 3.05 (l)
τ= 3.08Our numerical analysis reveals that the dynamic transition exhibits similar trend under other coupling strengths when the coupled network experiences Neimark-Sacker bifurcation. Figure13 demonstrates various network dynamical behaviors in the phase diagrams of one node. Like the nonchaos-chaos dynamic transition in Figure 12, the dynamics in Figure 13 sequently undergoes the alternation between invariant circles and tori, then the periodic window, and finally the chaos. However, the attractors of invariant circles, tori, and chaos exhibit different appearance with that in previous figures. This suggests the diversity and complexity of the predator-prey dynamics in the coupled network.Figure 13
Phase diagrams for different values of parameterτ corresponding to Figure 11(c), demonstrating the predator-prey dynamics at one node of the coupled network.
(a)
τ= 2.72 (b)
τ= 2.729 (c)
τ= 2.73 (d)
τ= 2.74 (e)
τ= 2.755 (f)
τ= 2.78 (g)
τ= 2.785 (h)
τ= 2.8Figure14 shows the Neimark-Sacker bifurcation diagram of the coupled dynamical network when the value of N is given as 20. Compared with Figure 4, the bifurcation critical point and first chaotic dynamics point here both move slightly toward the origin. Compared with Figure 11, the chaotic region in Figure 14 shows no windows for invariant circle and period orbit. Generally, the dynamic transition along the τ increasing direction in the bifurcation diagrams presents similar but simpler trend. Three regions in these diagrams can be divided. In the region of straight line, the network dynamics is attracted to fixed point. In the region of bullet shape, invariant circles prevail. In the third region, the coupled network exhibits tori and chaotic dynamics. At the intersection between the second and the third regions, new attractors for invariant circles or tori can often appear, such as the examples given in Figure 15. These attractors are also a reflection of rich dynamics of the coupled network.Figure 14
Neimark-Sacker bifurcation diagram of the coupled dynamical network with different coupling strengths. Globally coupled network is applied,N=20. Other parametric conditions are the same with Figure 4.
(a)
δ = 0.0015 (b)
δ = 0.002 (c)
δ = 0.005Figure 15
Phase diagrams corresponding to Figures14(b) and 14(c), demonstrating the predator-prey dynamics of the coupled network. (a)-(d) are plotted with (u,v) states at one node.
(a)
δ = 0.002,τ =2.82 (b)
δ = 0.002,τ =2.823 (c)
δ = 0.005,τ =2.783 (d)
δ = 0.005,τ =2.784


## 5.1. Bifurcations and Complex Dynamics of the Discrete Predator-Prey System

To exhibit the flip bifurcation of the discrete predator-prey system, we give a group of parameter values asa=0.1, b=0.001, c=10, d=0.05, and e=0.037 and choose parameter τ as bifurcation parameter. At the beginning, the occurrence of the flip bifurcation under such parametric conditions is demonstrated. The fixed point of the discrete predator-prey system is determined as (0.0125, 0.9897), and the critical point for flip bifurcation is calculated as τ∗=0.2052. When τ<τ∗, the fixed point (u3,v3) of the system is stable. Taking τ=0.2 for example, the two eigenvalues of the fixed point (u3,v3) are λ1=-0.9489, λ2=0.9466, the modulus of which are both less than one, suggesting (u3,v3) is a stable node. When τ=τ∗, we get the two eigenvalues of (u3,v3) as λ1=-1, λ2=0.9452. Then map (26) restricted to the center manifold is got as F:w~→-w~+4.0288w~2-9.7433w~τ~+57.8564w~2τ~-15.1681μ5w~3+o((w~+τ~)4). On the basis of the one-dimensional map F, the two discriminatory quantities η1 and η2 are got as η1=-9.7443, η2=1.0635. Hence, according to the flip bifurcation theorem, the discrete predator-prey system indeed undergoes flip bifurcation at τ=τ∗. Furthermore, because η2>0, we know that a period-2 orbit bifurcating from (u3,v3) is stable and the fixed point (u3,v3) becomes unstable. Taking τ=0.21, for example, the two eigenvalues of (u3,v3) are λ1=-1.0463, λ2=0.9439, respectively, therefore (u3,v3) is a saddle point. Simultaneously, the discrete predator-prey system is attracted to the two stable periodic points (0.0128, 0.8841) and (0.0131, 1.0923).Figure1(a) shows the flip bifurcation diagram of the discrete predator-prey system with parameter τ ranging in [0.2,0.3] under the above given parametric conditions. Explicitly, a period-doubling cascade is observed and the dynamic transition from periodic orbits to chaotic behaviors is found. Figures 1(b) and 1(c) are local amplifications of Figure 1(a), demonstrating the detail of periodic windows with period-6 and period-5 orbits. In the periodic windows, period-doubling cascade also occurs, resulting to period-10, period-20, period-40 orbits and period-12, period-24, period-48 orbits, etc. Moreover, other periodic windows can also be found in the chaotic region. Figure 2 shows the maximum Lyapunov exponents corresponding to the three bifurcation diagrams in Figure 1. It can be seen in Figure 2 that the discrete predator-prey system first goes into chaotic region at about τ=0.2642.Figure 1
(a) Flip bifurcation diagram of the discrete predator-prey system corresponding to variation ofτ ranging in [0.2,0.3]. a=0.1, b=0.001, c=10, d=0.05, and e=0.037. (b) and (c) Local amplifications of (a) in [0.27,0.2712] and [0.2812,0.2822] showing the detail of periodic windows.
(a) (b) (c)Figure 2
Diagrams of maximum Lyapunov exponent corresponding to the three graphs of Figure1.
(a) (b) (c)Figure3 shows different attractors in the phase portraits corresponding to the flip bifurcation diagram in Figure 1(a). Figures 3(a)–3(d) display the period-1, period-2, period-4, period-8 orbits. The doubling of periodic points in these orbits coincides with the period-doubling process in Figure 1(a). Figures 3(e), 3(g), and 3(i) show a cloud of random points, indicating chaotic behaviors. The maximum Lyapunov exponents of these dynamical behaviors are 0.0756, 0.1334, and 0.3265, respectively, confirming that these behaviors are indeed chaotic. Figures 3(f) and 3(h) display the period-6 and period-5 orbits; in fact, they correspond to the periodic windows as shown in Figures 1(b) and 1(c). The dynamic transition from period to chaos demonstrated in Figure 3 suggests a disordering process occurring in the discrete predator-prey system.Figure 3
Phase diagrams for different values of parameterτ corresponding to Figure 1(a).
(a)
τ = 0.2 (b)
τ = 0.24 (c)
τ = 0.26 (d)
τ = 0.263 (e)
τ = 0.267 (f)
τ = 0.2703 (g)
τ =0.2715 (h)
τ = 0.2814 (i)
τ = 0.295Figures4–6 demonstrate the dynamics of Neimark-Sacker bifurcation of the discrete predator-prey system. Before describing these figures, an example is provided to illustrate the occurrence of Neimark-Sacker bifurcation. Given parameter values as a=0.1, b=1.5, c=0.5, d=3.5, and e=0.05, the fixed point (u3,v3) is obtained as (0.8313, 0.2518), and the critical point for Neimark-Sacker bifurcation is determined as τ0=2.8234. When τ<τ0, taking τ=2.6, for instance, we get the two eigenvalues of (u3,v3) are λ1,2=-0.7138±0.4682i~ and λ1,2=0.8537; therefore (u3,v3) is a stable focus. When we take τ=τ0, the two eigenvalues are λ1=-0.8611±0.5085i~ and λ1,2=1. Besides, through calculation, we have γ=0.6592 and a=-16.8134. It is verified that the discrete predator-prey system undergoes Neimark-Sacker bifurcation at τ=τ0. Moreover, on account of γ>0 and a<0, we have an attracting invariant closed curve bifurcating from the fixed point for τ>τ0.Figure 4
(a) Neimark-Sacker bifurcation diagram of the discrete predator-prey system corresponding to variation ofτ ranging in [2.8,3.2]. a=0.1, b=1.5, c=0.5, d=3.5, and e=0.05. (b) and (c) Local amplifications of (a) in [3.045,3.07] and [3.065,3.11] showing the detail of periodic windows.
(a) (b) (c)Figure 5
Diagrams of maximum Lyapunov exponent corresponding to the three graphs of Figure4.
(a) (b) (c)Figure 6
Phase diagrams for different values of parameterτ corresponding to Figure 4(a).
(a)
τ= 2.815 (b)
τ= 2.9 (c)
τ= 3.015 (d)
τ= 3.037 (e)
τ= 3.06 (f)
τ= 3.07 (g)
τ= 3.08 (h)
τ= 3.15 (i)
τ= 3.19Figure4 shows the Neimark-Sacker bifurcation diagram with parameter τ ranging in [2.8, 3.2]. Under the above given parametric conditions, when τ<τ0, a locally asymptotically stable fixed point dominates in the phase space, whereas when τ>τ0 and the value of τ is at the nearby of τ0, the system dynamics converges to invariant circles. Via amplifying the bifurcation diagram in Figure 4(a), periodic windows with period-26 and period-7 orbits are found, as shown in Figures 4(b) and 4(c). Moreover, as the value of parameter τ continuously increases, the discrete predator-prey system will present chaotic behaviors. As demonstrated by the maximum Lyapunov exponent in Figure 5, the occurrence of chaos begins at about τ=3.0648. And in the chaotic region, periodic windows can repeatedly emerge (see Figures 5(b) and 5(c)).Figure6 shows a series of representative dynamical behaviors corresponding to the Neimark-Sacker bifurcation diagram. Figure 6(a) shows a locally asymptotically stable focus, Figures 6(b) and 6(d) display two different invariant circles, Figures 6(f), 6(h), and 6(i) all exhibit a cloud of random points, which suggest the emergence of chaos, and Figures 6(c), 6(e), and 6(g) demonstrate a few periodic windows occurring in-between the regions of invariant circles and chaos. The variation of these dynamical behaviors actually reveals the dynamic transition on the route to chaos induced by Neimark-Sacker bifurcation.


## 5.2. Bifurcations and Complex Dynamics of the Coupled Dynamical Network

To illustrate the nonchaos-chaos transition in the coupled dynamical network, here we present one example via applying the globally coupled network. Considering a globally coupled network withN nodes, the coupling matrix of which can be described as (60)Λ=-N+111⋯11-N+11⋯1⋮⋱⋱⋯1111⋯1111⋯-N+1.Through solving the solutions ofλE-Λ=0, the eigenvalues of the coupling matrix are determined as λ1=0 and λ2=λ3=⋯=λN=-N, respectively. According to (57), the maximum Lyapunov exponent of the coupled dynamical network is μN(λN)=h0+ln⁡1-cλN. When μN(λN)=h0+ln⁡1+cN>0, the coupled dynamical network exhibits chaotic behaviors. This condition can be further converted via (59) to determine the range of parameterτ where chaotic dynamics emerges in the coupled dynamical network (4a) and (4b).Figure7 exhibits the flip bifurcation diagram of the coupled dynamical network when N=5. Comparing with Figure 1, the bifurcation diagram here shows a few variations under the coupling effects of the discrete predator-prey system. First, the chaotic region extends. Following the calculation in previous section, the coupled network first experiences chaotic dynamics at about τ=0.2517, which comes earlier than that in the single discrete predator-prey system. Likewise, the periodic windows in Figure 1 also shrink or even disappear in the bifurcation diagrams of the coupled dynamical network. Second, the dynamics at the bifurcation critical point is greatly changed. As shown in Figures 7 and 8(b), invariant circles emerge when the value of parameter τ stays around τ∗. Third, new periodic window may appear in the chaotic region. As demonstrated in Figures 7(b) and 8(f), the new periodic window shows the existence of period-8 orbit, the occurrence of which results from that maximum Lyapunov exponent of the discrete predator-prey system is very small, and therefore the network dynamics is not chaotic and stays at the periodic orbit of the discrete predator-prey system.Figure 7
Flip bifurcation diagram of the coupled dynamical network with different coupling strengths. Globally coupled network is applied,N=5. Other parametric conditions are the same with Figure 1.
(a)
δ = 0.01 (b)
δ = 0.015 (c)
δ = 0.02Figure 8
Phase diagrams for different values of parameterτ corresponding to Figure 7(b), demonstrating the predator-prey dynamics at one node of the coupled network.
(a)
τ = 0.2 (b)
τ = 0.2052 (c)
τ = 0.252 (d)
τ = 0.256 (e)
τ = 0.262 (f)
τ = 0.263 (g)
τ = 0.265 (h)
τ = 0.27Figure8 displays the predator-prey dynamics in the phase diagrams at one node of the coupled network. As demonstrated in the figures, it can be found that the coupled dynamical network shows more complex dynamical behaviors than single discrete predator-prey system. At the bifurcation critical point and around, complex invariant circles occur and replace periodic behaviors (Figure 8(b)). Like the discrete predator-prey system, the coupled network also experiences period-doubling process, generating period-2, period-4, period-8 orbits. When the network dynamics goes into period-8 orbit zone, chaotic dynamics comes sequently, exhibiting various chaotic attractors (Figures 8(c)-8(d) and 8(f)–8(h)). In the chaotic region, the network dynamics once shortly falls back to periodic behaviors (Figure 8(e)).When the value ofN is given as 20, the flip bifurcation diagram of the coupled dynamical network remains similar with thatN = 5 (see Figure 9). However, with the increase of nodes, the phase diagrams of the coupled dynamical network may present more complex and interesting configurations (Figure 10). Around the bifurcation critical point τ∗, the phase diagrams of one node just show two simple invariant circles (Figures 10(a)-10(b)).Figure 9
Flip bifurcation diagram of the coupled dynamical network with different coupling strengths. Globally coupled network is applied,N=20. Other parametric conditions are the same with Figure 1.
(a)
δ = 0.001 (b)
δ = 0.002 (c)
δ = 0.003Figure 10
Phase diagrams corresponding to Figure9(c), demonstrating the predator-prey dynamics of the coupled network. (a)-(c) are plotted with the (u,v) states at one node.
(a)
τ = 0.206 (b)
τ = 0.207 (c)
τ = 0.25 (d)
τ = 0.254When the discrete predator-prey system undergoes Neimark-Sacker bifurcation, the bifurcation dynamics of the coupled network will show another situation. Comparing with Figure4, the change of the Neimark-Sacker bifurcation diagram under network coupling effect mainly reflects in the following aspects. First, the bifurcation critical point where dynamics transits from fixed point to attracted invariant circles occurs at a smaller value of parameter τ. In Figure 11, for δ=0.005, τ0=2.803; for δ=0.02, τ0=2.7465; for δ=0.05, τ0=2.6372. It can be found that the stronger the coupling effect, the smaller the bifurcation critical point. Second, theτvalue where the coupled network first goes into chaotic region comes earlier. Via the calculation in (59), the chaotic dynamics first takes place at about τ=2.8082, 2.7918, and 2.7847 for the three cases in Figure 11. Third, even though the dynamics of coupled network goes into the chaotic region, periodic windows can also emerge. Moreover, in Figure 11, there exists another dynamic window, where the network dynamics returns to invariant circles.Figure 11
Neimark-Sacker bifurcation diagram of the coupled dynamical network with different coupling strengths. Globally coupled network is applied,N=5. Other parametric conditions are the same with Figure 4.
(a)
δ = 0.005 (b)
δ = 0.02 (c)
δ = 0.02The Neimark-Sacker bifurcation of the coupled dynamical network can induce rich and complex dynamical behaviors as demonstrated in Figures12. When the value of parameter τ comes across the bifurcation critical point, the predator-prey dynamics at each node of the coupled network follows the behavior of invariant circle (Figure 12(a)). As the τ value grows, the network dynamics will suddenly become chaotic (Figure 12(b)). Following this chaotic region, the network dynamics returns to be nonchaotic, since the maximum Lyapunov exponent is smaller than zero. As shown in Figures 12(c)–12(g), different tori dominate in the phase diagrams. With further increasing τ value, dynamic windows of invariant circle and periodic orbit sequently take place (Figures 12(h)-12(i)). Finally, after the periodic window, the network dynamics again transits to chaos. Figures 12(j)–12(l) exhibit several chaotic attractors of the coupled dynamical network.Figure 12
Phase diagrams for different values of parameterτ corresponding to Figure 11(a), demonstrating the predator-prey dynamics at one node of the coupled network.
(a)
τ= 2.805 (b)
τ= 2.81 (c)
τ= 2.815 (d)
τ= 2.819 (e)
τ= 2.82 (f)
τ= 2.825 (g)
τ= 2.84 (h)
τ= 2.92 (i)
τ= 3.015 (j)
τ = 3.02 (k)
τ= 3.05 (l)
τ= 3.08Our numerical analysis reveals that the dynamic transition exhibits similar trend under other coupling strengths when the coupled network experiences Neimark-Sacker bifurcation. Figure13 demonstrates various network dynamical behaviors in the phase diagrams of one node. Like the nonchaos-chaos dynamic transition in Figure 12, the dynamics in Figure 13 sequently undergoes the alternation between invariant circles and tori, then the periodic window, and finally the chaos. However, the attractors of invariant circles, tori, and chaos exhibit different appearance with that in previous figures. This suggests the diversity and complexity of the predator-prey dynamics in the coupled network.Figure 13
Phase diagrams for different values of parameterτ corresponding to Figure 11(c), demonstrating the predator-prey dynamics at one node of the coupled network.
(a)
τ= 2.72 (b)
τ= 2.729 (c)
τ= 2.73 (d)
τ= 2.74 (e)
τ= 2.755 (f)
τ= 2.78 (g)
τ= 2.785 (h)
τ= 2.8Figure14 shows the Neimark-Sacker bifurcation diagram of the coupled dynamical network when the value of N is given as 20. Compared with Figure 4, the bifurcation critical point and first chaotic dynamics point here both move slightly toward the origin. Compared with Figure 11, the chaotic region in Figure 14 shows no windows for invariant circle and period orbit. Generally, the dynamic transition along the τ increasing direction in the bifurcation diagrams presents similar but simpler trend. Three regions in these diagrams can be divided. In the region of straight line, the network dynamics is attracted to fixed point. In the region of bullet shape, invariant circles prevail. In the third region, the coupled network exhibits tori and chaotic dynamics. At the intersection between the second and the third regions, new attractors for invariant circles or tori can often appear, such as the examples given in Figure 15. These attractors are also a reflection of rich dynamics of the coupled network.Figure 14
Neimark-Sacker bifurcation diagram of the coupled dynamical network with different coupling strengths. Globally coupled network is applied,N=20. Other parametric conditions are the same with Figure 4.
(a)
δ = 0.0015 (b)
δ = 0.002 (c)
δ = 0.005Figure 15
Phase diagrams corresponding to Figures14(b) and 14(c), demonstrating the predator-prey dynamics of the coupled network. (a)-(d) are plotted with (u,v) states at one node.
(a)
δ = 0.002,τ =2.82 (b)
δ = 0.002,τ =2.823 (c)
δ = 0.005,τ =2.783 (d)
δ = 0.005,τ =2.784


## 6. Conclusions

Generally, the dynamics of coupled network is far richer and more complex than single dynamic system considered without system coupling. Previous research has been made, revealing that even though the dynamic system at each node holds parametric conditions which generate nonchaotic dynamic behaviors, the coupled dynamical network can still exhibit chaotic dynamics [4, 14]. On the basis of former research works, this research investigates the complex dynamic behaviors and nonlinear characteristics in a coupled network of discrete predator-prey system. New improvements are obtained on the bifurcations and nonchaos-chaos transition of the coupled dynamical network. Via comparing the dynamics between single discrete predator-prey system and the coupled dynamical network, how the coupled network influences the predator-prey dynamics is determined.The complex behaviors of the discrete predator-prey system are described firstly. The system may have four fixed points, among which the stable unique positive fixed point reflects the predator and prey can stably coexist. Moreover, the discrete predator-prey system can undergo flip bifurcation and Neimark-Sacker bifurcation at this fixed point. The two bifurcations both trigger a route to chaos, which is also a disordering process from period to chaos. On the routes to chaos, the system displays complex behaviors and various attractors, including period-doubling cascades in orbits of period-2, 4, 8, periodic windows with period-5, 6, 7, 10, 12, 20, 24, 26, 40, 48 orbits, invariant circles, and chaotic sets, which imply that the predator and prey can coexist in many possible complicated ways.In comparing with the discrete predator-prey system, the coupled dynamical network can exhibit much more interesting behaviors. To illustrate the complex dynamics and nonchaos-chaos transition in the coupled network, globally coupled network is taken for investigation. The flip bifurcation and Neimark-Sacker bifurcation also emerge and induce dynamic transitions from fixed point to chaotic dynamics. Very complicated behaviors are found with the increase of varying parameter. With the findings of this research, the network dynamics can be concluded as follows:(1)
Compared with the flip bifurcation of single discrete predator-prey system, the flip bifurcation of the coupled network also triggers period-doubling cascade with period-2, 4, 8 orbits. However, the network dynamics at the bifurcation critical point turns to be invariant circles. Furthermore, the coupled network goes into chaotic dynamics earlier and the chaotic region is extended. New periodic window (such as period-8 orbit) may appear in the chaotic region.(2)
Compared with the Neimark-Sacker bifurcation of single discrete predator-prey system, the bifurcation critical point where dynamic behaviors transit from fixed point to attracted invariant circle and the point where the dynamics firstly becomes chaotic both comes earlier in the coupled dynamical network.(3)
In the coupled network, the dynamic transition on the routes to chaos exhibits complex process. When flip bifurcation occurs, the network dynamics may sequently experience invariant circles, period-doubling cascade, chaotic sets, and periodic windows and finally stay at the chaotic region. When Neimark-Sacker bifurcation takes place, the coupled network may sequently undergo invariant circles, chaotic sets, tori, invariant circles, periodic orbits, and eventually chaotic dynamics.



---

*Source: 2583730-2019-06-13.xml*