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

Authors: Lei Zhang; Zhibin Li
Journal: Abstract and Applied Analysis
 (2014)

---

## 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.


---

## Body

## 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]
(1)
d
N
d
t
=
N
f
(
N
)
-
m
P
g
(
N
,
P
)
,
d
P
d
t
=
P
[
c
m
g
(
N
,
P
)
-
d
]
,
where N and P stand for prey and predator quantity, respectively, f
(
N
) is the per capita rate of increase of the prey in absence of predation, d is the food-independent death rate of predator, g
(
N
,
P
) 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, m
g
(
N
,
P
) is the amount of prey consumed per predator per unit time, and c
m
g
(
N
,
P
) is the predator production per capita with predation.In population dynamics, a functional responseg
(
N
,
P
) 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 [6–9]. In the history of population ecology, both ecologists and mathematicians have a great interest in the Holling-type predator-prey models [3, 8, 10–21], 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 g
(
N
,
P
) in (1) is a function only related to prey N. The classical expression of Holling-type II functional response is g
(
N
,
P
)
=
m
N
/
(
1
+
b
N
), and g
(
N
,
P
)
=
m
N
2
/
(
1
+
a
N
2
) is called Holling-type III. The Holling-type IV functional response is written as follows:
(2)
g
(
N
,
P
)
=
m
N
1
+
b
N
+
a
N
2
.
Function (2) is called Monod-Haldane-type functional response too [25]. In addition, when b
=
0, a simplified form g
(
N
,
P
)
=
m
N
/
(
1
+
a
N
2
) 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
(3)
d
N
d
t
=
r
N
(
1
-
N
K
)
-
m
N
P
1
+
b
N
+
a
N
2
,
d
P
d
t
=
P
(
-
q
+
c
m
N
1
+
b
N
+
a
N
2
)
,
where r
>
0 stands for maximum per capita growth rate of the prey, m
>
0 is the capture rate, c
>
0 is the conversion rate of prey captured by predator, q
>
0 is the food-independent death rate of predator, K
>
0 is the carrying capacity, and a
>
0 is the so-called half-saturation constant; b
>
-
2
a such that the denominator of above system does not vanish for nonnegative N.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, 31–56]. 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(4)
∂
N
∂
t
=
r
N
(
1
-
N
K
)
-
m
N
P
1
+
b
N
+
a
N
2
+
d
1
∇
2
N
,
∂
P
∂
t
=
P
(
-
q
+
c
m
N
1
+
b
N
+
a
N
2
)
+
d
2
∇
2
P
,
where d
1 and d
2 are the diffusion coefficients, respectively, and ∇
2
=
∂
/
∂
x
2
+
∂
/
∂
y
2 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, 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 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:
(5)
∂
N
∂
t
=
r
N
(
1
-
N
K
)
-
m
N
P
1
+
b
N
+
a
N
2
+
A
sin
(
ω
t
)
+
d
1
∇
2
N
,
∂
P
∂
t
=
P
(
-
q
+
c
m
N
1
+
b
N
+
a
N
2
)
+
η
(
r
,
t
)
+
d
2
∇
2
P
,
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:
(6)
∂
η
(
r
,
t
)
∂
t
=
-
1
τ
η
(
r
,
t
)
+
1
τ
ξ
(
r
,
t
)
,
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:
(7)
〈
ξ
(
r
,
t
)
〉
=
0
,
〈
ξ
(
r
,
t
)
ξ
(
r
′
,
t
′
)
〉
=
2
ɛ
δ
(
r
-
r
′
)
δ
(
t
-
t
′
)
,
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
(8)
η
(
r
,
t
)
=
η
(
r
,
0
)
e
-
t
/
τ
+
1
τ
e
-
t
/
τ
∫
0
t
‍
e
s
/
τ
ξ
(
r
,
s
)
d
s
.
The mean value of the colored noise is
(9)
〈
η
(
r
,
t
)
〉
=
〈
η
(
r
,
0
)
〉
e
-
t
/
τ
+
1
τ
e
-
t
/
τ
×
∫
0
t
e
s
/
τ
〈
ξ
(
r
,
s
)
〉
d
s
=
〈
η
(
r
,
0
)
〉
e
-
t
/
τ
,
and the correlation function of the colored noise is given by
(10)
〈
η
(
r
,
t
)
η
(
r
′
,
t
′
)
〉
=
〈
η
(
r
,
0
)
〉
〈
η
(
r
′
,
0
)
〉
e
-
(
t
+
t
′
)
/
τ
+
1
τ
2
e
-
(
t
+
t
′
)
/
τ
×
∫
0
t
‍
∫
0
t
′
‍
e
(
s
+
s
′
)
/
τ
〈
ξ
(
r
,
s
)
ξ
(
r
′
,
s
′
)
〉
d
s
d
s
′
=
〈
η
(
r
,
0
)
〉
〈
η
(
r
′
,
0
)
〉
e
-
(
t
+
t
′
)
/
τ
+
ɛ
τ
2
e
-
(
t
+
t
′
)
/
τ
δ
(
r
-
r
′
)
×
∫
0
t
‍
∫
0
t
′
‍
e
(
s
+
s
′
)
/
τ
δ
(
t
-
t
′
)
d
s
d
s
′
=
〈
η
(
r
,
0
)
〉
〈
η
(
r
′
,
0
)
〉
e
-
(
t
+
t
′
)
/
τ
+
ɛ
τ
δ
(
r
-
r
′
)
×
(
e
-
(
t
+
t
′
)
/
τ
-
2
e
-
t
/
τ
+
e
-
(
t
-
t
′
)
/
τ
)
.
Let t
→
+
∞; then
(11)
〈
η
(
r
,
t
)
η
(
r
′
,
t
′
)
〉
⟶
ɛ
τ
e
-
(
t
-
t
′
)
/
τ
δ
(
r
-
r
′
)
.The colored noiseη
(
r
,
t
) 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
(12)
〈
η
(
r
,
t
)
η
(
r
′
,
t
′
)
〉
=
ɛ
τ
e
-
|
t
-
t
′
|
/
τ
δ
(
r
-
r
′
)
,
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.


## 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: (
0,0
) (total extinct) is a saddle; (
K
,
0
) (extinct of the predator or prey-only) is a attracting node if q
>
c
m
K
/
(
1
+
K
b
+
a
K
2
), a saddle if q
<
c
m
K
/
(
1
+
K
b
+
a
K
2
), or a saddle-node if q
=
c
m
K
/
(
1
+
K
b
+
a
K
2
). When
(13)
(
a
,
b
,
c
,
m
,
q
,
r
,
K
)
∈
E
1
,
here
,
E
1
=
{
-
m
2
c
2
+
2
q
b
m
c
+
q
a
K
m
c
-
q
2
a
K
b
+
2
q
2
a
-
q
2
b
2
-
m
c
+
q
b
+
q
a
K
(
a
,
b
,
c
,
m
,
q
,
r
,
K
)
∣
m
c
>
q
b
,
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
q
2
a
<
(
m
c
-
q
b
)
2
,
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
m
c
q
a
(
m
c
-
q
b
)
2
-
4
q
2
a
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
>
(
-
m
2
c
2
+
2
q
b
m
c
+
q
a
K
m
c
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
iiiiii
-
q
2
a
K
b
+
2
q
2
a
-
q
2
b
2
)
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
×
(
-
m
c
+
q
b
+
q
a
K
)
-
1
>
0
,
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
b
a
-
m
c
q
a
+
K
<
0
-
m
2
c
2
+
2
q
b
m
c
+
q
a
K
m
c
-
q
2
a
K
b
+
2
q
2
a
-
q
2
b
2
-
m
c
+
q
b
+
q
a
K
}
,
there exists unique stationary coexistence state (
N
1
*
,
P
1
*
), where
(14)
N
1
*
=
1
2
-
q
b
+
m
c
-
A
q
a
,
P
1
*
=
c
r
(
(
-
m
c
+
b
q
+
q
a
K
)
N
1
*
+
q
)
a
q
2
K
.
On the other hand, when
(15)
(
a
,
b
,
c
,
m
,
q
,
r
,
K
)
∈
E
2
,
here
,
E
2
=
{
-
m
2
c
2
+
2
q
b
m
c
+
q
a
K
m
c
-
q
2
a
K
b
+
2
q
2
a
-
q
2
b
2
-
m
c
+
q
b
+
q
a
K
(
a
,
b
,
c
,
m
,
q
,
r
,
K
)
∣
m
c
>
q
b
,
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
q
2
a
<
(
m
c
-
q
b
)
2
,
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
m
c
q
a
(
m
c
-
q
b
)
2
-
4
q
2
a
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
>
-
(
-
m
2
c
2
+
2
q
b
m
c
+
q
a
K
m
c
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
iiiiii
iii
-
q
2
a
K
b
+
2
q
2
a
-
q
2
b
2
)
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
×
(
-
m
c
+
q
b
+
q
a
K
)
-
1
>
0
,
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
b
a
-
m
c
q
a
+
K
>
0
-
m
2
c
2
+
2
q
b
m
c
+
q
a
K
m
c
-
q
2
a
K
b
+
2
q
2
a
-
q
2
b
2
-
m
c
+
q
b
+
q
a
K
}
,
there exists another unique stationary coexistence state (
N
2
*
,
P
2
*
) implying
(16)
N
2
*
=
1
2
-
q
b
+
m
c
+
A
q
a
,
P
2
*
=
c
r
(
(
-
m
c
+
b
q
+
q
a
K
)
N
2
*
+
q
)
a
q
2
K
.It is worth mentioning that equilibria(
N
1
*
,
P
1
*
) and (
N
2
*
,
P
2
*
) cannot coexist. In this paper, we mainly focus on the dynamics of (
N
1
*
,
P
1
*
) and rewrite it as (
N
*
,
P
*
). The dynamics behavior of (
N
2
*
,
P
2
*
) is similar to that of (
N
1
*
,
P
1
*
).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:
(17)
N
(
r
,
t
)
~
N
*
e
λ
t
e
i
k
→
·
r
,
P
(
r
,
t
)
~
P
*
e
λ
t
e
i
k
→
·
r
,
iiiiiiiiiiiiiiiiiiiiii
iii
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 whenIm
⁡
(
λ
)
≠
0and Re
(
λ
)
=
0, at k
=
0; Im
⁡
(
λ
) is the imaginary part, Re
(
λ
) is the real part, and k is the wave number. So we get the Hopf bifurcation surface:
(18)
H
=
{
(
a
,
b
,
c
,
m
,
q
,
r
,
K
)
∣
det
⁡
(
J
0
)
>
0
,
i
i
i
i
i
i
i
trace
(
J
0
)
=
0
}
,
where
(19)
det
⁡
(
J
0
)
=
-
(
r
-
2
r
N
*
K
)
q
+
m
q
P
*
+
c
m
(
r
-
2
(
r
N
*
/
K
)
)
N
*
(
1
+
b
N
*
+
a
N
*
2
)
-
m
q
N
*
P
*
(
b
+
2
a
N
*
)
(
1
+
b
N
*
+
a
N
*
2
)
2
,
trace
(
J
0
)
=
r
-
2
r
N
*
K
-
q
+
m
(
-
P
*
+
c
N
*
+
a
N
*
2
P
*
+
b
c
N
*
2
+
c
N
*
3
a
)
(
1
+
b
N
*
+
a
N
*
2
)
2
,
the frequency of periodic oscillations in time ω
H satisfies ω
H
=
Im
⁡
(
λ
)
=
det
⁡
(
J
0
), and the corresponding wavelength λ
H satisfies λ
H
=
2
π
/
ω
H
=
2
π
/
det
⁡
(
J
0
). In particular, we take K as the bifurcation parameter and can get the critical value of Hopf bifurcation from (18):
(20)
K
H
=
(
-
(
a
q
2
(
5
m
c
-
3
q
b
)
-
(
3
m
c
-
q
b
)
(
m
c
-
q
b
)
2
)
i
i
i
i
i
×
(
m
c
-
q
b
)
2
-
4
q
2
a
-
4
q
4
a
2
+
q
2
(
m
c
-
q
b
)
i
i
i
i
i
×
(
11
m
c
-
5
q
b
)
a
-
(
3
m
c
-
q
b
)
(
m
c
-
q
b
)
3
)
i
i
i
×
(
(
-
a
q
(
(
(
2
m
c
-
q
b
)
(
m
c
-
q
b
)
-
2
q
2
a
)
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
×
(
m
c
-
q
b
)
2
-
4
q
2
a
-
2
a
q
2
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
×
(
3
m
c
-
2
q
b
)
+
(
2
m
c
-
q
b
)
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
×
(
m
c
-
q
b
)
2
)
)
)
-
1
.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 k
=
k
T
≠
0,” and the wavenumber k
T satisfies k
T
2
=
det
⁡
(
J
0
)
/
d
1
d
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 λ
T
=
2
π
/
k
T. And the Turing bifurcation surface is given by
(21)
T
=
{
(
a
,
b
,
c
,
m
,
q
,
r
,
d
1
,
d
2
,
K
)
∣
det
⁡
(
J
k
)
=
0
,
i
i
i
trace
(
J
k
)
=
0
}
,
where
(22)
det
⁡
(
J
k
)
=
-
(
r
-
2
r
N
*
K
-
d
1
k
2
)
(
q
+
d
2
k
2
)
×
(
(
q
+
d
2
k
2
)
m
P
*
mmmm
+
(
r
-
2
(
r
N
*
K
)
-
d
1
k
2
)
c
m
N
*
(
q
+
d
2
k
2
)
)
×
(
1
+
b
N
*
+
a
N
*
2
)
-
1
-
m
(
b
+
2
a
N
*
)
(
q
+
d
2
k
2
)
N
*
P
*
(
1
+
b
N
*
+
a
N
*
2
)
2
,
(23)
trace
(
J
k
)
=
r
-
2
r
N
*
K
-
q
-
(
d
1
+
d
2
)
k
2
+
m
(
-
P
*
+
c
N
*
+
a
N
*
2
P
*
+
b
c
N
*
2
+
c
N
*
3
a
)
(
1
+
b
N
*
+
a
N
*
2
)
2
,
and the critical value of Turing bifurcation can be obtained from (21) as follows:
(24)
K
T
=
F
1
F
2
,
where
(25)
F
1
=
r
(
(
4
q
2
a
(
2
m
c
-
q
b
)
i
i
i
i
i
i
i
i
-
(
3
m
c
-
q
b
)
(
m
c
-
q
b
)
2
)
i
i
i
i
i
i
i
i
×
(
m
c
-
q
b
)
2
-
4
q
2
a
i
i
i
i
i
i
i
i
+
(
(
m
c
-
q
b
)
2
-
4
q
2
a
)
i
i
i
i
i
i
·
(
(
3
m
c
-
q
b
)
(
m
c
-
q
b
)
-
2
q
2
a
)
)
×
(
(
m
c
-
q
b
)
2
-
4
q
2
a
(
3
m
c
-
q
b
)
(
m
c
-
q
b
)
-
4
q
2
a
-
(
3
m
c
+
q
b
)
s
s
s
s
s
s
×
(
m
c
-
q
b
)
2
-
4
q
2
a
)
d
2
,
F
2
=
q
a
(
(
m
c
-
q
b
)
2
-
4
q
2
a
2
m
c
(
(
m
c
-
q
b
)
(
m
c
-
q
b
)
2
-
4
q
2
a
iiiiiiiiiiiiiiiiiii
+
4
q
2
a
-
(
m
c
-
q
b
)
2
(
m
c
-
q
b
)
2
-
4
q
2
a
)
B
iiiiiiiii
+
(
(
m
c
-
q
b
)
2
-
4
q
2
a
(
3
m
c
-
q
b
)
d
2
iiiiiiiiiiiiiii
·
(
2
m
2
c
2
-
3
q
b
m
c
+
q
2
b
2
-
4
q
2
a
)
r
iiiiiiiiiiiiiii
+
2
q
d
1
m
c
(
m
c
-
q
b
)
2
(
m
c
-
q
b
)
2
-
4
q
2
a
)
i
i
i
i
i
i
i
i
×
(
(
m
c
-
q
b
)
2
-
4
q
2
a
)
i
i
i
i
i
i
i
i
+
(
-
2
d
2
(
(
3
m
c
-
q
b
)
(
2
m
c
-
q
b
)
(
m
c
-
q
b
)
2
iiiiiiiiiiiiiiiiiiiiiiiiii
-
2
q
2
a
(
-
4
q
2
a
+
3
q
2
b
2
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
ss
iii
-
12
q
b
m
c
+
11
m
2
c
2
)
)
r
i
i
i
i
i
i
i
i
i
i
i
i
i
-
4
q
c
m
(
m
c
-
q
b
)
d
1
(
(
m
c
-
q
b
)
2
-
4
q
2
a
)
)
iiiiiiiii
×
(
m
c
-
q
b
)
2
-
4
q
2
a
iiiiiiiii
+
(
q
2
b
2
-
2
q
b
m
c
+
m
2
c
2
-
4
q
2
a
)
iiiiiiiii
×
(
(
m
c
-
q
b
)
2
-
4
q
2
a
(
2
m
c
-
q
b
)
d
2
(
3
m
2
c
2
-
4
q
b
m
c
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
+
q
2
b
2
-
4
q
2
a
)
r
iiiiiiiiiiiii
+
2
m
c
q
d
1
(
(
m
c
-
q
b
)
2
-
4
q
2
a
)
(
2
m
c
-
q
b
)
d
2
(
3
m
2
c
2
-
4
q
b
m
c
(
m
c
-
q
b
)
2
-
4
q
2
a
)
(
m
c
-
q
b
)
2
-
4
q
2
a
)
,
B
=
(
(
m
c
-
q
b
)
2
-
4
q
2
a
(
m
c
-
q
b
)
2
-
4
q
2
a
-
2
d
1
q
(
(
(
m
c
-
q
b
)
2
(
d
1
b
q
2
-
q
d
1
m
c
-
q
r
b
d
2
+
2
r
m
c
d
2
)
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
×
(
4
q
2
a
-
(
m
c
-
q
b
)
2
)
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
-
r
m
c
d
2
(
m
c
-
q
b
)
2
)
iiiiiiiiiiiiiiiiiii
·
(
m
c
-
q
b
)
2
-
4
q
2
a
iiiiiiiiiiiiiiiiiii
+
(
8
q
4
d
2
r
-
8
q
5
d
1
)
a
2
iiiiiiiiiiiiiiiiiii
+
4
(
m
c
-
q
b
)
a
q
2
r
m
c
d
2
iiiiiiiiiiiiiiiiiii
+
(
3
r
m
c
d
2
-
q
d
1
m
c
+
d
1
b
q
2
-
q
r
b
d
2
)
iiiiiiiiiiiiiiiiiii
×
(
(
m
c
-
q
b
)
3
-
6
(
m
c
-
q
b
)
a
q
2
)
)
(
m
c
-
q
b
)
2
-
4
q
2
a
)
1
/
2
.Linear stability analysis yields the bifurcation diagram withr
=
1, a
=
0.125, b
=
1, c
=
0.7, m
=
0.625, q
=
0.18, and d
2
=
0.2 as shown in Figure 1(a). In this case, parameters (
a
,
b
,
c
,
m
,
q
,
r
,
K
)
∈
E
1, and (
N
*
,
P
*
) 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 K
=
2.8
>
K
H
=
2.279, which is located in domain II; one can see that when k
=
0, Re
(
λ
(
k
)
)
>
0 and Im
⁡
(
λ
(
k
)
)
≠
0. Figure 1(c) displays the case of the critical value of Turing bifurcation K
=
K
T
=
3.499; in this case, Re
(
λ
(
k
)
)
=
0 and Im
⁡
(
λ
(
k
)
)
=
0 at k
=
k
T
=
2.080. When K
=
4.0, parameters are located in domain III; Figure 1(d) indicates that, at k
=
0, Re
(
λ
(
k
)
)
>
0 and Im
⁡
(
λ
(
k
)
)
≠
0.(a)K
-
d
1 Bifurcation diagram for model (4) with r
=
1, a
=
0.125, b
=
1, c
=
0.7, m
=
0.625, q
=
0.18, and d
2
=
0.2. Hopf and Turing bifurcation lines separate the parameter space into three domains. The other parameters in (b)–(d) are d
1
=
0.02; the bifurcation parameter K equals (b) 2.8
>
K
H
=
2.279; (c) 3.499
=
K
T; (d) 4.0
>
K
T
>
K
H. The real parts Re
(
λ
) and the imaginary parts Im
⁡
(
λ
) are shown by solid curves and dashed curves, respectively.
(a)
(b)
(c)
(d)


## 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 200
×
200 space units. The parameters are r
=
1, a
=
0.125, b
=
1, c
=
0.7, m
=
0.625, q
=
0.18, d
1
=
0.02, d
2
=
0.2, and K
=
2.8 or K
=
4.0, which satisfy (
a
,
b
,
c
,
m
,
q
,
r
,
K
)
∈
E
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 (
N
*
,
P
*
) with small random perturbation superimposed; in each, the initial condition is always a small amplitude random perturbation (
±
5
×
1
0
-
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 Δ
t
=
1
/
24 and space stepsize (lattice constant) of Δ
x
=
Δ
y
=
1. When the system reached a periodic oscillatory state, we took a snapshot with white corresponding to the high value of prey N 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 preyN.


### 3.1. Pattern Formation of Model (4)

Figure2 shows the evolution of the spatial patterns of prey N at t
=
0, 100, 300, 500, 1000, and 2000, with random small perturbation of the equilibrium (
N
*
,
P
*
)
=
(
0.748,2.132
) of model (4) with 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.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 N and periodic oscillating in time finally; (h) exhibits the fact that a limit cycle arises, which is caused by the Hopf bifurcation.Grey-scaled snapshots of spatiotemporal pattern of the preyN of model (4) with K
=
2.8. (a) t
=
0, (b) t
=
100, (c) t
=
300, (d) t
=
500, (e) t
=
1000, and (f) t
=
2000. The lower panels show the corresponding (g) time-series plots and (h) phase portraits.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)WhenK
=
4.0
>
K
T
>
K
H, in this case, parameters in domain III (Figure 1(a)) and both Hopf and Turing instabilities occur. The nontrivial stationary state is (
N
*
,
P
*
)
=
(
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.Grey-scaled snapshots of spatiotemporal pattern of the preyN of model (4) with K
=
4.0. (a) t
=
0, (b) t
=
100, (c) t
=
300, (d) t
=
500, (e) t
=
1000, and (f) t
=
2000. The lower panels show the corresponding (g) time-series plots and (h) phase portraits.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Comparing this situation (Figure3) 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 K
=
2.8, the wavelength λ
=
3.100 while, at K
=
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 K
=
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 K
=
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 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 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 N) should have a positive effect on the creation rate of both activator (prey N) and inhibitor (prey P) species. Second, an increment in inhibitor species should 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-prey model (4), with Holling-type IV functional response and predators diffusing faster than prey (i.e., d
2
>
d
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 N and the predator P 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 P on the inside and prey N 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 predatorP of stochastic model (5). In this case, A
=
0; that is, the periodic forcing is not present.Figure4 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 t
=
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 N 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.Figure 4
Dynamics of model (5), for the following noise intensity. (a) ɛ
=
0.0001; (b) ɛ
=
0.01; (c) ɛ
=
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.


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

In the previous subsection, we have shown the effect of noise on the predatorP of model (5). An interesting question is whether such noise-sustained 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–61]. That is, without noise, the spatially homogeneous oscillation does not respond to the external periodic forcing when the amplitude A is below a threshold whose value depends on the external period T
in
=
2
π
/
ω. Above the threshold, model (5) may produce oscillations about period T
out with respect to external period T
in, which is called frequency locking or resonant response. That is, the model produces one spike within each of the M
=
T
out
/
T
in periods of the external force, called M
:
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 T
out is defined as follows: T
i is the time interval between the ith spike and (
i
+
1
)th spike. m spikes are taken into account and the average value of them is T
out
=
∑
i
=
1
m
-
1
T
i
/
(
m
-
1
).As an example, with the amplitudeA
=
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.External periodic forcing induced frequency locking of model (4). The solid curve is time series of prey U; the dash curve is the corresponding external periodic forcing. Other parameters are the same as those in Figure 3.
(a)
(b)
(c)
(d)


### 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 A and angular frequency ω. Figures 6(b) and 6(d) illustrate the spiral pattern at t
=
2000 corresponding to (a) and (c), respectively. In contrast, we change one of the parameters of Figure 6(c)  A
=
0.001 to A
=
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 A 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.Dynamics of model (5) with both noise and periodic forcing. (b, d, and f) are snapshots at t
=
2000 corresponding to the left hand side resonant response. The other parameters are the same as those in Figure 3.
(a)
(b)
(c)
(d)
(e)
(f)In Figure7, we have shown a typical pattern formation process in the 5 : 1 frequency locking regime with A
=
0.001 and ω
=
0.2
π. From t
=
1870 (a) to t
=
1920 (f), the pattern formation of prey N 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 N does not exactly repeat that due to a small fluctuation of the phase difference.Typical pattern formation of the forced noisy prey in the 5 : 1 locking region atA
=
0.001 and ɛ
=
0.0001 corresponding to Figure 6(a). The lower panel shows the time series of the prey N (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.
(a)
(b)
(c)
(d)
(e)
(f)
(g)


## 3.1. Pattern Formation of Model (4)

Figure2 shows the evolution of the spatial patterns of prey N at t
=
0, 100, 300, 500, 1000, and 2000, with random small perturbation of the equilibrium (
N
*
,
P
*
)
=
(
0.748,2.132
) of model (4) with 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.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 N and periodic oscillating in time finally; (h) exhibits the fact that a limit cycle arises, which is caused by the Hopf bifurcation.Grey-scaled snapshots of spatiotemporal pattern of the preyN of model (4) with K
=
2.8. (a) t
=
0, (b) t
=
100, (c) t
=
300, (d) t
=
500, (e) t
=
1000, and (f) t
=
2000. The lower panels show the corresponding (g) time-series plots and (h) phase portraits.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)WhenK
=
4.0
>
K
T
>
K
H, in this case, parameters in domain III (Figure 1(a)) and both Hopf and Turing instabilities occur. The nontrivial stationary state is (
N
*
,
P
*
)
=
(
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.Grey-scaled snapshots of spatiotemporal pattern of the preyN of model (4) with K
=
4.0. (a) t
=
0, (b) t
=
100, (c) t
=
300, (d) t
=
500, (e) t
=
1000, and (f) t
=
2000. The lower panels show the corresponding (g) time-series plots and (h) phase portraits.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Comparing this situation (Figure3) 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 K
=
2.8, the wavelength λ
=
3.100 while, at K
=
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 K
=
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 K
=
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 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 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 N) should have a positive effect on the creation rate of both activator (prey N) and inhibitor (prey P) species. Second, an increment in inhibitor species should 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-prey model (4), with Holling-type IV functional response and predators diffusing faster than prey (i.e., d
2
>
d
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 N and the predator P 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 P on the inside and prey N 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 predatorP of stochastic model (5). In this case, A
=
0; that is, the periodic forcing is not present.Figure4 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 t
=
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 N 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.Figure 4
Dynamics of model (5), for the following noise intensity. (a) ɛ
=
0.0001; (b) ɛ
=
0.01; (c) ɛ
=
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.


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

In the previous subsection, we have shown the effect of noise on the predatorP of model (5). An interesting question is whether such noise-sustained 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–61]. That is, without noise, the spatially homogeneous oscillation does not respond to the external periodic forcing when the amplitude A is below a threshold whose value depends on the external period T
in
=
2
π
/
ω. Above the threshold, model (5) may produce oscillations about period T
out with respect to external period T
in, which is called frequency locking or resonant response. That is, the model produces one spike within each of the M
=
T
out
/
T
in periods of the external force, called M
:
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 T
out is defined as follows: T
i is the time interval between the ith spike and (
i
+
1
)th spike. m spikes are taken into account and the average value of them is T
out
=
∑
i
=
1
m
-
1
T
i
/
(
m
-
1
).As an example, with the amplitudeA
=
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.External periodic forcing induced frequency locking of model (4). The solid curve is time series of prey U; the dash curve is the corresponding external periodic forcing. Other parameters are the same as those in Figure 3.
(a)
(b)
(c)
(d)


## 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 A and angular frequency ω. Figures 6(b) and 6(d) illustrate the spiral pattern at t
=
2000 corresponding to (a) and (c), respectively. In contrast, we change one of the parameters of Figure 6(c)  A
=
0.001 to A
=
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 A 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.Dynamics of model (5) with both noise and periodic forcing. (b, d, and f) are snapshots at t
=
2000 corresponding to the left hand side resonant response. The other parameters are the same as those in Figure 3.
(a)
(b)
(c)
(d)
(e)
(f)In Figure7, we have shown a typical pattern formation process in the 5 : 1 frequency locking regime with A
=
0.001 and ω
=
0.2
π. From t
=
1870 (a) to t
=
1920 (f), the pattern formation of prey N 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 N does not exactly repeat that due to a small fluctuation of the phase difference.Typical pattern formation of the forced noisy prey in the 5 : 1 locking region atA
=
0.001 and ɛ
=
0.0001 corresponding to Figure 6(a). The lower panel shows the time series of the prey N (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.
(a)
(b)
(c)
(d)
(e)
(f)
(g)


## 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 parameterK, 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 N, while it begins from the predator P if K 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.



---

*Source: 675378-2014-06-01.xml*