Complex dynamics of a Holling-type IV predator-prey model
aa r X i v : . [ q - b i o . P E ] J a n Complex dynamics of a Holling-type IV predator-prey model
Lei Zhang,
1, 2
Weiming Wang, ∗ Yakui Xue, and Zhen Jin Institute of Nonlinear Analysis, School of Mathematics and Information Science,Wenzhou University, Wenzhou,Zhejiang, 325035 Department of Mathematics, North University of China,Taiyuan, Shan’xi 030051, P.R. China (Dated: October 29, 2018)
Abstract
In this paper, we focus on a spatial Holling-type IV predator-prey model which contains someimportant factors, such as diffusion, noise (random fluctuations) and external periodic forcing. Bya brief stability and bifurcation analysis, we arrive at the Hopf and Turing bifurcation surface andderive the symbolic conditions for Hopf and Turing bifurcation in the spatial domain. Based onthe stability and bifurcation analysis, we obtain spiral pattern formation via numerical simulation.Additionally, we study the model with colored noise and external periodic forcing. From thenumerical results, we know that noise or external periodic forcing can induce instability and enhancethe oscillation of the species, and resonant response. Our results show that modeling by reaction-diffusion equations is an appropriate tool for investigating fundamental mechanisms of complexspatiotemporal dynamics.
PACS numbers: 87.23.Cc, 82.40.Ck, 05.40.Ca, 47.54.-rKeywords: Functional response; Hopf bifurcation; Turing instability; Pattern formation; ∗ Electronic address: [email protected] . INTRODUCTION Predation, a complex natural phenomena, exists widely in the world, e.g., the sea, theplain, the forest, the desert and so on [7]. To model this phenomenon, predator-prey modelhas been suggested for a long time since the pioneer work of Lotka and Volterra [22]. Thepredator-prey model is a kind of “pursuit and evasion” system in which the prey try toevade the predator and the predator try to catch the prey if they interact [36]. Pursuitmeans the predator try to shorten the spatial distance between the predator and the prey,while evasion means the prey try to widen this spatial distance. In fact, the predator-preymodel is a mathematical method to approximate some part of our real world. And thedynamic behavior of the predator-prey model has long been and will continue to be one ofthe dominant themes in both ecology and mathematical ecology due to its universal existenceand importance [6, 25].Generally, a classical predator-prey model can be written as the form [3, 4]: dNdt = N f ( N ) − mP g ( N, P ) , dPdt = P [ cmg ( N, P ) − d ] , (1)where N and P stand for prey and predator quantity, respectively, f ( N ) the per capitarate of increase of the prey in absence of predation, d the food-independent death rate ofpredator, g ( N, P ) the functional response, the prey consumption rate by an average singlepredator, which obviously increases with the prey consumption rate and can be influenced bythe predator density, and which refers to the change in the density of prey attached per unittime per predator as the prey density changes, mg ( N, P ) the amount of prey consumed perpredator per unit time, and cmg ( N, P ) the predator production per capita with predation.In population dynamics, a functional response g ( N, P ) describes the relationship be-tween the predator and their prey, and the predator-prey model is always named after thecorresponding functional response for its key position [1, 3, 4, 42]. In the history of popula-tion ecology, both ecologists and mathematicians have a great interest in the Holling-typepredator-prey models [1, 8, 11, 12, 20, 24, 26, 34, 36, 39, 47, 55, 56, 58], including Hollingtype I-III, originally due to Holling [18, 19], and Holling type IV, suggested by Andrews [2].The Holling-type functional responses are so-called “prey-dependent” type [1], for g ( N, P )in Eq.(1) is a function only related to prey N . The classical expression of Holling-type IIfunctional response is g ( N, P ) = mN bN , and g ( N, P ) = mN aN is called Holling-type III. The2olling-type IV functional response is written as follows: g ( N, P ) = mN bN + aN . (2)Fuction (2) is called Monod-Haldane-type functional response too [38]. In addition, when b = 0, a simplified form g ( N, P ) = mN aN is proposed by Sokol and Howell [45], and somescholars also called it as Holling-type-IV [38, 42]. In this paper, we focus on the Holling-typeIV functional response which takes the form (2), and the corresponding predator-prey modelwhich takes the form: dNdt = rN (cid:0) − NK (cid:1) − mNP bN + aN , dPdt = P (cid:0) − q + cmN bN + aN (cid:1) , (3)where r > m > c > q > K > a > b > − √ a such that the denominator of above system does not vanish for non-negative N .On the other hand, the real world we lived in is a spatial world, and spatial patternsare ubiquitous in nature, which modify the temporal dynamics and stability properties ofpopulation density on range of spatial scales, whose effects must be incorporated in temporalecological models that do not represent space explicitly [37]. And the spatial component ofecological interactions has been identified as an important factor in how ecological commu-nities are shaped. Empirical evidence suggests that the spatial scale and structure of theenvironment can influence population interactions and the composition of communities [7].Reaction-diffusion model is a typical spatial extended model. It involves not only timebut also space and consists of several species which react with each other and diffuse withinthe spatial domain. It also involves a pair of partial differential equations, and represents thetime course of reacting and diffusing process. In the spatial extended predator-prey model,the interaction between the predator and the prey is the reaction item, and the diffusionitem comes into being for the predator’s “pursuit” and the prey’s “evasion”. Diffusion is aspatial process, and the whole model describes the evolution of the predator and the preygoing with time.Decades after Turing [48] demonstrated that spatial patterns could arise from the inter-action of reactions or growth processes and diffusion, reaction-diffusion models have been3tudied in ecology to describe the population dynamics of predator-prey model for a longtime since Segel and Jackson applied Turing’s idea [43]. Since then, a new field of ecology,pattern formation, came into being. The problem of pattern and scale is the central prob-lem in ecology, unifying population biology and ecosystem science, and marrying basic andapplied ecology [27]. And the study of spatial patterns in the distribution of organisms is acentral issue in ecology, geology, chemistry, physics and so on [5, 7, 9, 11, 13, 14, 15, 16, 17,21, 23, 26, 29, 30, 31, 32, 34, 35, 36, 38, 40, 41, 49, 51, 52, 53, 54, 57]. Theoretical work hasshown that spatial and temporal pattern formation can play a very important role in ecolog-ical and evolutionary systems. Patterns can affect, for example, stability of ecosystems, thecoexistence of species, invasion of mutants and chaos. Moreover, the patterns themselvesmay interact, leading to selection on the level of patterns, interlocking ecoevolutionary timescales, evolutionary stagnation and diversity.Based on the above discussions, the spatial extended Holling-type IV predator-prey modelwith reaction diffusion takes the form: ∂N∂t = rN (cid:0) − NK (cid:1) − mNP bN + aN + d ∇ N, ∂P∂t = P (cid:0) − q + cmN bN + aN (cid:1) + d ∇ P, (4)where d and d are the diffusion coefficients respectively, ∇ = ∂∂x + ∂∂y the usual Laplacianoperator in two-dimensional space, and other parameters are the same definition as those inmodel (3).Easy to know that, when a spatial extended predator-prey model is considered, the evo-lution 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, andthe external source is the variability of environment. Some of the variability is periodic, suchas temperature, water, food supply of the prey and mating habits. It is necessary and impor-tant to consider models with periodic ecological parameters or perturbations which mightbe quite naturally exposed [10]. These periodic factors are regarded as the external peri-odic forcing in the predator-prey systems. The external forcing can affect the population ofthe predator and prey, respectively, which would go extinct in a deterministic environment.And some of the variability is irregular, such as the the seasonal changes of the weather,food supply of the prey, mating habits, and the affects of this variability are called “noise”.Ecological population dynamics are inevitably “noisy” [22]. In the predator-prey systems,the random fluctuations are also undeniably arising from either environmental variability or4nternal species. To quantify the relationship between fluctuations and species’ concentra-tion with spatial degrees of freedom, the consideration of these fluctuations supposes to dealwith noisy quantities whose variance might at time be a sizable fraction of their mean levels.For example, the birth and death processes of individuals are intrinsically stochastic fluctu-ations which becomes especially pronounced when the number of individuals is small [34].Moreover, there are many other stochastically factors causing predator-prey population tochange, such as effects of spatial structure of the habitat on the predator-prey ecosystem.The interaction between the predator and prey, which are far from being uniformly dis-tributed, also introduce randomness. And these processes can be regarded as parameterfluctuates irregularly with spaces and time.The induced effects of the external forcing and noise in population dynamics, such aspattern formation, stochastic resonance, delayed extinction, enhanced stability, quasi peri-odic oscillations and so on, have been investigated with an increasing interest in the pastdecades [13, 28, 33, 34, 40, 44, 46, 50, 54, 57]. And noise cannot systematically be neglectedin models of population dynamics [50]. Zhou and Kurths [57] concluded these periodic vari-abilities as external forcing, and investigated the interplay among noise, excitability, mixingand external forcing in excitable media advected by a chaotic flow, in a two-dimensionalFitzHugh-Nagumo model described by a set of reaction-advection-diffusion equations. AndSi et al [44] studied the propagation of traveling waves in sub-excitable systems driven, andLiu et al [28] considered a spatial extended phytoplankton-zooplankton system with addi-tive noise and periodic forcing. Following these models they considered, the Holling-type IVpredator-prey model with external periodic forcing and colored noise is as follows: ∂N∂t = rN (cid:0) − NK (cid:1) − mNP bN + aN + A sin ( ωt ) + d ∇ N, ∂P∂t = P (cid:0) − q + cmN bN + aN (cid:1) + η ( r , t ) + d ∇ P, (5)where A sin( ωt ) denotes the periodic forcing with amplitude A and angular frequency ω . Thecolored noise term η ( r , t ) ( r = ( x, y ) ) is introduced additively in space and time, referring tothe fluctuations in the predator increase rate, which partially results from the environmentalfactors such as epidemics, weather and nature disasters, and it is the Ornstein-Uhlenbechprocess that obeys the following linear stochastic partial differential equation ∂η ( r , t ) ∂t = − τ η ( r , t ) + 1 τ ξ ( r , t ) , (6)5here ξ ( r , t ) is a Gaussian white noise or the so called Markovian random telegraph processin both space and time with zero mean and correlation: h ξ ( r , t ) i = 0 , h ξ ( r , t ) ξ ( r ′ , t ′ ) i = 2 εδ ( r − r ′ ) δ ( t − t ′ ) , where h·i denotes mean value with respect to the noise ξ ( r , t ), and δ the Dirac delta-function, δ ( r − r ′ ) the spatial correlation function of the Gaussian white noise ξ ( r , t ) .Integrating equation (6) with respect to time t , we get η ( r , t ) = η ( r , e − tτ + 1 τ e − tτ Z t e sτ ξ ( r , s ) ds. The mean value of the colored noise is h η ( r , t ) i = h η ( r , i e − tτ + 1 τ e − tτ Z t e sτ h ξ ( r , s ) i ds = h η ( r , i e − tτ , and the correlation function of the colored noise is given by h η ( r , t ) η ( r ′ , t ′ ) i = h η ( r , ih η ( r ′ , i e − t + t ′ τ + 1 τ e − t + t ′ τ Z t Z t ′ e s + s ′ τ h ξ ( r , s ) ξ ( r ′ , s ′ ) i dsds ′ (7)= h η ( r , ih η ( r ′ , i e − t + t ′ τ + ετ e − t + t ′ τ δ ( r − r ′ ) Z t Z t ′ e s + s ′ τ δ ( t − t ′ ) dsds ′ (8)= h η ( r , ih η ( r ′ , i e − t + t ′ τ + ετ δ ( r − r ′ )( e − t + t ′ τ − e − tτ + e − t − t ′ τ ) . (9)Let t → + ∞ , then h η ( r , t ) η ( r ′ , t ′ ) i → ετ e − t − t ′ τ δ ( r − r ′ ) . The colored noise η ( r , t ) generated in this way represents a simple spatiotemporal struc-tured noise that can be used in real mimic situations, which is temporally correlated andwhite in space, and satisfies h η ( r , t ) η ( r ′ , t ′ ) i = ετ e − | t − t ′| τ δ ( r − r ′ ) , where the temporal memory of the stochastic process controlled by τ and ε is the intensityof noise. In this paper, we set τ = 1.Based on these discussion above, in this paper, we mainly focus on the spatiotemporaldynamics of model (4) and (5). And the organization is as follows: In section 2, we employthe method of stability analysis to derive the symbolic conditions for Hopf and Turingbifurcation in the spatial domain. In section 3, we give the complex dynamics of model (4)and (5), involving pattern formation, phase portraits, time series plots and resonant responseand so on, via numerical simulation. Then, in the last section, we give some discussions andremarks. 6 I. HOPF AND TURING BIFURCATION
The non-spatial model (3) has at least two equilibria (steady states) which correspond tospatially homogeneous equilibria of the model (4) and (5), in the positive quadrant: (0 , K,
0) (extinct of the predator, or prey-only) is a attracting nodeif q > cmK Kb + aK , a saddle if q < cmK Kb + aK or a saddle-node if q = cmK Kb + aK . When( a, b, c, m, q, r, K ) ∈ E , here, E = { ( a, b, c, m, q, r, K ) | mc > qb, q a < ( mc − qb ) , p ( mc − qb ) − q a > − m c +2 qbmc + qaKmc − q aKb +2 q a − q b − mc + qb + qaK > , ba − mcqa + K < } , there exists unique stationary coexistent state ( N ∗ , P ∗ ), where N ∗ = − qb + mc − Aqa , P ∗ = cr (( − mc + bq + qaK ) N ∗ + q ) aq K . On the other hand, when( a, b, c, m, q, r, K ) ∈ E , here, E = { ( a, b, c, m, q, r, K ) | mc > qb, q a < ( mc − qb ) , p ( mc − qb ) − q a > − − m c +2 qbmc + qaKmc − q aKb +2 q a − q b − mc + qb + qaK > , ba − mcqa + K > } , there exists another unique stationary coexistent state ( N ∗ , P ∗ ) implying: N ∗ = − qb + mc + Aqa , P ∗ = cr (( − mc + bq + qaK ) N ∗ + q ) aq K . It is worth mentioning that equilibria ( N ∗ , P ∗ ) and ( N ∗ , P ∗ ) cannot coexist. In thispaper, we mainly focus on the dynamics of ( N ∗ , P ∗ ) and rewrite it as ( N ∗ , P ∗ ). The dynamicbehavior of ( N ∗ , P ∗ ) is similar to those of ( N ∗ , P ∗ ).To perform a linear stability analysis, we linearize model (3) around the stationary state( N ∗ , P ∗ ) for small space- and time-dependent fluctuations and expand them in Fourier space N ( r , t ) ∼ N ∗ e λt e i~k · r , P ( r , t ) ∼ P ∗ e λt e i~k · r , r = ( x, y ) , ~k = ( k x , k y ) . 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 Im( λ ) = 0 , Re( λ ) = 0 at k = 0,where Im( λ ) is the imaginary part, Re( λ ) the real part and k the wave number. So we getthe Hopf bifurcation surface H = { ( a, b, c, m, q, r, K ) | det( J ) > , trace( J ) = 0 } , (10)7here det( J ) = − (cid:0) r − rN ∗ K (cid:1) q + mqP ∗ + cm ( r − rN ∗ K ) N ∗ ( bN ∗ + aN ∗ ) − mqN ∗ P ∗ ( b +2 aN ∗ ) ( bN ∗ + aN ∗ ) , trace( J ) = r − rN ∗ K − q + m ( − P ∗ + cN ∗ + aN ∗ P ∗ + bcN ∗ + cN ∗ a )( bN ∗ + aN ∗ ) , the frequency of periodic oscillations in time ω H satisfies ω H = Im( λ ) = p det( J ), and thecorresponding wavelength λ H satisfies λ H = πω H = π √ det( J ) . Especially, we take K as thebifurcation parameter, and can get the critical value of Hopf bifurcation from Eq. (10): K H = ( − ( aq (5 mc − qb ) − (3 mc − qb )( mc − qb ) ) p ( mc − qb ) − q a − q a + q ( mc − qb )(11 mc − qb ) a − (3 mc − qb )( mc − qb ) ) / ( − aq (((2 mc − qb )( mc − qb ) − q a ) p ( mc − qb ) − q a − aq (3 mc − qb ) + (2 mc − qb )( mc − qb ) )) . (11)Turing instability is induced only by “pursuit and evasion” if the predator can catch theprey by pursuit. We call the critical state of Turing instability as Turing bifurcation. Turingbifurcation occurs when Im( λ ) = 0 , Re( λ ) = 0 at k = k T = 0 , and the wavenumber k T satisfies k T = q det( J ) d d . In addition, at the Turing threshold, the spatial symmetry of thesystem is broken and the patterns are stationary in time and oscillatory in space with thewavelength λ T = πk T . And the Turing bifurcation surface implies: T = { ( a, b, c, m, q, r, d , d , K ) | det( J k ) = 0 , trace( J k ) = 0 } , (12)where det( J k ) = − (cid:0) r − rN ∗ K − d k (cid:1) ( q + d k ) (( q + d k ) mP ∗ + ( r − rN ∗ K − d k ) cmN ∗ ) bN ∗ + aN ∗ − m ( b +2 aN ∗ ) ( q + d k ) N ∗ P ∗ (1+ bN ∗ + aN ∗ ) , trace( J k ) = r − rN ∗ K − q − ( d + d ) k + m ( − P ∗ + cN ∗ + aN ∗ P ∗ + bcN ∗ + cN ∗ a )( bN ∗ + aN ∗ ) , the critical value of Turing bifurcation can be obtained from Eq. (12) as follows: K T = F F , (13)where F = r ((4 q a (2 mc − qb ) − (3 mc − qb )( mc − qb ) ) p ( mc − qb ) − q a + (( mc − qb ) − q a ) · ((3 mc − qb )( mc − qb ) − q a ))((3 mc − qb )( mc − qb ) − q a − (3 mc + qb ) p ( mc − qb ) − q a ) d , = qa (2 mc (( mc − qb ) p ( mc − qb ) − q a + 4 q a − ( mc − qb ) ) B + ((3 mc − qb ) d · (2 m c − qbmc + q b − q a ) r + 2 qd mc ( mc − qb ) )(( mc − qb ) − q a )+( − d ((3 mc − qb )(2 mc − qb )( mc − qb ) − q a ( − q a + 3 q b − qbmc +11 m c )) r − qcm ( mc − qb ) d (( mc − qb ) − q a )) p ( mc − qb ) − q a +( q b − qbmc + m c − q a )((2 mc − qb ) d (3 m c − qbmc + q b − q a ) r +2 mcqd (( mc − qb ) − q a ))) ,B = ( − d q ((( d bq − qd mc − qrbd + 2 rmcd )(4 q a − ( mc − qb ) ) − rmcd ( mc − qb ) ) · p ( mc − qb ) − q a + (8 q d r − q d ) a + 4 ( mc − qb ) aq rmcd + (3 rmcd − qd mc + d bq − qrbd )(( mc − qb ) − mc − qb ) aq ))) . Linear stability analysis yields the bifurcation diagram with r = 1, a = 0 . b = 1, c = 0 . m = 0 . q = 0 .
18 and d = 0 . a, b, c, m, q, r, K ) ∈ E , and ( N ∗ , P ∗ ) is the unique stationary coexistent state. FromFigure 1(a), one can see that the Hopf bifurcation line and the Turing bifurcation curveseparate the parametric space into three distinct domains. In domain I, located below alltwo bifurcation lines, the steady state is the only stable solution of the model. Domain II isthe region of pure Hopf instability. When the parameters correspond to domain III, locatedabove all two bifurcation lines, both Hopf and Turing instability occur. Figure 1(b) displaysthe dispersion relations showing unstable Hopf mode, transition of Turing mode from stableto unstable for model (4), e.g., as d decreased. Figure 1(c) illustrates the relation betweenthe real and the imaginary parts of the eigenvalue λ with K = 2 . > K H = 2 . k = 0, Re( λ ( k )) > λ ( k )) = 0. Figure 1(d)displays the case of the critical value of Turing bifurcation K = K T = 3 . λ ( k )) = 0 and Im( λ ( k )) = 0 at k = k T = 2 . K = 4 .
0, located in domain III,figure 1(e) indicates that at k = 0, Re( λ ( k )) >
0, Im( λ ( k )) = 0. III. SPATIOTEMPORAL DYNAMICS OF THE MODELS
In this section, we perform extensive numerical simulations of the spatially extendedmodel (4) and (5) in two-dimensional spaces, and the qualitative results are shown here. Allour numerical simulations employ the Zero-flux Neumann boundary conditions with a system9
IG. 1: (a) K − d Bifurcation diagram for model (4) with r = 1, a = 0 . b = 1, c = 0 . m = 0 . q = 0 .
18 and d = 0 .
2. Hopf and Turing bifurcation lines separate the parameterspace into three domains. P (0 . , . P (0 . , . P (0 . , .
0) are the different selectionscorresponding to (c), (d), (e), respectively. (b) Dispersion relations showing unstable Hopf mode,transition of Turing mode from stable to unstable for the system model (4) as d decreased. Theother parameters in (c)–(e) are: d = 0 .
02, the bifurcation parameter K equals: (c) 2 . > K H =2 . .
499 = K T ; (e) 4 . > K T > K H . The real parts Re( λ ) and the imaginary parts Im( λ )are shown by solid curves and dashed curves, respectively. size of 200 ×
200 space units. The parameters are r = 1, a = 0 . b = 1, c = 0 . m = 0 . q = 0 . d = 0 . d = 0 . K = 2 . K = 4 .
0, which satisfy ( a, b, c, m, q, r, K ) ∈ E .Model (4) and (5) are integrated initially in two-dimensional space from the homogeneoussteady state, i.e., we start with the unstable uniform solution ( N ∗ , P ∗ ) with small randomperturbation superimposed, in each, the initial condition is always a small amplitude randomperturbation ( ± × − ), using a finite difference approximation for model (4) or Fouriertransform method for model (5) for the spatial derivatives and an explicit Euler method forthe time integration with a time stepsize of ∆ t = 1 /
24 and space stepsize (lattice constant)∆ x = ∆ y = 1. We have taken some snapshots with white (black) corresponding to the high(low) value of prey N . 10n the numerical simulations, different types of dynamics are observed and we have foundthat the distributions of the predator and prey are always of the same type. Consequently,we can restrict our analysis of pattern formation to one distribution. In this section, weshow the distribution of prey N , for instance. A. Pattern formation of model (4)
Figure 2 shows the evolution of the spatial patterns of prey N at t = 0, 100, 300, 500,1000, 2000, with random small perturbation of the equilibrium ( N ∗ , P ∗ ) = (0 . , . K = 2 .
8, located in domain II, more than the Hopf bifurcation threshold K H = 2 .
279 and less than the Turing bifurcation threshold K T = 3 . N , periodic oscillating in time finally, (h) exhibitsthat a limit cycle arises, which is caused by the Hopf bifurcation.When K = 4 . > K T > K H , in this case, parameters in domain III (figure 1(A)),both Hopf and Turing instabilities occur. The nontrivial stationary state is ( N ∗ , P ∗ ) =(0 . , . K = 2 .
8, the wavelength λ = 3 .
100 while at K = 4 . λ = 3 . K = 4 .
0, the time series plots (c.f., figure 3(g))11
IG. 2: Grey-scaled snapshots of spatiotemporal pattern of the prey N of model (4) with K = 2 . t = 0, (b) t = 100, (c) t = 300, (d) t = 500, (e) t = 1000, (f) t = 2000. The lower panels showthe corresponding (g) time series plots and (h) phase portraits. indicate that when Turing instability occurs, the solution of model (4) is strongly oscillatoryin time while with K = 2 . N . And from figure 3(h), one can see that a quasilimit cycle emerges while in figure 2(h), it is a cycle. Although there are some different pointsbetween figure 2 and figure 3, but we can know that Turing instability can’t give birth todifferent type patterns. In our previous work [51], we find that Turing instability can changepattern type. This may be a very point between the Holling-type IV with Michael-Mentonfunctional response of predator-prey model.On the other hand, the basic idea of diffusion-driven instability in a reaction-diffusion sys-tem can be understood in terms of an activator-inhibitor system or predator-prey model (4).12 IG. 3: Grey-scaled snapshots of spatiotemporal pattern of the prey N of model (4) with K = 4 . t = 0, (b) t = 100, (c) t = 300, (d) t = 500, (e) t = 1000, (f) t = 2000. The lower panels showthe corresponding (g) time series plots and (h) phase portraits. The functioning of this mechanism is based on three points [3]. First, a random increase ofactivator species (prey, N ) should have a positive effect on the creation rate of both acti-vator (prey, N ) and inhibitor (prey, P ) species. Second, an increment in inhibitor speciesshould have a negative effect on formation rate of both species. Finally, inhibitor species P must diffuse faster than activator species N . Certainly, the reaction-diffusion predator-preymodel (4), with Holling-type IV functional response and predators diffusing faster than prey(i.e., d > d ), provides this mechanism. And spirals and curves are the most fascinatingclusters to emerge from the predator-prey model. A spiral will form from a wave front whenthe prey line (which is leading the front) overlaps the pursuing line of predator [17]. Theprey on the extreme end of the line stops moving as there is no predator in the immediatevicinity. However, the prey N and the predator P in the center of the line continue moving13orward. This forms a small trail of prey at one (or both) ends of the front. These preystart breeding and the trailing line of prey thickens and attracts the attention of predatorat the end of the fox line that turn towards this new source of prey. Thus a spiral formswith predator P on the inside and prey N on the outside. If the original overlap of preyoccurs at both ends of the line a double spiral will form. Spirals can also form as a preyblob collapses after the predator eats into it. This is the reason why the pattern formationof model (4) is spiral wave. B. The effect of noise only of model (5)
Now, we turn our focus on the effect of noise on the predator P of stochastic model (5).In this case, A = 0, i.e., the periodic forcing is not at presence.Figure 4 shows the dynamics of model (5) with noise on the predator. The first row offigure 4, i.e., (a), ε = 0 . ε = 0 .
01; the third row, (c), ε = 0 .
05; andthe last row of figure 4, (d), ε = 0 .
1. And the first collum of figure 4, marked as (i), showsthe snapshots of spatiotemporal pattern of model (5) at t = 2000 with different intensity ofnoise, respectively. In this case, one can see that the pattern formation turns into spatialchaotic from spiral wave with the increase of noise intensity ε . And the second collum offigure 4, marked as (ii), displays the phase portraits of model (5) with different intensity ofnoise, respectively. We can see that, as noise intensity ε increasing, the symmetry of thelimit cycle is broken and gives rise to chaos. The last collum of figure 4, (iii), illustrates thetime-series plots of prey N with different intensity of noise, respectively. One can see thatnoise breaks the periodic oscillations in time and gives rise to drastically ruleless oscillationsin time. C. The effect of periodic forcing of model (5)
In the previous subsection, we have shown the effect of noise on the predator P ofmodel (5). An interesting question is whether such noise-sustained oscillations can be en-trained by a weak external forcing, in this case, ε = 0, which is investigated here.When model (5) is noise free, there is a phenomenon of frequency locking or resonantresponse [28, 33, 44, 57]. That is, without noise, the spatially homogeneous oscillation does14 IG. 4: Dynamics of model (4), for the following noise intensity. (a) ε = 0 . ε = 0 . ε = 0 .
05; (d) ε = 0 .
1. (i) snapshots of pattern formation at time 2000; (ii) phase portraits; (iii)time-series plots. A = 0 and the other parameters are the same as those in figure 2. not respond to the external periodic forcing when the amplitude A is below a thresholdwhose value depends on the external periodic T in = πω . Above the threshold, model (5)may produce oscillations about period T out with respect to external period T in , which iscalled frequency locking or resonant response. That is, the model produces one spike withineach of the M = T out T in periods of the external force, called M : 1 resonant response [44, 57].The phenomenon of coherent resonance is of great importance [33]. Following Si [44], in thepresent paper, the output periodic T out is defined as follows: T i is the time interval betweenthe i th spike and ( i + 1)th spike. m spikes are taken into account and the average value of15hem is T out = m P i =1 T i / ( m − A = 0 . ω = 0 . π (a) and ω = 0 . π (c), respectively. And figure 5(b) and (d) are the phase portraitscorresponding to (a) and (c). We can see that when ω = 0 . π , there exists a periodic orbit,while ω = 0 . π , a periodic-2 orbit of model (5) emerges. Obviously, different ω can formthe same resonant response, and different phase orbits, i.e., different numerical solution ofmodel (5), may correspond to the same resonant response. FIG. 5: External periodic forcing induced frequency locking of model(4). The solid curve is timeseries of prey U , the dash curve is the corresponding external periodic forcing. Other parametersare the same as those in figure 3. D. The effect of noise and periodic forcing of model (5)
Now, we consider the dynamics about resonant response of model (5) with both noiseand periodic forcing. As depicted in figure 6, the prey can generate 5 : 1 (a), 4 : 1 (c) lockedoscillations, depending on the amplitude A and angular frequency ω . Figures 6 (b) and16d) illustrate the spiral pattern at t = 2000 corresponding to (a) and (c), respectively. Forcontrast, we change one of the parameters of figure 6(c) A = 0 .
001 to A = 0 .
01 (e), one cansee that the resonant response vanishes, the corresponding spiral pattern (f) is similar to(b). It indicates that the amplitude A is a control factor for pattern formation. In addition,comparing figures 6(b) with (d), one can see that the pattern formations are determined bynoise intensity ε , too. FIG. 6: Dynamics of model (5) with both noise and periodic forcing. (b)(d)(f) are snapshots at t = 2000 corresponding to the left hand side resonant response. The other parameters are the sameas those in figure 3. In Figure 7, we have shown a typical pattern formation process in the 5 : 1 frequencylocking regime with A = 0 .
001 and ω = 0 . π . From t = 1870 (a) to t = 1920 (f), the pattern17 IG. 7: Typical pattern formation of the forced noisy prey in the 5 : 1 locking region at A = 0 . ε = 0 . N (the solid curve) and the corresponding external periodic forcing (the dash curve) correspondingto the snapshots of the patterns. The grey scale, from black to white, is in [0 , .
5] for all thesnapshots. formation of prey N is spiral wave and some small excitations already develop. One can seethat, during the second periodic of the forcing, the prey is almost fully synchronized andrelaxes slowly back to the state at moment (f). Obviously, the external periodic forcing atmoment (e) repeats that at moment (a). However, the prey N does not exactly repeat dueto a small fluctuation of the phase difference. IV. CONCLUSIONS AND REMARKS
In this paper, we present spatial Holling-type IV predator-prey model containing someimportant factors, such as noise (random fluctuations), the external periodic forcing and18iffusion processes. And the numerical simulations were consistent with the predictionsdrawn from the bifurcation analysis, that is, Hopf bifurcation and Turing bifurcation.If the parameter K , the carrying capacity, located in the domain II of figure 1(a), the Hopfinstability occurs, the destruction of the pattern begins from the prey N , while it beginsfrom the predator P if K located in the domain III, both Hopf and Turing instabilitiesoccur.Furthermore, we demonstrate that noise and the external periodic forcing play a keyrole in the predator-prey model (5) with the numerical simulations. We provoke qualitativetransformations of the response of the model by changing noise intensity, and noise canenhance the oscillation of the species density, and format large clusters in the space. Theperiodic oscillations appear when the spatial noise and external periodic forcing are turnedon. It has also been realized that model (5) is very sensitive to external periodic forcingthrough the natural annual variation of prey growth. In conclusion, we have shown thatthe cooperation between noise and external periodic forcing inherent to the deterministicdynamics of periodically driven models gives rise to the appearance of resonant response.Significantly, model (5) exhibits oscillations when both noise and external forcing arepresent. This means that the predator population may be partly due to the external forcingand stochastic factors instead of deterministic factors. Therefore, the model for spatiallyextended systems composed by two species could be useful to explain spatiotemporal be-havior of populations whose dynamics is strongly affected by noise and the environmentalphysical variables, and the results of this paper are an important step toward providing thetheoretical biology community with simple practical numerical methods, for investigatingthe key dynamics of realistic predator-prey models. Acknowledgments
This work was supported by the National Natural Science Foundation of China(10471040) and the Youth Science Foundation of Shanxi Provence (20041004). [1] Abrams P. A. and Ginzburg L. R., 2000. The nature of predation: prey dependent, ratiodependent or neither? Trends Ecol. Evol. 15, 337-341.
2] Andrews J. F., 1968. A mathematical model for the continuous culture of microorganismsutilizing inhibitory substrates, Biot. Bioe. 10, 707-723.[3] Alonso D., Bartumeus F. and Catalan J., 2002. Mutual interference between predators cangive rise to Turing spatial patterns, Ecology 83, 28-34.[4] Arditi R. and Ginzburg L. R., 1989. Coupling in predator-prey dynamics: ratio-dependence,J. Theo. Biol. 139, 311-326.[5] Baurmann M., Gross T. and Feudel U., 2007. Instabilities in spatially extended predator-preysystems: Spatio-temporal patterns in the neighborhood of Turing-Hopf bifurcations, J. Theo.Biol. 245, 220-229.[6] Berryman A., 1992. The origins and evolution of predator-prey theory, Ecology 75, 1530-1535.[7] Cantrell R., and Cosner C., 2003. Spatial Ecology via Reaction-Diffusion Equations, JohnWiley & Sons, Ltd., Chichester, England.[8] Chen Y., 2004. Multiple periodic solutions of delayed predatorCprey systems with type IVfunctional responses, Nonl. Anal.: Real World Appl. 5, 45-53.[9] Cross M. C. and Hohenberg P. C., 1993. Pattern formation outside of equilibrium, Rev. Mod.Phys. 65, 851-1112.[10] Cushing J., 1977. Periodic time-dependent predator-prey system, SIAM J. Appl. Math 32,82-95.[11] Estep D. and Neckels D., 2007. Fast methods for determining the evolution of uncertainparameters in reaction-diffusion equations, Computer Methods in Applied Mechanics andEngineering 196, 3967-3979.[12] Gakkhar S. and Singh B., 2007. The dynamics of a food web consisting of two preys and aharvesting predator, Chaos, Solitons and Fractals 34, 1346-1356.[13] Garcia-Ojalvo J. and Schimansky-Geier L., 1999. Noise-induced spiral dynamics in excitablemedia, Europhys. Lett. 47, 298-303.[14] Garvie M., 2007. Finite-difference schemes for reaction-diffusion equations modelling predator-prey interactions in matlab, Bull. Math. Biol. 69, 931-956.[15] Gierer A. and Meinhardt H., 1972. A theory of biological pattern formation, Kybernetik 12,30-39.[16] Griffith D. A. and Peres-Neto P. R., 2006. Spatial modeling in ecology: the flexibility ofeigenfunction spatial analyses, Ecology 87, 2603-2613.
17] Hawick, K., James, H. and Scogings, C., 2006. A zoology of emergent life patterns in apredator-prey simulation model, Technical Note CSTN-015 and in Proc. IASTED Interna-tional Conference on Modelling, Simulation and Optimization, September 2006, Gabarone,Botswana. 507-115.[18] Holling C. S., 1959. The components of predation as revealed by a study of small mammalpredation of the european pine sawfly, Cana. Ento. 91, 293-320.[19] Holling C. S., 1959. Some characteristics of simple types of predation and parasitism, Cana.Ento. 91, 385-395.[20] Huang J. and Xiao D., 2004. Analyses of Bifurcations and Stability in a Predator-prey Systemwith Holling Type-IV Functional Response, Acta Mathematicae Applicatae Sinica, EnglishSeries 20, 167-178.[21] Katsuragi H., 2006. Diffusion-induced spontaneous pattern formation on gelation surfaces.Europhys. Lett. 73, 793-799.[22] Kendall B., 2001. Cycles, chaos, and noise in predator-prey dynamics, Chaos Solitons Fractals12, 321-332.[23] Klausmeier C.A., 1999. Regular and Irregular Patterns in Semiarid Vegetation, Science 284,1826-1828.[24] Ko W. and Ryu K., 2007. Coexistence states of a predator-prey system with non-monotonicfunctional response, Nonl. Anal.: Real World Appl. 8, 769-786.[25] Kuang Y. and Beretta E., 1998. Global qualitative analysis of a ratio-dependent predator-preysystem, J. Math. Biol. 36, 389-406.[26] Lepp¨anen T., 2004. Coputational studies of pattern formation in Turing systems, PhD-Thesis,Helsinki University of Technology.[27] Levin S. A., 1992. The Problem of Pattern and Scale in Ecology, Ecology 73, 1943-1967.[28] Liu Q., Li B., Jin Z., 2007. Resonant patterns and frequency-locked induced by additive noiseand periodically forced in phytoplankton-zooplankton system, Preprint: arXiv:0705.3724.[29] Li Z., Gao M., Hui C., Han X. and Shi H., 2005. Impact of predator pursuit and prey evasionon synchrony and spatial patterns in metapopulation, Ecol. Model. 185, 245-254.[30] Madzvamuse A., Thomas R. D. K. , Maini P. K. and Wathen A. J., 2002. A NumericalApproach to the Study of Spatial Pattern Formation in the Ligaments of Arcoid Bivalves,Bull. Math. Biol. 64, 501-530.
31] Maini P. K., Baker R. E. and Chuong C., 2006. The Turing model comes of molecular age,Science 314, 1397-1398.[32] Maionchi D. O., Reis S. F. and Aguiar M. A. M., 2006. Chaos and pattern formation in aspatial tritrophic food chain, Ecol. Model. 191, 291-303.[33] Mankin R., Laas T., Sauga A. and Ainsaar A., 2006. Colored-noise-induced Hopf bifurcationsin predator-prey communities, Phys. Rev. E 74, 021101.[34] Malchow H., Hilker F. M., Petrovskii S. V., 2004. Noise and productivity dependence ofspatiotemporal pattern formation in a prey-predator system, Disc. Cont. Dyna. Syst. B 4,705-711.[35] Medvinsky A. B., Petrovskii S. V., Tikhonova I. A., Malchow H. and Li B.-L., 2002. Spa-tiotemporal Complexity of Plankton and Fish Dynamics, SIAM Review 44, 311-370.[36] Murray J., 2003. Mathematical biology. II. Spatial models and biomedical applications, 3rdedn., Interdisciplinary Applied Mathematics 18, Springer, New York.[37] Neuhauser C., 2001. Mathematical Challenges in Spatial Ecology, Notices of the AmericanMathematical Society 47, 1304-1314.[38] Pang P. Y. H. and Wang M. X., 2004. Non-constant positive steady positive steady statesofa predator-prey system with non-monotonic functional responctionlresponse and diffusion,Proc. Lond. Math. Soc. 88, 135-157.[39] Page K., Maini P. K. and Monk N. A. M., 2003. Pattern formation in spatially heterogeneousTuring reaction-diffusion models, Phys. D 181, 80-101.[40] Riaz S. S., Dutta S., Kar S. and Ray D. S., 2005. Pattern formation induced by additive noise:a moment-based analysis, Eur. Phys. J. B 47, 255-263.[41] Riaz S. S., Banarjee S., Kar S. and Ray D. S., 2006. Pattern formation in reaction-diffusionsystem in crossed electric and magnetic fields, Eur. Phys. J. B 53, 509-515.[42] Ruan S. G. and Xiao D. G., 2001. Global Analysis in a Predator-Prey System with Nonmono-tonic Functional Response, SIAM J. Appl. Math. 61, 1445-1472.[43] Segel L. A. and Jackson J. L., 1972. Dissipative structure: An explanation and an ecologicalexample, J. Theo. Biol. 37, 545-559.[44] Si F., Liu Q., Zhang J. and Zhou L., 2008. Propagation of travelling waves in sub-excitablesystems driven by noise and periodic forcing, Eur. Phys. J. B 60, 507-513[45] Sokol W. and Howell J. A., 1987. Kinetics of phenol oxidation by washed cell, Biot. Bioe. 30,21-927.[46] Spagnolo B., Valenti D. and Fiasconaro A., 2004. Noise in ecosystems: a short review, Math.Bios. Engi. 1, 185-211.[47] Skalski G. T. and Gilliam J. F., 2001. Functional responses with predator interference: Viablealternatives to the Holling type II model, Ecology 82, 3083-3092.[48] Turing A. M., 1952. The chemical basis of morphogenisis, Phil. Tran. R. Soc. London B 237,7-72.[49] Uriu K. and Iwasa Y., 2007. Turing Pattern Formation with Two Kinds of Cells and a DiffusiveChemical, Bull. Math. Biol. 69, 2515-2536.[50] Vilar J. M. G. and Sole R. V., 1998. Effects of Noise in Symmetric Two-Species Competition,Phys. Rev. Lett. 80, 4099-4102.[51] Wang W., Liu Q. and Jin Z., 2007. Spatiotemporal complexity of a ratio-dependent predator-prey system, Phys. Rev. E 75, 051913.[52] Yang L., Dolnik M., Zhabotinsky A. and Epstein, I., 2002. Pattern formation arising frominteractions between turing and wave instabilities, J. Chem. Phys. 117, 7259-7265.[53] Yang L., Milos D., Zhabotinsky A. M. and Epstein I. R., 2002. Spatial Resonances and Super-position Patterns in a Reaction-Diffusion Model with Interacting Turing Modes, Phys. Rev.Lett. 88, 208303.[54] Yang L., Zhabotinsky A. M. and Epstein I. R., 2004. Stable Squares and Other OscillatoryTuring Patterns in a Reaction-Diffusion Model, Phys. Rev. Lett. 92, 198303.[55] Zhang S., Tan D. and Chen L., 2006. Chaos in periodically forced Holling type IV predator-Cprey system with impulsive perturbations, Chaos, Soli. Frac. 27, 980-990.[56] Zhang W., Zhu D. and Bi P., 2007. Multiple positive periodic solutions of a delayed discretepredator-prey system with type IV functional responses, Appl. Math. Lett. 20, 1031-1038.[57] Zhou C., Kurths J., 2005. Noise-sustained and controlled synchronization of stirred excitablemedia by external forcing, New J. Phys. 7, 18.[58] Zhu H., Campbell S. A. and Wolkowicz G., 2002. Bifurcation analysis of a predator-preysystem with nonmonotonic functional response, SIAM J. Appl. Math. 63, 636-682.