Abstract
Predator-prey models describe biological phenomena of pursuit-evasion interaction. And this interaction exists widely in the world for the necessary energy supplement of species. In this paper, we have investigated a ratio-dependent spatially extended food chain model. Based on the bifurcation analysis (Hopf and Turing), we give the spatial pattern formation via numerical simulation, that is, the evolution process of the system near the coexistence equilibrium point , and find that the model dynamics exhibits complex pattern replication. For fixed parameters, on increasing the control parameter , the sequence “holes holes-stripe mixtures stripes spots-stripe mixtures spots” pattern is observed. And in the case of pure Hopf instability, the model exhibits chaotic wave pattern replication. Furthermore, we consider the pattern formation in the case of which the top predator is extinct, that is, the evolution process of the system near the equilibrium point , and find that the model dynamics exhibits stripes-spots pattern replication. Our results show that reaction-diffusion model is an appropriate tool for investigating fundamental mechanism of complex spatiotemporal dynamics. It will be useful for studying the dynamic complexity of ecosystems.
1. Introduction
Predator-prey models are studied in detail in the focus on equilibria, stability, asymptotic behavior, persistence, bifurcation, chaos, and so on [1–8]. In the past 40 years, with the idea of Turing [9], spatial extended models, in which not only the species evolve through time but also distribute in space, and pattern formation are one of the hot spots [6, 8, 10–20].
Food web models describe the same phenomena as predator-prey models, but the former description is more actual than the latter since our real world is so complex. Until recently, food webs models are widely studied as predator-prey models [12–14, 17, 21–29]. But as far as we know, spatially extended models seem rare and not regarded. In fact, we live in a spatial world, and the spatial component of ecological interactions has been identified as an important factor in how ecological communities are shaped. Understanding the role of space is challenging both theoretically and empirically [30]. And the issue of spatial and spatiotemporal pattern formation in biological communities is probably one of the most exciting problems in modern biology and ecology [31, 32]. And the food web models with spatial distribution will do better job than the classical models.
In general, a classical food chain model with the nondimensional form can be written as follows: where stands for prey density, is the top predator density, and —the density of the intermediate predator—describes the predator of and the prey of ; is the per capita rate of increase of the prey in the absence of predation. And all coefficients are positive constants, and are the maximum ingestion rates of intermediate predator and top predator, is the conversion factor of prey to intermediate predator, is the conversion factor of intermediate predator to top predator, is the food-independent death rate of the intermediate predator, and is the food-independent death rate of the top predator. is the functional response. The functional response is the prey consumption rate by an average single predator. It can be influenced by the prey consumption rate and the predator density. is the amount of prey consumed per predator per unit time, is the predator production per capita with predation.
In this paper, we focus on the following ratio-dependent food chain model [28]: The necessary condition of the persistence of and is and , respectively.
When all the species distribute randomly in the space, model (2) can be rewritten with a supplement: where , , and are the diffusion coefficients of the three species, respectively, is the usual Laplacian operator in two-dimensional space, and other parameters have the same definitions as those above.
Model (3) is to be analyzed under the nonzero initial condition and Neumann, or zero flux, boundary conditions: In the above, is the outward unit normal vector of the boundary which we will assume is smooth. The main reason for choosing such boundary conditions is that we are interested in the self-organization of pattern; zero-flux conditions imply no external input [17].
This paper is organized as follows. In the next section, we give a local stability analysis of model (3). Then, we present the pattern formation of model (3) via numerical simulations, which is followed by Section 2. Finally, we give some discussions in Section 4.
2. Linear Stability Analysis
There are two equilibria (steady states) in model (2), which correspond to spatially homogeneous equilibria of model (3): corresponding to top-predator extinction when . Consider corresponding to coexistence of prey and predators when or or where
And in the presence of diffusion, set , , , and the standard linear analysis predicts exponentially growing solutions of model (3) in the form where , , and are the wave-number and frequency, respectively.
And the eigenvalue equation then reads where the diffusion matrix , and the Jacobian matrix
Then we can obtain the eigenvalues as functions of the wave number as the roots of where
And one type of bifurcation will break one type of symmetry of a system; that is, in the bifurcation point, two equilibrium states intersect and exchange their stability. Biologically speaking, this bifurcation corresponds to a smooth transition between equilibrium states [33]. The reaction-diffusion systems have led to the characterization of two basic types of symmetry-breaking bifurcations—Hopf and Turing bifurcation, which are responsible for the emergence of spatiotemporal patterns.
The onset of Hopf instability corresponds to the case when a pair of imaginary eigenvalues cross the real axis from the negative to the positive side. And this situation occurs only when the diffusion vanishes. Mathematically speaking, the Hopf bifurcation occurs when , at the wavenumber . For unstable steady states to heterogeneous perturbations leading to Turing patterns, the real part of the eigenvalue, , has to be greater than zero. Mathematically speaking, the Turing bifurcation occurs when , at the wavenumber .
Here, we take as the bifurcation parameter; linear stability analysis yields the bifurcation diagram with , , , , , , , and is a variational parameter (c.f., Figure 1). In Figure 1, the spotted curve is critical state in which above the spotted curve, the three species cannot both be positive; under the spotted curve, they are both positive. The bifurcation diagram shows the two bifurcation curves separate the coexistence space into four domains. In domain I, located below all two bifurcation lines, the uniform steady state is the only stable solution of the model. Domain II is the region of pure Turing instability. Domain III is the region of pure Hopf instability. When the parameters correspond to domain IV, which is located above all two bifurcation lines, both Hopf instability and Turing instability occur.
In Figure 1, the stationary state in the parameter domains II and IV (sometimes called the “Turing space”) is unstable only to a nonuniform perturbation. As expected, this domain exists only when the inhibitor species (for predator-prey system, predator ) diffuses faster than the activator species (for predator-prey system, prey ) and the area of this Turing space increases with .
3. Pattern Formation
In this section, we perform extensive numerical simulations of the spatially extended model (3) in two-dimensional spaces, and the qualitative results are shown here. The parameters are , , , , , , and . Model (3) is integrated initially in two-dimensional space from the homogeneous steady state; that is, we start with the unstable uniform solution with small random perturbation superimposed; in each, the initial condition is always a small amplitude random perturbation , using an explicit Euler method for the time integration with a time stepsize of . We use the standard five-point approximation for the Laplacian operator with the Zero-flux boundary conditions and the system size is space units with space stepsize (lattice constant) , discretized through and . The form of the Laplacian operator is taken as follows:
The concentrations at the moment at the mesh position are given by
When the evolution processes reached steady state, we took a snapshot with white corresponding to the high value of prey while black corresponding to the low one.
In the numerical simulations, different types of dynamics are observed and it is found that the distributions of predator and prey are always of the same type. Consequently, we can restrict our analysis of pattern formation to one distribution. In this section, we show the distribution of prey , for instance.
From the bifurcation diagram in the above section (cf., Figure 1), the results of numerical simulations show that the type of the system dynamics is determined by the values of and . And for different sets of parameters, the features of the spatial patterns become essentially different if exceeds the bifurcation curves which depend on .
First, we consider the pattern formation for the parameters located in domain II (c.f., Figure 1); the region of pure Turing instability occurs while Hopf stability occurs. As an example, we show the time evolution of three typical patterns when . With the parameters set, one can conclude that the critical value of Hopf bifurcation is and Turing bifurcation value is . So, the values of that we adopt are between and .
As an example, in Figure 2, we show the time evolution of holes pattern of prey at 0, 20000, 60000, and 200000 iterations for . In this case, one can see that for model (1), the pattern takes a long time to settle down, starting with a homogeneous state (c.f., Figure 5(a)), and the random initial distribution leads to the formation of regular holes (c.f., Figure 2(d)). This pattern (c.f., Figure 2(d)) consists of black (minimum density of ) hexagons on a white (maximum density of ) background, that is, isolated zones with low population densities. Baurmann et al. [33] called this type pattern “cold spots” and von Hardenberg et al. [34] called it “holes.” In this paper, we adopt the name “holes.”
(a)
(b)
(c)
(d)
When increasing to , a few of stripes emerge, and the remainder of the holes pattern remains time independent (Figure 3(a)). And while increasing to , model dynamics exhibits a transition from stripe-hole growth to stripes replication; that is, holes decay and the stripes pattern emerges (Figure 3(b)).
(a)
(b)
Next, we consider the pattern formation in domain IV in Figure 1; both Hopf and Turing instability occur in this domain. We adopt and —the maximized value of the coexistence of prey and their predators. The model dynamics exhibits two typical pattern formations.
In Figure 4, with the increasing of to 2.1, a few of white hexagons (i.e., spots, associate with high population densities) fill in the stripes; that is, the stripes-spots pattern emerges (c.f., Figure 4(a)). And while increasing to , model dynamics exhibits a transition from stripe-spots growth to spots replication; that is, stripes decay and the spots pattern emerges(c.f., Figure 4(b)).
(a)
(b)
(a)
(b)
(c)
(d)
From Figures 2–4, one can see that, with fixed parameters, on increasing the control parameter , the sequence “holes holes-stripes mixtures stripes spots-stripes mixtures spots” pattern is observed.
In addition, we consider the pattern formation when locates in domain III in Figure 1, pure Hopf instability occurs. Figure 5 shows the evolution of the chaotic wave pattern of prey at 0, 50000, 100000, and 200000 iterations with . With these fixed parameters, the critical value of Hopf bifurcation is and the Turing bifurcation values equal . In order to make it clearer, in Figure 6, we show oscillate time series plots of (c.f., Figures 6(a), 6(b), and 6(c)), respectively. And phase portrait (c.f., Figure 6(d)) shows that there exhibits the “local” phase plane of the system obtained in a fixed point inside the region invaded by the irregular spatiotemporal oscillations.
(a)
(b)
(c)
(d)
Furthermore, we restrict our attention to the case when the top predator vanishes. Extinction of the top predator is studied by Chiu and coworkers; they gave a criterion for the extinction of top predator [35]. Here, we will illustrate the pattern formation about this case.
According to food chain model, describes extinction of the top predator. With the same method and the same parameters in Section 2, the bifurcation diagram is shown in Figure 7. In Figure 7, the spotted curve is critical state in which the domain above the spotted curve is noncoexistence space; the domain under the spotted curve is coexistence space. Only Turing curve intersects with the spotted curve, and it separates the coexistence space into two domains. When locates in domain I, under the Turing curve, the steady state is only stable solution of model (3); when locates in domain II in Figure 7, pure Turing instability occurs. That is to say, domain II is the “Turing space” only.
Figure 8 shows the evolution of the spatial pattern of prey at 0, 10000, 100000, and 300000 iterations with and ; that is, point locates in domain II in Figure 7. The random initial distribution around the steady state leads to the formation of stripes-spots pattern (c.f., Figure 8(d)).
(a)
(b)
(c)
(d)
4. Conclusions and Remarks
In summary, we have investigated a ratio-dependent spatially extended food chain model. Based on the bifurcation analysis (Hopf and Turing), we give the spatial pattern formation via numerical simulation. For the coexistence equilibrium point , we find that the model dynamics exhibits complex pattern replication, such as holes, holes-stripes, stripes, spots-stripes, spots, and chaotic wave pattern. And for the extinction of the top predator equilibrium point , we find that the model dynamics exhibits stripes-spots pattern replication.
In fact, in our world, every day, hundreds of species are extinct, and the extinction of a species is a fearful thing. And the top predator is extinct because there is a balance between the prey and the intermediate predator . In the case we considered, the density of the intermediate predator is not small, but very big. The intermediate predator is strong enough to fight back the top predator .
On the other hand, in the analysis of bifurcations (i.e., Hopf and Turing), we find that huge-sized computations are required, so we have to obtain more help via computers. In fact, computer-aided analysis is useful for nonlinear analysis. And computers have played an important role throughout the history of ecology. Today, numerical simulations also play an important role in spatial ecology. There are some international mathematical softwares, such as Matlab, Maple, and Mathematica, all of which have powerful function library and can provide scientific calculation and programming with friendly platform. We have finished all our symbolic computations in Maple and obtained our pattern snapshots (i.e., numerical simulations) in Matlab as Maple is more superior in symbolic computations while Matlab is more superior in numerical computations.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This research was supported by NSFC no. 11071273.