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.
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. Section 2 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: in which and represent population densities of prey and predator at the time of t respectively; R1, R2, A, B, m, D, E, and K are positive constants; R1 and are the intrinsic growth rate and the carrying capacity of the prey; the constants , , and are the saturating Crowley-Martin type functional response parameters, where measures the magnitude of interference among prey, expresses the interference among the predator, and represents the maximum consumption rate; R2 describes the growth rate of predator, is the maximum value which per capita reduction rate of can attain, and 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, . Then (1a) and (1b) can be transformed intowhere , , , , and .
Applying the forward Euler scheme to (2a) and (2b), the discrete-time Leslie-Gower predator-prey system can be expressed by the following equations:in which represents the sequence number of iterations and describes the time step size. It should be noticed that .
Consider a coupled dynamical network of 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]where and denote the sequence number of the nodes in the coupled dynamical network, represents the coupling strength of the network, and denotes the coupling configuration of the network. If there is a connection between node and node , then ; otherwise, . Simultaneously, let , , where is the degree of node and can be defined to be the number of its outreaching connections, 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.,
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:
Direct calculation obtains four nonnegative fixed points, i.e.,where is a positive real root of the following cubic equation:
If the inequalities are satisfied then the map (6) has an unique positive equilibrium ,
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
The four fixed points , , , and are substituted into the Jacobian matrix (12), respectively, and the eigenvalues and are calculated. If and are found, the corresponding fixed point is stable; if or , an unstable fixed point emerges. The stability of the four fixed points can be determined as follows.(1)The two eigenvalues of are and . Obviously, we have and , then we know is unstable.(2)The two eigenvalues of are and . Hence, the conditions for that is stable can be determined as(3)The two eigenvalues of are and . Since , the fixed point is also unstable.(4)For , the corresponding Jacobian matrix is where Calculating the two eigenvalues of as where
By calculating and , the stable conditions for are determined as
3.2. Flip Bifurcation Analysis
Flip bifurcation is performed around the fixed point 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 is equal to −1. In such case, we have , and direct calculation yieldswhere
When condition (19) establishes, the other eigenvalue of is obtained as . For the occurrence of flip bifurcation, we need , namely,
Let and we can translate the fixed point to the origin. Simultaneously, via doing Taylor expansion near the fixed point, the map (6) can be rewritten as where , , , and are described in (15) andIn (24), all the τ values satisfy .
On the basis of map (23), we make a reversible transformation as follows: Then the map (23) is transformed to a norm form in which
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 existing at the fixed point (0,0,0), which can be approximated byIn (28), is assumed to beAccording to , we can get
Using the balance of (30) to determine the values of , , , and , we get the following results:
Consider the dynamics of map (26) restricted in the central manifold . Hence, a one-dimensional map is obtained, i.e.,
The coefficients in map (32) are expressed as follows:
If the map (32) experiences flip bifurcation, two determinant values, and , should be nonzero, i.e., It can be obtained by direct calculation,
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 when the parameter varies in a small neighborhood of . Moreover, if , the period-2 orbit bifurcates from which is stable; if , the period-2 orbit bifurcates from 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 and , which lead to
When parametric conditions satisfy (36a) and (36b), we translate the fixed point to the origin by the translation Applying the Taylor expansion to map (6), we obtain a new mapwhere the coefficients in map (38) are described in (24), but with . 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 aswhere , and explicitly, .
At the same time, the occurrence of Neimark-Sacker bifurcation needs that the derivation of to parameter is not zero at . With a calculation, we haveOn the other hand, to avoid that is a real number or a pure imaginary number, we also need that Condition (41) is equivalent to . Since , we know . Hence, condition (41) becomes , i.e.,
Then we study the norm form of map (38). With the invertible transformationthe map (38) can be transformed into the following new map:where
Calculate the two- and three-order derivatives of and at and . The following results are obtained: