Spatial Complexity of a Predator-Prey Model with Holling-Type Response

and Applied Analysis 3 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, 58–63]. 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 excitablemedia advected by a chaotic flow, in a twodimensional 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 thesemodels they considered, the Holling-type IV predator-prey model with external periodic forcing and colored noise is as follows: ∂N ∂t = rN(1 − N K ) − mNP 1 + bN + aN + A sin (ωt) + d 1 ∇ 2 N, ∂P ∂t = P(−q + cmN 1 + bN + aN ) + η (r, t) + d 2 ∇ 2 P, (5) where A sin(ωt) denotes the periodic forcing with amplitude A and angular frequency ω. The colored noise term η(r, t) (r = (x, y)) 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: ∂η (r, t) ∂t = − 1 τ η (r, t) + 1 τ ξ (r, t) , (6) where ξ(r, t) is a Gaussian white noise or the so-called Markovian random telegraph process in both space and time with zero mean and correlation: ⟨ξ (r, t)⟩ = 0, ⟨ξ (r, t) ξ (r󸀠, t󸀠)⟩ = 2εδ (r − r󸀠) δ (t − t󸀠) , (7) where ⟨⋅⟩ denotes averaging with respect to the noise ξ(r, t) and δ the Dirac delta-function and δ(r − r) is the spatial correlation function of the Gaussian white noise ξ(r, t). Integrating (6) with respect to time t, we get η (r, t) = η (r, 0) e−t/τ + 1 τ e −t/τ ∫ t 0 e s/τ ξ (r, s) ds. (8) The mean value of the colored noise is ⟨η (r, t)⟩ = ⟨η (r, 0)⟩ e−t/τ + 1 τ e −t/τ


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, predatorprey model is a mathematical method to approximate some part of our real world.And the dynamic behavior of predatorprey 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.
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,.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  1 and  2 are the diffusion coefficients, respectively, and ∇ 2 = / 2 +/ 2 is the usual Laplacian operator in twodimensional space; other parameters are the same definition as those in model (3).
It is easy to know that when a spatially extended predatorprey 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,[58][59][60][61][62][63].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 twodimensional 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  sin() denotes the periodic forcing with amplitude  and angular frequency .The colored noise term (r, ) (r = (,)) 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 (r, ) 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 (r, ) and  the Dirac delta-function and (r − r  ) is the spatial correlation function of the Gaussian white noise (r, ).
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 (r, ) 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  = 1.
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.
Hopf bifurcation is an instability induced by the transformation of the stability of a focus.Mathematically speaking, Hopf bifurcation occurs when Im() ̸ = 0 and Re() = 0, at  = 0; Im() is the imaginary part, Re() 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   = Im() = √det( 0 ), and the corresponding wavelength   satisfies   = 2/  = 2/√det( 0 ).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 "Im() = 0 and Re() = 0, at  =   ̸ = 0," and the wavenumber   satisfies  2  = √det( 0 )/ 1  2 .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   = 2/  .And the Turing bifurcation surface is given by where trace and the critical value of Turing bifurcation can be obtained from (21) as follows: where Linear stability analysis yields the bifurcation diagram with  = 1,  = 0.125,  = 1,  = 0.7,  = 0.625,  = 0.18, and  2 = 0.2 as shown in Figure 1(a).In this case, parameters (, , , , , , ) ∈  1 , 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  = 2.8 >   = 2.279, which is located in domain II; one can see that when  = 0, Re(()) > 0 and Im(()) ̸ = 0.

Spatiotemporal Dynamics of the Models
In this section, we perform extensive numerical simulations of the spatially extended models ( 4) and ( 5) in twodimensional space, and the qualitative results are shown here.All our numerical simulations employ the zero-flux Neumann boundary conditions with a system size of 200×200 space units.The parameters are  = 1,  = 0.125,  = 1,  = 0.7,  = 0.625,  = 0.18,  1 = 0.02,  2 = 0.2, and  = 2.8 or  = 4.0, which satisfy (, , , , , , ) ∈  1 .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 (±5 × 10 −4 ), 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 Δ = 1/24 and space stepsize (lattice constant) of Δ = Δ = 1.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 .(4). Figure 2 shows the evolution of the spatial patterns of prey  at  = 0, 100, 300, 500, 1000, and 2000, with random small perturbation of the equilibrium ( * ,  * ) = (0.748, 2.132) of model ( 4) with  = 2.8, located in domain II, more than the Hopf bifurcation threshold   = 2.279 and less than the Turing bifurcation threshold   = 3.499.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.

Pattern Formation of Model
When  = 4.0 >   >   , in this case, parameters in domain III (Figure 1(a)) and both Hopf and Turing instabilities occur.The nontrivial stationary state is ( * ,  * ) = (0.748, 2.365).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  = 2.8, the wavelength  = 3.100 while, at  = 4.0,  = 3.021.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  = 4.0, 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  = 2.8 (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.,  2 >  1 ), 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  Abstract and Applied Analysis  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.(5).Now, we turn our focus on the effect of noise on the predator  of stochastic model (5).In this case,  = 0; 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),  = 0.0001; the second row, (b),  = 0.01; the third row, (c),  = 0.05; and the last row of Figure 4, (d),  = 0.1.And the first column of Figure 4, marked as (i), shows the snapshots of spatiotemporal pattern of model ( 5) at  = 2000 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.(5).In the previous subsection, we have shown the effect of noise on the predator  of model ( 5).An interesting question is whether such noisesustained oscillations can be entrained by a weak external forcing, in this case,  = 0.This is investigated here.When model ( 5) is noise free, there is a phenomenon of frequency locking or resonant response [56,[58][59][60][61].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  in = 2/.Above the threshold, model (5) may produce oscillations about period  out with respect to external period  in , which is called frequency locking or resonant response.That is, the model produces one spike within each of the  =  out / in periods of the external force, called  : 1 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  out is defined as follows:   is the time interval between the th spike and ( + 1)th spike. spikes are taken into account and the average value of them is  out = ∑ −1 =1   /( − 1).As an example, with the amplitude  = 0.001, Figure 5 shows 5 : 1 resonant response with  = 0.2 (a) and  = 0.02 (c), respectively.And Figures 5(b) and 5(d) are the phase portraits corresponding to (a) and (c).We can see that when  = 0.2, there exists a periodic orbit, while,  = 0.02,   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.(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 .In Figure 7, we have shown a typical pattern formation process in the 5 : 1 frequency locking regime with  = 0.001 and  = 0.2.From  = 1870 (a) to  = 1920 (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.

Conclusions and Remarks
In this paper, we present a spatial Holling-type IV predatorprey 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 4).The solid curve is time series of prey ; the dash curve is the corresponding external periodic forcing.Other parameters are the same as those in Figure 3.
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.Figure 7: Typical pattern formation of the forced noisy prey in the 5 : 1 locking region at  = 0.001 and  = 0.0001 corresponding to Figure 6(a).The lower panel shows the time series of the prey  (the solid curve) and the corresponding external periodic forcing (the dash curve) corresponding to the snapshots of the patterns.The grey scale, from black to white, is in [0, 3.5] for all the snapshots.
Figures 6(b) and 6(d) illustrate the spiral pattern at  = 2000 corresponding to (a) and (c), respectively.In contrast, we change one of the parameters of Figure 6(c)  = 0.001 to  = 0.01 (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.

Figure 5 :
Figure 5: External periodic forcing induced frequency locking of model (4).The solid curve is time series of prey ; the dash curve is the corresponding external periodic forcing.Other parameters are the same as those in Figure3.

Figure 6 :
Figure 6: Dynamics of model (5) with both noise and periodic forcing.(b, d, and f) are snapshots at  = 2000 corresponding to the left hand side resonant response.The other parameters are the same as those in Figure 3.