Abstract

We focus on a spatially extended Holling-type IV predator-prey model that contains some important factors, such as noise (random fluctuations), external periodic forcing, and diffusion processes. By a brief stability and bifurcation analysis, we arrive at the Hopf and Turing bifurcation surface and derive the symbolic conditions for Hopf and Turing bifurcation on the spatial domain. Based on the stability and bifurcation analysis, we obtain spiral pattern formation via numerical simulation. Additionally, we study the model with a color noise and external periodic forcing. From the numerical results, we know that noise or external periodic forcing can induce instability and enhance the oscillation of the species density, and the cooperation between noise and external periodic forces inherent to the deterministic dynamics of periodically driven models gives rise to the appearance of a rich transport phenomenology. Our results show that modeling by reaction-diffusion equations is an appropriate tool for investigating fundamental mechanisms of complex spatiotemporal dynamics.

1. Introduction

Predation, a complex natural phenomenon, exists widely in the world, for example, the sea, the plain, the forest, the desert, and so on [1]. To model this phenomenon, the predator-prey model has been suggested for a long time since the pioneering works of Kendall [2]. Predator-prey model is a kind of “pursuit and evasion” system in which the prey trie to evade the predator and the predator tries to catch the prey if they interact [3]. Pursuit means the predator tries to shorten the spatial distance between the predator and the prey; evasion means the prey tries to widen this spatial distance. In fact, predator-prey model is a mathematical method to approximate some part of our real world. And the dynamic behavior of predator-prey model has long been and will continue to be one of the dominant themes in both ecology and mathematical ecology due to its universal existence and importance [4, 5].

In general, a classical predator-prey model can be written as the form [6, 7] where and stand for prey and predator quantity, respectively, is the per capita rate of increase of the prey in absence of predation, is the food-independent death rate of predator, is the functional response, the prey consumption rate by an average single predator, which obviously increases with the prey consumption rate and can be influenced by the predator density, which refers to the change in the density of prey attached per unit time per predator as the prey density changes, is the amount of prey consumed per predator per unit time, and is the predator production per capita with predation.

In population dynamics, a functional response describes the relationship between the predator and their prey, and the predator-prey model is always named after the corresponding functional response for its key position [69]. In the history of population ecology, both ecologists and mathematicians have a great interest in the Holling-type predator-prey models [3, 8, 1021], including Holling-types I–III, originally due to Holling [22, 23], and Holling-type IV, suggested by Andrews [24]. The Holling-type functional responses are the so-called “prey-dependent” type [8], for in (1) is a function only related to prey . The classical expression of Holling-type II functional response is , and is called Holling-type III. The Holling-type IV functional response is written as follows: Function (2) is called Monod-Haldane-type functional response too [25]. In addition, when , a simplified form is proposed by Sokol [26], and some scholars also called it as Holling-type IV [9, 25]. In this paper, we focus on the Holling-IV functional response taking the form (2), and the corresponding predator-prey model takes the form where stands for maximum per capita growth rate of the prey, is the capture rate, is the conversion rate of prey captured by predator, is the food-independent death rate of predator, is the carrying capacity, and is the so-called half-saturation constant; such that the denominator of above system does not vanish for nonnegative .

On the other hand, the real world we live in is a spatial world, and spatial patterns are ubiquitous in nature, which modify the temporal dynamics and stability properties of population density at a range of spatial scales, whose effects must be incorporated in temporal ecological models that do not represent space explicitly [27]. And the spatial component of ecological interactions has been identified as an important factor in how ecological communities are shaped. Empirical evidence suggests that the spatial scale and structure of the environment can influence population interactions and the composition of communities [1].

The reaction-diffusion model is a typical spatially extended model. It considers not only time but also space and consists of several species which react with each other and diffuse within the spatial domain. It involves a pair of partial differential equations and represents the time course of reacting and diffusing process. In spatially extended predator-prey model, the interaction between the predator and the prey is the reaction item, and the diffusion item comes to being for the predator’s “pursuit” and the prey’s “evasion.” Diffusion is a spatial process, and the whole model describes the evolution of the predator and the prey going with time.

Decades after Turing [28] demonstrated that spatial patterns could arise from the interaction of reactions or growth processes and diffusion; reaction-diffusion models have been studied in ecology to describe the population dynamics of predator-prey model for a long time since Segel and Jackson applied Turing’s idea [29]. Since then, a new field of ecology, pattern formation, came into being. The problem of pattern and scale is the central problem in ecology, unifying population biology and ecosystems science and marrying basic and applied ecology [30]. The study of spatial patterns in the distribution of organisms is a central issue in ecology, geology, chemistry, physics, and so on [1, 3, 11, 15, 16, 25, 3156]. Theoretical work has shown that spatial and temporal pattern formation can play a very important role in ecological and evolutionary systems. Patterns can affect, for example, stability of ecosystems, the coexistence of species, invasion of mutants, and chaos. Moreover, the patterns themselves may interact, leading to selection on the level of patterns, interlocking eco-evolutionary time scales, evolutionary stagnation, and diversity.

Based on the above discussions, the spatially extended Holling-type IV predator-prey model with reaction diffusion takes as the form where and are the diffusion coefficients, respectively, and is the usual Laplacian operator in two-dimensional space; other parameters are the same definition as those in model (3).

It is easy to know that when a spatially extended predator-prey model is considered, the evolution of the model is decided by two sorts of sources (internal source and external source) which act together. The internal source is the dynamics of the individuals of the model, and the external source is the variability of environment. Some of the variability is periodic, such as temperature, water, food supply of the prey, and mating habits. It is necessary and important to consider models with periodic ecological parameters or perturbations which might be quite naturally exposed [57]. These periodic factors are regarded as external periodic forcing in the predator-prey systems. The external forcing can affect the population of predator and prey, respectively, which would go extinct in a deterministic environment. And some of the variability is irregular, such as the seasonal changes of the weather, food supply of the prey, and mating habits, and the effects of this variability are the so-called “noise.” Ecological population dynamics are inevitably “noisy” [2]. In the predator-prey systems, the random fluctuations also are undeniably arising from either environmental variability or internal species. To quantify the relationship between fluctuations and species’ concentration with spatial degrees of freedom, the consideration of these fluctuations supposes to deal with noisy quantities whose variance might at times be a sizable fraction of their mean levels. For example, the birth and death processes of individuals are intrinsically stochastic fluctuations which become especially pronounced when the number of individuals is small [16]. Moreover, there are many other stochastically factors causing predator-prey populations to change, such as effects of spatial structure of the habitat on the predator-prey ecosystem. The interactions between the predator and prey, which are far from being uniformly distributed, also introduce randomness. And these processes can be regarded as a parameter that fluctuates irregularly in space and time.

External forcing and noise induce effects in population dynamics, such as pattern formation, stochastic resonance, delayed extinction, enhanced stability, and quasiperiodic oscillations, which have been investigated with increasing interest in the past decades [16, 34, 48, 56, 5863]. And noise cannot systematically be neglected in models of population dynamics [63]. Zhou and Kurths [56] concluded these periodic variabilities as external forcing and investigated the interplay among noise, excitability, and mixing and external forcing in excitable media advected by a chaotic flow, in a two-dimensional FitzHugh-Nagumo model described by a set of reaction-advection-diffusion equations. And Si et al. [61] studied the propagation of traveling waves in subexcitable systems driven; Liu et al. [59] considered a spatially extended phytoplankton-zooplankton system with additive noise and periodic forcing. Following these models they considered, the Holling-type IV predator-prey model with external periodic forcing and colored noise is as follows: where denotes the periodic forcing with amplitude and angular frequency . The colored noise term is introduced additively in space and time, referring to the fluctuations in the predator death rate, which partially results from the environmental factors such as epidemics, weather, and nature disasters and it is the Ornstein-Uhlenbeck process that obeys the following linear stochastic partial differential equation: where is a Gaussian white noise or the so-called Markovian random telegraph process in both space and time with zero mean and correlation: where denotes averaging with respect to the noise and the Dirac delta-function and is the spatial correlation function of the Gaussian white noise .

Integrating (6) with respect to time , we get The mean value of the colored noise is and the correlation function of the colored noise is given by Let ; then

The colored noise generated in this way represents a simple spatiotemporal structured noise that can be used to real mimic situations, which is temporally correlated and white in space, satisfying where the temporal memory of the stochastic process is controlled by and is the intensity of noise. In this paper, we set .

Based on these discussions above, in this paper, we mainly focus on the spatiotemporal dynamics of models (4) and (5). And the organization is as follows. In Section 2, we employ the method of stability analysis to derive the symbolic conditions for Hopf and Turing bifurcation in the spatial domain. In Section 3, we give the complex dynamics of models (4) and (5), involving pattern formation, phase portraits, time-series plots and resonant response, and so on, via numerical simulation. Then, in the last section, we give some discussions and remarks.

2. Hopf and Turing Bifurcation

The nonspatial model (3) has at least two equilibria (steady states) which correspond to spatially homogeneous equilibria of models (4) and (5), in the positive quadrant: (total extinct) is a saddle; (extinct of the predator or prey-only) is a attracting node if , a saddle if , or a saddle-node if . When there exists unique stationary coexistence state , where On the other hand, when there exists another unique stationary coexistence state implying

It is worth mentioning that equilibria and cannot coexist. In this paper, we mainly focus on the dynamics of and rewrite it as . The dynamics behavior of is similar to that of .

To perform a linear stability analysis, we linearize model (3) around the stationary state for small space- and time-dependent fluctuations and expand them in Fourier space: where is the eigenvalue of the Jacobian matrix of model (3).

Hopf bifurcation is an instability induced by the transformation of the stability of a focus. Mathematically speaking, Hopf bifurcation occurs when and , at ; is the imaginary part, is the real part, and is the wave number. So we get the Hopf bifurcation surface: where the frequency of periodic oscillations in time satisfies , and the corresponding wavelength satisfies . In particular, we take as the bifurcation parameter and can get the critical value of Hopf bifurcation from (18):

Turing instability is induced only by “pursuit and evasion” if the predator can catch the prey by pursuit. We call the critical state of Turing instability as Turing bifurcation. Turing bifurcation occurs when “ and , at ,” and the wavenumber satisfies . In addition, at the Turing threshold, the spatial symmetry of the system is broken and the patterns are stationary in time and oscillatory in space with the wavelength . And the Turing bifurcation surface is given by where and the critical value of Turing bifurcation can be obtained from (21) as follows: where

Linear stability analysis yields the bifurcation diagram with , , , , , , and as shown in Figure 1(a). In this case, parameters , and is the unique stationary coexistence state. From Figure 1(a), one can see that the Hopf bifurcation line and the Turing bifurcation curve separate the parametric space into three distinct domains. In domain I, all two bifurcation lines are located below; the uniform steady state is the only stable solution of the model. Domain II is the region of pure Hopf instability. When the parameters correspond to domain III, which is located above all two bifurcation lines, both Hopf and Turing instability occur. Figure 1(b) illustrates the relation between the real and the imaginary parts of the eigenvalue with , which is located in domain II; one can see that when , and . Figure 1(c) displays the case of the critical value of Turing bifurcation ; in this case, and at . When , parameters are located in domain III; Figure 1(d) indicates that, at , and .

3. Spatiotemporal Dynamics of the Models

In this section, we perform extensive numerical simulations of the spatially extended models (4) and (5) in two-dimensional space, and the qualitative results are shown here. All our numerical simulations employ the zero-flux Neumann boundary conditions with a system size of space units. The parameters are , , , , , , , , and or , which satisfy . Models (4) and (5) are 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 a finite difference approximation for model (4) or Fourier transform method for model (5) for the spatial derivatives and an explicit Euler method for the time integration with a time stepsize of and space stepsize (lattice constant) of . When the system reached a periodic oscillatory 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 we have 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, for instance, the distribution of prey .

3.1. Pattern Formation of Model (4)

Figure 2 shows the evolution of the spatial patterns of prey at , 100, 300, 500, 1000, and 2000, with random small perturbation of the equilibrium of model (4) with , located in domain II, more than the Hopf bifurcation threshold and less than the Turing bifurcation threshold . In this case, pure Hopf instability occurs. One can see that, for model (4), the random initial distribution (cf. Figure 2(a)) leads to the formation of macroscopic spiral patterns (cf. Figures 2(d) to 2(f)). In other words, in this situation, spatially uniform steady-state predator-prey coexistence no longer exists. Small random fluctuations will be strongly amplified by diffusion, leading to nonuniform population distributions. From the analysis in Section 2, we find that, with these parameters in domain II, the spiral pattern arises from the Hopf instability. The lower panel in Figure 2 shows the corresponding (g) time series and (h) phase portraits. Figure 2(g) illustrates the evolution process of prey and periodic oscillating in time finally; (h) exhibits the fact that a limit cycle arises, which is caused by the Hopf bifurcation.

When , in this case, parameters in domain III (Figure 1(a)) and both Hopf and Turing instabilities occur. The nontrivial stationary state is . As an example, the formation of a regular macroscopic two-dimensional spatial pattern is shown in Figure 3. The lower panel in Figure 3 shows the corresponding (g) time-series plots and (h) phase portraits.

Comparing this situation (Figure 3) with the one above (Figure 2), it is easy to see that the pattern formations are all spiral wave. From the analysis in Section 2, we know that when , the wavelength while, at , . And the frequency of periodic oscillations in time is as inverse proportion with wavelength, so we can know that Turing instability has positive effect on the frequency while it has negative effect on wavelength. This is the reason why the spiral curves are more dense in Figure 3(f) than in Figure 2(f). On the other hand, one can see that when , the time-series plots (cf. Figure 3(g)) indicate that when Turing instability occurs, the solution of model (4) is strongly oscillatory in time while with (pure Hopf bifurcation emerges) it is periodic (cf. Figure 2(g)). In addition, comparing Figure 2(g) with Figure 3(g), one can see that Turing instability has positive effects on the amplitude of prey . And from Figure 3(h), one can see that a quasilimit cycle emerges while, in Figure 2(h), it is a cycle. Although there is some difference points between Figures 2 and 3, we can know that Turing instability cannot give birth to different type pattern. In our previous work [51], we find that Turing instability can change pattern type. This may be an important difference between the Holling-type IV and the ratio-dependent functional response of predator-prey model.

On the other hand, the basic idea of diffusion-driven instability in a reaction-diffusion system can be understood in terms of an activator-inhibitor system or predator-prey model (4). The functioning of this mechanism is based on three points [6]. First, a random increase of activator species (prey ) should have a positive effect on the creation rate of both activator (prey ) and inhibitor (prey ) species. Second, an increment in inhibitor species should have a negative effect on formation rate of both species. Finally, inhibitor species must diffuse faster than activator species . Certainly, the reaction-diffusion predator-prey model (4), with Holling-type IV functional response and predators diffusing faster than prey (i.e., ), provides this mechanism. And spirals and curves are the most fascinating clusters to emerge from the predator-prey model. A spiral will form from a wave front when the prey line (which is leading the front) overlaps the pursuing line of predator [38]. The prey on the extreme end of the line stops moving as there is no predator in their immediate vicinity. However the prey and the predator in the center of the line continue moving forward. This forms a small trail of prey at one (or both) end of the front. This prey starts breeding and the trailing line of prey thickens and attracts the attention of predator at the end of the fox line that turns towards this new source of prey. Thus a spiral forms with predator on the inside and prey on the outside. If the original overlap of prey occurs at both ends of the line a double spiral will form. Spirals can also form as prey blob collapses after predator eats into it. This is the reason why the pattern formation of model (4) is spiral wave.

3.2. The Effect of Noise Only of Model (5)

Now, we turn our focus on the effect of noise on the predator of stochastic model (5). In this case, ; that is, the periodic forcing is not present.

Figure 4 shows the dynamics of model (5) with noise on the predator. The first row of Figure 4, that is, (a), ; the second row, (b), ; the third row, (c), ; and the last row of Figure 4, (d), . And the first column of Figure 4, marked as (i), shows the snapshots of spatiotemporal pattern of model (5) at with different intensity of noise, respectively. In this case, one can see that the pattern formation turns into spatial chaotic from spiral wave with the increase of noise intensity . And the second column of Figure 4, marked as (ii), displays the phase portraits of model (5) with different intensity of noise, respectively. We can see that, as noise intensity increasing, the symmetry of the limit cycle is broken and gives rise to chaos. The last column of Figure 4, (iii), illustrates the time-series plots of prey with different intensity of noise, respectively. One can see that noise breaks the periodic oscillations in time and gives rise to drastically ruleless oscillations in time.

3.3. The Effect of Periodic Forcing of Model (5)

In the previous subsection, we have shown the effect of noise on the predator of model (5). An interesting question is whether such noise-sustained oscillations can be entrained by a weak external forcing, in this case, . This is investigated here.

When model (5) is noise free, there is a phenomenon of frequency locking or resonant response [56, 5861]. That is, without noise, the spatially homogeneous oscillation does not respond to the external periodic forcing when the amplitude is below a threshold whose value depends on the external period . Above the threshold, model (5) may produce oscillations about period with respect to external period , which is called frequency locking or resonant response. That is, the model produces one spike within each of the periods of the external force, called resonant response [56, 61]. The phenomenon of coherence resonance is of great importance [60]. Following Si et al. [61], in the present paper, the output period is defined as follows: is the time interval between the th spike and th spike. spikes are taken into account and the average value of them is .

As an example, with the amplitude , Figure 5 shows resonant response with (a) and (c), respectively. And Figures 5(b) and 5(d) are the phase portraits corresponding to (a) and (c). We can see that when , there exists a periodic orbit, while, , a periodic-2 orbit of model (5) emerges. Obviously, different can emerge from the same resonant response, and different phase orbits, that is, different numerical solution of model (5), may correspond to the same resonant response.

3.4. The Effect of Noise and Periodic Forcing of Model (5)

Now, we consider the dynamics about resonant response of model (5) with both noise and periodic forcing. As depicted in Figure 6, the prey can generate 5 : 1 (a) and 4 : 1 (c) locked oscillations, depending on the amplitude and angular frequency . Figures 6(b) and 6(d) illustrate the spiral pattern at corresponding to (a) and (c), respectively. In contrast, we change one of the parameters of Figure 6(c)   to (e); one can see that the resonant response vanishes and the corresponding spiral pattern (f) is similar to (b). It indicates that the amplitude is a control factor for pattern formation. In addition, comparing Figure 6(b) with 6(d), one can see that the pattern formations are determined by noise intensity , too.

In Figure 7, we have shown a typical pattern formation process in the 5 : 1 frequency locking regime with and . From (a) to (f), the pattern formation of prey is spiral wave and some small excitations already develop. One can see that, during the second period of the forcing, the prey is almost fully synchronized and relaxes slowly back to the state at moment (f). Obviously, the external periodic forcing at moment (e) repeats that at moment (a). However, the prey does not exactly repeat that due to a small fluctuation of the phase difference.

4. Conclusions and Remarks

In this paper, we present a spatial Holling-type IV predator-prey model containing some important factors, such as noise (random fluctuations), the external periodic forcing, and diffusion processes. And the numerical simulations were consistent with the predictions drawn from the bifurcation analysis, that is, Hopf bifurcation and Turing bifurcation.

If the parameter , the carrying capacity, is located in domain II of Figure 1(a), the Hopf instability occurs and the destruction of the pattern begins from the prey , while it begins from the predator if is located in domain III and both Hopf and Turing instabilities occur. From an ecological viewpoint, it shows that the initial and relatively rapid invasion of prey by predators can be followed by two subsequent invasions.

Furthermore, we demonstrate that noise and the external periodic forces play a key role in the predator-prey model (5) with the numerical simulations. We provoke qualitative transformations of the response of the model by changing noise intensity; noise can enhance the oscillation of the species density and format large clusters in the space. Periodic oscillations appear when the spatial noise and external periodic forcing are turned on; it also has been realized that model (5) is very sensitive to external periodic forcing through the natural annual variation of prey growth. In conclusion, we have shown that the cooperation between noise and external periodic forces inherent to the deterministic dynamics of periodically driven models gives rise to the appearance of a rich transport phenomenology.

Significantly, model (5) exhibits oscillations when both noise and external forces are present. This means that the dynamics of the predator population may be partly determined not only by the deterministic factors but also by the external forcing and the stochastic factors. Therefore, the model for spatially extended systems composed of two species could be useful to explain spatiotemporal behaviors of populations whose dynamics are strongly affected by noise and the environmental physical variables, and the results of this paper are an important step toward providing the theoretical biology community with simple practical numerical methods, for investigating the key dynamics of realistic predator-prey models.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This research was supported by NSFC no. 11071273.