Emergence of spatiotemporal chaos driven by far-field breakup of spiral
  waves in the plankton ecological systems
Emergence of spatiotemporal chaos driven by far-field breakup of spiral waves in the
plankton ecological systems
Quan-Xing Liu,1 Gui-Quan Sun,1 Bai-Lian Li,2 and Zhen Jin1, ∗
Department of Mathematics, North University of China,
Taiyuan, Shan’xi 030051, People’s Republic of China
Ecological Complexity and Modeling Laboratory, Department of Botany and Plant Sciences,
University of California, Riverside, CA 92521-0124, USA
(Dated: October 25, 2018)
Alexander B. Medvinsky et al [A. B. Medvinsky, I. A. Tikhonova, R. R. Aliev, B.-L. Li, Z.-S. Lin,
and H. Malchow, Phys. Rev. E 64, 021915 (2001)] and Marcus R. Garvie et al [M. R. Garvie and C.
Trenchea, SIAM J. Control. Optim. 46, 775-791 (2007)] shown that the minimal spatially extended
reaction-diffusion model of phytoplankton-zooplankton can exhibit both regular, chaotic behavior,
and spatiotemporal patterns in a patchy environment. Based on that, the spatial plankton model
is furtherly investigated by means of computer simulations and theoretical analysis in the present
paper when its parameters would be expected in the case of mixed Turing-Hopf bifurcation region.
Our results show that the spiral waves exist in that region and the spatiotemporal chaos emerge,
which arise from the far-field breakup of the spiral waves over large ranges of diffusion coefficients
of phytoplankton and zooplankton. Moreover, the spatiotemporal chaos arising from the far-field
breakup of spiral waves does not gradually involve the whole space within that region. Our results
are confirmed by means of computation spectra and nonlinear bifurcation of wave trains. Finally,
we give some explanations about the spatially structured patterns from the community level.
PACS numbers: 87.23.Cc, 82.40.Ck, 82.40.Bj, 92.20.jm
Keywords: Spiral waves; Spatio-temporal pattern; Plankton dynamics; Reaction-diffusion system
I. INTRODUCTION
There is a growing interest in the spatial pattern dy-
namics of ecological systems [1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
11, 12, 13]. However, many mechanisms of the spatio-
temporal variability of natural plankton populations are
not known yet. Pronounced physical patterns like ther-
moclines, upwelling, fronts and eddies often set the frame
for the biological process. Measurements of the underwa-
ter light field are made with state-of-the-art instruments
and used to calculate concentrations of phytoplankton
biomass (as chlorophyll) as well as other forms of organic
matter. Very high diffusion of the marine environment
would prevent the formation of any stable patch spatial
distribution with much longer life-time than the typical
time of biodynamics. Meanwhile, in addition to very
changeable transient spatial patterns, there also exist
other spatial patterns in marine environment, much more
stable spatial structure associated with ocean fronts, spa-
tiotemporal chaos [10, 11, 14], cyclonic rings, and so
called meddies [15]. In fact, it is significant to create
the biological basis for understanding spatial patterns
of plankton [16]. For instance, the impact of space on
the persistence of enriched ecological systems was proved
in laboratory experiments [17]. Recently, it has been
shown both in laboratory experiments [18] and theoreti-
cally [14, 19, 20, 21] that the existence of a spatial struc-
ture makes a predator-prey system less prone to extinc-
∗Corresponding author; Electronic address: jinzhn@263.net
tion. This is due to the temporal variations of the density
of different sub-populations can become asynchronous
and the events of local extinction can be compensated
due to re-colonization from other sites in the space [22].
During a long period of time, all the spiral waves have
been widely observed in diverse physical, chemical, and
biological systems [23, 24, 25, 26]. However, a quite lim-
ited number of documents [11, 12, 27, 28, 29] concern
the spiral wave pattern and its breakup in the ecological
systems.
The investigation of transition from regular patterns
to spatiotemporally chaotic dynamics in spatially ex-
tended systems remains a challenge in nonlinear sci-
ence [14, 23, 30, 31]. In a nonlinear ecology system, the
two most commonly seen patterns are spiral waves and
turbulence (spatio-temporal chaos) for the level of the
community [32]. It has been recently shown that sponta-
neous spatiatemoporal pattern formation is an instrinsic
property of a predator-prey system [11, 14, 33, 34, 35, 36]
and spatiotemporal structures play an important role in
ecological systems. For example, spatially induced speci-
ation prevents the extinction of the predator-prey mod-
els [11, 12, 37]. So far, plankton patchiness has been ob-
served on a wide range of spatial temporal scales [38, 39].
There exist various, often heuristic explanations of the
spatial patterns phenomenon for these systems. It should
be noted that, although conclusive evidence of ecological
chaos is still to be found, there is a growing number of
indications of chaos in real ecosystems [40, 41, 42, 43].
Recently developed models show that spatial self-
structuring in multispecies systems can meet both cri-
teria and provide a rich substrate for community-level
http://arxiv.org/abs/0704.0322v3
mailto:jinzhn@263.net
section and a major transition in evolution. In present
paper, the scenario in the spatially extended plankton
ecological system is observed by means of the numeri-
cal simulation. The system has been demonstrated to
exhibit regular or chaostic, depending on the initial con-
ditions and the parameter values [10, 29]. We find that
the far-field breakup of the spiral wave leads to complex
spatiotemporal chaos (or a turbulentlike state) in the spa-
tially extended plankton model (1). Our results show
that regular spiral wave pattern shifts into spatiotempo-
ral chaos pattern by modulating the diffusion coefficients
of the species.
II. MODEL
In this paper we study the spatially extended nutrient-
phytoplankton-zooplankton-fish reaction-diffusion sys-
tem. Following Scheffer’s minimal approach [44], which
was originally formulated as a system of ordinary diff-
ential equation (ODEs) and later developed models [10,
11, 29, 45, 46], as a further investigation, we study a
two-variable phytoplankton and zooplankton model on
the level of the community to describe pattern formation
with the diffusion. The dimensionless model is written
= rp(1 − p)−
1 + bp
h+ dp∇
2p, (1a)
1 + bp
h−mh− f
n2 + h2
+ dh∇
2h, (1b)
where the parameters are r, a, b, m, n, dp, dh, and
f which refer to work in Refs. [10, 11]. The explana-
tion of model (1) relates to the nutrient-phytoplankton-
zooplankton-fish ecological system [see Refs. [10, 29, 44]
for details]. The local dynamics are given by
g1(p, h) = rp(1− p)−
1 + bp
h, (2a)
g2(p, h) =
1 + bp
h−mh− f
n2 + h2
. (2b)
From the earlier results [45] about non-spatial system
of model (1) by means of numerical bifurcation analysis
show that the bifurcation and bistability can be found in
the system (1) when the parameters are varied within a
realistic range. For the fixed parameters (see the caption
of Fig. 1 and 2), we can see that the f controls the dis-
tance from Hopf bifurcation. For larger f , there exists
only one stable steady state. As f is decreased further,
the homogeneous steady state undergoes a saddle node
bifurcation (SN), that is fSN = 0.658. In this case, a
stable and an unstable steady state become existence.
Moreover, the bistability will emerge when the parame-
ter f lies the interval fSN > f > fc = 0.445 (this value is
more than the Hopf onset, fH = 0.3397). There are three
steady states: with these kinetics A and C are linearly
stable while B is unstable. Outside this interval, the sys-
tem (1) has unique nontrivial equilibrium. Recent stud-
ies [11, 29] shown that the systems (1) can well-develop
the spiral waves in the oscillation regime, but where the
authors only consider the special case, i.e., dp = dh. A
few important issue have not yet been properly addressed
such as the spatial pattern if dp 6= dh.
Here we report the result that emergence of spatiotem-
poral chaos due to breakup in the system under the
dh 6= dp case. We may now use the f and diffusion
ratio, ν = dh/dp, as control parameters to evaluate
the region for the spiral wave. Turing instability in
reaction-diffusion can be recast in terms of matrix sta-
bility [47, 48]. Such with the help of Maple software
assistance algebra computing, we obtain the parameters
space (f, ν) bifurcation diagrams of the spiral waves as
showing Fig. 2, in which two lines are plotted, Hopf line
(solid) and Turing lines (dotted) respectively. In domain
I, located above all three bifurcation lines, the homo-
geneous steady states is the only stable solution of the
system. Domain II are regions of homogeneous oscilla-
tion in two dimensional spaces [49]. In domain III, both
Hopf and Turing instabilities occur, (i.e., mixed Turing-
Hopf modes arise), in which the system generally pro-
duces the phase waves. Our results show that the system
has spiral wave in this regions. One can see that a Hopf
bifurcation can occur at the steady when the parameter
f passes through a critical values fH while the diffusion
coefficients dp = dh = 0 and the bifurcation periodic so-
lutions are stable. From our analysis (see Fig. 2), one
could also see that the diffusion can induce Turing type
instability for the spatial homogeneous stable periodic
solutions and the spatially extended model (1) exhibit
spatio-temporal chaos patterns. These spatial pattern
formation arise from interaction between Hopf and Tur-
ing modes, and their subharmonics near hte codimension-
two Hopf-Turing bifucation point. Special, it is interest-
ing that spiral wave and travelling wave will appear when
the parameters correspond to the Turing-Hopf bifurca-
tion region III in the spatially extended model (1), i.e.,
the Turing instability and Hopf bifurcation occur simul-
taneously.
III. NUMERICAL RESULTS
The simulation is done in a two-dimensional (2D)
Cartesian coordinate system with a grid size of 600×600.
The fourth order Runger-Kutta integrating method is
applied with a time step ∆t = 0.005 time unit and a
space step ∆x = ∆y = 0.20 length unit. The results
remain the same when the reaction-diffusion equations
were solved numerically in one and two spatial dimen-
sions using a finite-difference approximation for the spa-
tial derivatives and an explicit Euler method for the time
integration. Neumann (zero-flux) boundary conditions
FIG. 1: The sketch map for the bistability and the Hopf bi-
furcation in the system (2) with r = 5.0, a = 5.0, b = 5.0,
m = 0.6, and n = 0.4. The black curve is the g1(p, h). The
colored curves are g2(p, h) with different values of f . The red
curve: f = 0.3; the blue: f = 0.445; the green: f = 0.5; and
the cyan: f = 0.658.
5 10 15
Turing instability
FIG. 2: The sketch map of parameter space (f, ν) bifurcation
diagrams for the spatially extended system (1) with r = 5.0,
a = 5.0, b = 5.0, m = 0.6, dp = 0.05, and n = 0.4.
were emmployed in our simulation. The diffusion terms
in Eqs. (1a) and (1b) often describe the spatial mixing
of species due to self-motion of the organism. The typi-
cal diffusion coefficient of plankton patterns dp is about
0.05, based on the parameters estimatie of Refs [50, 51]
using the relationship between turbulent diffusion and
the scale of the space in the sea. In the previous stud-
ies [10, 11, 29, 45, 46], the authors provided a valueable
insight into the role of spatial pattern for the system (1)
if dp = dh. From the biological meaning, the diffusion
coefficients should satisfy dh ≥ dp. However, in nature
waters it is turbulent diffusion that is supposed to domi-
nate plankton mixing [52], when dh < dp is allowed. The
other reason for choosing such parameter is that it is well-
known new patterns, such as Turing patterns, can emerge
in reaction-diffusion systems in which there is an imbal-
ance between the diffusion coefficients dp and dh [23, 53].
Therefore, we set ν = dh/dp, and investigated whether a
spiral wave would break up into complex spatiotemporal
chaos when the diffusion ratio was varied. Throughout
this paper, we fix dp = 0.05 and dh is a control parameter.
In the following, we will show that the dynamic behav-
ior of the spiral wave qualitatively change as the control
parameter dh increases from zero, i.e., the diffusion ra-
tio ν increases from zero, to more than one. For large
ν (ν > 1), the outwardly rotating spiral wave is com-
pletely stable everywhere, and fills in the space when the
proper parameters are chosen, as shown in Fig. 3(A). Fig-
ure 3(A) shows a series of snapshots of a well-developed
single spiral wave formed spontaneously for the variable
p in system (1). The spiral is initiated on a 600×600 grid
by the cross-field protocol (the initial distribution chosen
in the form of allocated “constant-gradient” perturbation
of the co-existence steady state) and zero boundary con-
ditions are employed for simulations in the two dimen-
sions. From Fig. 3(A) we can see that the well-developed
spiral waves are formed firstly by the evolution. Inside
the domain, new waves emerge, but are evolved by the
spiral wave growing from the center. The spiral wave
can steadily grow and finally prevail over the whole do-
main (a movie illustrating the dynamical evolution for
this case [54] [partly movie−1, movie−2, and movie−3
for dh = 0.2]). Fig. 3(B) shows that the spiral wave
first break up far away from the core center and even-
tually relatively large spiral fragments are surrounded
by a ‘turbulent’ bath remain. The size of the surviv-
ing part of the spiral does not shrink when dh is further
decreasing until finally dh equals to 0, which is different
from phenomenon that is observed previous in the two-
dimensional space Belousov-Zhabotinsky and FitzHugn-
Nagumo oscillatory system [30, 31, 55, 56, 57], in which
the breakup gradually invaded the stable region near the
core center, and finally the spiral wave broke up in the
whole medium. Figure 3(C) is the time sequences (ar-
bitrary units) of the variables p and h at an arbitrary
spatial point within the spiral wave region, from which
we can see that the spiral waves are caused by the ac-
cepted as “phase waves” with substantially group veloc-
ity, phase velocity and sinusoidal oscillation rather than
the relaxational oscillation with large amplitude. This
breakup scenario is similar to the breakup of rotating
spiral waves observed in numerical simulation in chemi-
cal systems [30, 31, 55, 56, 57], and experiments in BZ
systems [58, 59], which shows that spiral wave breakup
in these systems was related to the Eckhaus instability
and more important, the absolute instability.
The corresponding trajectories of the spiral core and
the spiral arm (far away from the core center) at y = 300
are shown in Fig. 4, respectively. From Fig. 4, we can
see that the spiral core is not completely fixed, but oscil-
lates with a large amplitude. However, as dh decreases
to a critical value, an unstable modulation develops in
200 220 240 260 280 300
(D) t (arb. units)
FIG. 3: Well developed spiral waves and some properties of
them. The figures show simulations of the system (1) with
r = 5, a = 5, b = 5, m = 0.6, n = 0.4, dp = 0.05, and
f = 0.3. (A)Well developed spiral waves shown at subsequent
snapshot in time, dh = 0.2. (B) Far-field breakup of the spiral
waves shown at subsequent snapshot in time, dh = 0.002.
The white (black) areas correspond to maximum (minimum)
values of p [Additional movie format available from Ref. [54]].
(C) Oscillations of the variable p and h at an arbitrary spatial
point within the regular spiral wave region for both scenarios.
Each figure is ran the long time until it spatial patterns are
unchange.
regions which is far away from the spiral core (cf. the
middle column of the Fig. 4). These oscillations eventu-
ally grow large enough to cause the spiral arm far away
from the core to breakup into complex multiple spiral
waves, while the core region remains stable (the corre-
sponding movie can be viewed in the online supplemen-
tal in Ref. [54] [partly movie−1 and movie−2, and for
dh = 0.02]). Figures 3(B) and 4(B) show the dynamic
behavior for dh = 0.02, i.e., ν = 0.4. The regular tra-
jectories far away from the core are now the same as the
region of the spatial chaos (cf. the middle column of the
Fig. 4). It is shown that an decrease in the diffusion ra-
tio ν which leads to population oscillations of increasing
amplitude (cf. the left column of the Fig. 4). In the
tradition explain that the minimum value of the popula-
tion density decreases and population extinction becomes
more probable due to stochastic environmental perturba-
tions. However, from the spatial evolution of system (1)
(see Fig. 3), the temporal variations of the density of
different sub-population can become asynchronous and
the events of local extinction can be compensated due to
re-colonization (or diffusion) from other sites.
FIG. 4: The corresponding trajectories (from left to right)
for locations (300, 300), (250, 300), and (50, 300) respectively.
The parameters in (A), and (B) were the same as these in
Fig. 3(A) and (B), respectively.
Furthermore, it is well known that the basic arguments
in spiral stability analysis can be carried out by reducing
the system to one dimensional space [30, 31, 55, 56, 57].
Here we show some essential properties of the spiral
breakup resulting from the numerical simulation. In the
next section we will give the theoretical computation by
using the eigenvalue spectra. In this model, it is worth
noting that we do not neglect the oscillation of the dy-
namics in the core as shown in Fig. 4 due to the system
exhibiting spatial periodic wave trains when the model
is simulated in one-dimensional space. Breakup occurs
first far away from the core (the source of waves). The
spiral wave breaks towards the core until it gets to some
constant distance and then the surviving part of the spi-
ral wave stays stable. These minimal stable wavelengths
are called λmin. So the one-parameter family may be
described by a dispersion curve λ(dh) (see Fig. 5). The
minimal stable wavelength λmin of the spiral wave are
shown in Fig. 5 coming from the simulation in two di-
mensional space. The results of Fig. 5 can be interpreted
as follows: the minimal stable wavelengths decrease with
respect to the decrease of dh but eventually stay at a
relative constant value, which is that the stable spiral
waves are always existing for a larger region values of dh.
Space-time plots at different times are shown in Fig. 6
for two different dh, i.e., different ν, which display the
time evolution of the spiral wave along the cross section
in the two-dimensional images of Fig. 3(A) and (B). As
shown in Fig. 6(A) and (B) for dh = 0.2 and dh = 0.02
respectively, the waves far away from the core display
unstable modulated perturbation due to convective in-
stability [30, 31, 55, 56, 57], but this perturbation is
gradually advected to the left and right sides, and finally
disappears. The instability manifests itself to produce
the wave train breakup several waves from the far-field,
as shown in Figs. 6(B).
FIG. 5: Dependence of the wavelength λmin on the parameter
dh for the system (1) with r = 5.0, a = 5.0, b = 5.0, m = 0.6,
dp = 0.05, and n = 0.4. Note the log scale for dh.
IV. SPECTRA AND NONLINEAR
BIFURCATION OF THE SPIRAL WAVE
In this section, we concentrate on the linear stabil-
ity analysis of spiral wave by using the spectrum the-
ory [56, 60, 61, 62, 63]. From the results in Refs. [56, 62]
we know that the absolute spectrum must be computed
numerically for any given reaction-diffusion systems. In
practice, such computations only require discretization
in one-dimensional space and compare with computing
eigenvalues of the full stability problem on a large do-
main due to the spiral wave exhibitting traveling waves
in the plane (see Fig. 6 about the space-time graphes).
For spiral waves on the unbounded plane, the essential
FIG. 6: Space-time plots of variable p for different time and
dh. The parameters in (A), and (B) are the same as those in
Fig. 3(A) and (B), respectively.
spectrum is also required to compute, since it determined
only by the far-field wave trains of the spiral. The lin-
ear stability spectrum consists of point eigenvalues and
the essential spectrum that is a continuous spectrum for
spiral waves.
For sake of simplicity, the Eqs. (1a) and (1b) can been
written as following
= dp∇
2p+ g1(p, h), (3a)
= dh∇
2h+ g2(p, h). (3b)
Suppose that (p∗, h∗) are a solutions and refer to them
as steady spirals of Eq. (3) that rotate rigidly with a
constant angular velocity ω, and that are asymptotically
periodic along rays in the plane. In a coratating coordi-
nate frame, using the standardized analysis method for
the spiral waves [62, 63], the Eq. (3) is given by
= dp∇
ρ,θp+ ω
+ g1(p
∗, h∗), (4a)
= dh∇
ρ,θh+ ω
+ g2(p
∗, h∗), (4b)
where (ρ, θ) denote polar coordinates, spirals waves are
relative equilibria, then the statianry solutions p∗(ρ, θ)
and h∗(ρ, θ) both are 2π-periodic functions with θ = ϕ−
ωt. In Eqs. (4a) and (4b) the operator∇2ρ,θ denotes ∂ρρ+
A. Computation of spiral spectra
Next, we commpute the leading part of its linear stabil-
ity spectrum for the system (4). Consider the linearized
evolution equation in the rotating frame, the eigenvalue
problem of Eqs. (4a) and (4b) associated with the planar
spiral solutions p∗(ρ, θ) and h∗(ρ, θ) are given by
ρ,θp+ ω
∗, h∗)p+ gh1 (p
∗, h∗)h = λp, (5a)
ρ,θh+ ω
∗, h∗)p+ gh2 (p
∗, h∗)h = λh, (5b)
where g
1 , · · · , g
2 denote the derivatives of the nonlin-
ear functions and g
1(p, h) = r(1 − p) − rp −
(1+bp)2
, gh1 (p, h) = −
2(p, h) =
− abph
(1+bp)2
, and
gh2 (p, h) =
−m− 2fnh
n2+h2
+ 2fnh
(n2+h2)2
. We shall ignore
isolated eigenvalues that belong to the point spectrum,
instabilities caused by point eigenvalues lead to mean-
deringor drifting waves, or to an unstable tip motionin
in excitable media and oscillation media [56, 64, 65, 66].
This phenomenon is not shown in the present paper. In-
stead, we focus on the continuous spectrum that is re-
sponsible for the spiral wave breakup in the far field (see
Fig. 3(b)). By the results in Ref. [62], it turns out that
the boundary of the continuous spectrum depends only
on the limiting equation for ρ → ∞. Thus, we have that
λ is the boundary of the continuous spectrum if, and only
if the limiting equation
ρ,ρp+ ω
∗, h∗)p+ gh1 (p
∗, h∗)h = λp, (6a)
ρ,ρh+ ω
∗, h∗)p+ gh2 (p
∗, h∗)h = λh, (6b)
have solutions p(ρ, θ) and h(ρ, θ) for (ρ, θ) ∈ R+× [0, 2π],
which are bounded but does not decay as ρ → ∞. Since
spiral waves are rotating waves in the plane, the wave
train solutions have the form as u(t, x, y) = u(ρ, ϕ− ωt)
for an appropriate wave numbers k and temporal fre-
quency ω, where we assume that u is 2π-periodic in
its argument so that u(ξ) = u(ξ + 2π) for all ξ and
u = (p, h)T. Spiral waves converge to wave trains
u(ρ, ϕ − ωt) → uwt(kρ + ϕ − ωt) for ρ → ∞, which
are corresponding to asymptotically Archimedean in the
two-dimensional space. Assume that k 6= 0 and ω 6= 0,
and in this case, we can pass from the theoretical frame ρ
to the comoving frame ξ = kρ+ϕ−ωt (ξ ∈ R) in which
the eigenvalue equation (6) becomes
2∇2ξ,ξp+ωpξ+g
1(uwt(ξ))p+g
1 (uwt(ξ))h = λp, (7a)
2∇2ξ,ξh+ ωhξ + g
2(uwt(ξ))p + g
2 (uwt(ξ))h = λh.(7b)
Indeed, any nontrivial solution u(ξ) = (p(ξ), h(ξ))T cor-
responding to the linearization eigenvalue problem (7)
give a solution U(ρ, ·) of the eigenvalue problem for the
temporal period map of (3) in the corotating frame via
U(ρ, ·) = eλtu(kρ− ωt), U(ρ, T ) = eλTu(kρ− 2π).
We write the equations (7) as the first-order systems
= p1,
= h1,
= k−2d−1p
µp− ωp1 − g
1(uwt(ξ))p− g
1 (uwt(ξ))h
= k−2d−1
µh− ωh1 − g
2(uwt(ξ))p− g
2 (uwt(ξ))h
in the radial variable ρ. Then the spatial eigenvalues or
spatial Floquet exponents are deternined as the roots of
the Wronskian
A(λ, k) :=
0 0 1 0
0 0 0 1
(λ− g
1(uwt(ξ))) −
gh1 (uwt(ξ)) −
2(uwt(ξ))
(λ− gh2 (uwt(ξ))) 0 −
where k ∈ R. The function U(ρ, ·) = eλteikρu0(kρ− ωt)
satisfies the equation (3) when the spatial and temporal
exponents ik and λ satisfy the complex dispersion rela-
tion det(A(λ, k) − ik) = 0 for λ ∈ C. We call the ik
in spectrum of A(λ, k) as spatial eigenvalues or spatial
Floquet exponents.
The stability of the spiral waves state (p∗, h∗) on the
plane is determined by the essential spectrum given by
Σess = {λ ∈ C; det(A(λ, k) − ik) = 0 for some k ∈ R}.
Now, we compute the continuous spectrum with the
equation (9) that are parameterized by the wave num-
ber k. For each λ, there are infinitely many stable and
unstable spatial eigenvalues. We plot λ in the complex
plane associated spatial spectrum, see Fig. 7. By the ex-
plaination of Sandstede et al [60], one would know that
if the real part of the essentail spectra is positive, then
the associated eigenmodes grow exponentially toward the
boundary, i.e., they correspond to a far-field instability.
Note that we find the essentail spectra are not sensitive
to temporal frequency, ω.
Re(λ)
K30 K20 K10 0
Im(λ)
Re(λ)
K0.8 K0.6 K0.4 K0.2 0 0.2
Im(λ)
FIG. 7: The essentail spectra of wave trains are obtained by
using the algorithms outlined in Refs. [60, 61]. The param-
eters of (A) and (B) are corresponding to the values used in
the simulations of Fig. 3(A) and (B).
B. Existence and properties of wave trains
Suppose that a reaction-diffusion system on the one-
dimensional space such that the variables equal to a
homogeneous stationary solution. If the homogeneous
steady-state destabilizes, then its linearization accommo-
dates waves of the form ei(kx−ωt) for certain values k and
ω. Typically, near the transition to instability, small spa-
tially periodic travelling waves arise for any wave number
close to kc, which is the critical wavenumber. Their wave
speed is approximately equal to ωc
, where ωc is corre-
sponding to kc. In present paper, we focus exclusively on
the situation where ωc = 0 and kc 6= 0. The bifurcation
with ωc = 0 and kc 6= 0 is known as the Turing bifur-
cation, and the bifurcating spatially periodic steady pat-
terns are often referred to as Turing patterns. Another
class of moved patterns will appear when the instabilities
modulated by Hopf-Turing bifurcation, which is resem-
ble a travelling waves. Moreover, the common feature
of the spiral waves in one-dimensional space mentioned
above is the presence of wave trains which are spatially
periodic travelling waves of the form pwt(kx−ωt; k) and
hwt(kx − ωt; k), where pwt(φ; k) and hwt(φ; k) are 2π-
periodic about φ. Typically, the spatial wavenumber k
and the temporal frequency ω are related via the non-
linear dispersion relation ω = ω(k) so that the phase
velocity is given by
. (12)
A second quantity related to the nonlinear dispersion
relation is the group velocity, cg =
, of the wave
train which also play a central role in the spiral waves.
The group velocity cg gives the speed of propagation of
small localized wave-package perturbations of the wave
train [67]. Here, we are only concerned the existence of
travelling wave solution. In fact, the spiral waves move
at a constant speed outward from the core (see Fig. 6),
so that they have the mathematical form p(x, t) = P (z),
and h(x, t) = H(z) where z = x−cpt. Substituting these
solution forms into Eq. (3) gives the ODEs
+ g1(P,H) = 0, (13a)
+ g2(P,H) = 0. (13b)
Here, we investigate numerically the existence, speed
and wavelength of travelling wave patterns. Our ap-
proach is to use the bifurcation package Matcont 2.4 [68]
to study the pattern ODEs (13). To do this, the most
natural bifurcation parameters are the wave speed cp and
f , but they give no information about the stability of
travelling wave as solutions of the model PDEs (3).
Our starting point is the homogeneous steady state of
Eq. (13) with in the domain III of Fig. 2. The typical bi-
furcation diagrams are illustrated in Fig. 8, which shows
that steady spatially peroidic travelling waves exist for
the larger values of the speed cp, but it is unstable for
small values of cp. The changes in stability occur via
Hopf bifurcation, from which a branch of periodic orbits
emanate. Note that here we use the terms “stable” and
“unstable” as referring to the ODEs system (13) rather
than the model PDEs. Fig. 8(B) illustrates the max-
imun stable wavelength against the bifucation parame-
ter, speed cp, and the small amplitudes have very long
wavelength. It is known that cp =
, hence the tavelling
wave solution exist when the cp 6= 0, i.e., k 6= 0, ω 6= 0.
Using Matcont 2.4 package, it is possible to track the lo-
cus of the Hopf bifurcation points and the Limit point
(fold) bifurcation in a parameter plane, and a typical ex-
ample of this for the cp-f and cp-dh plane are illustrated
in Fig. 9. The travelling wave solutions exist for values
of cp and f lying in left of Hopf bifurcation locus (see
Fig. 9(A)). The same structure about the cp-dh plane
is shown in Fig. 9(B). These reuslts confirm our previ-
ous analysis coming from the algebra computation (see
Fig. 2) and the numerical results (see Fig. 6).
V. CONCLUSIONS AND DISCUSSION
We have investigated a spatially extended plank-
ton ecological system within two-dimensional space and
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
Speed, c
Hopf bifurcation point
FIG. 8: Typical bifurcation diagrams for the pattern
ODEs (13). (A) The spatially periodic travelling waves of
system (3) is existence. The changes in stability occur via
Hopf bifurcation, from which a branch of periodic orbits em-
anate. Thus unstable travelling waves appear. (B) Maxi-
mum stable wavelength along the bifurcation parametercp,
i.e., k 6= 0, ω 6= 0. The parameter values in (A) and (B) are
the same as Fig. 3(A).
found that its spatial patterns exhibit spiral waves dy-
namics and spatial chaos patterns. Specially, the sce-
nario of the spatiotemporal chaos patterns arising from
the far-field breakup is observed. Our research is based
on numerical analysis of a kinematic mimicking the dif-
fusion in the dynamics of marine organisms, coupled to a
two component plankton model on the level of the com-
munity. By increasing (decreasing) the diffusion ratio of
the two variables, the spiral arm first broke up into a
turbulence-like state far away from the core center, but
which do not invade the whole space. From the previous
studies in the Belousov-Zhabotinsky reaction, we know
the reason causing this phenomenon can be illuminated
theoretically by the M. Bär and L. Brusch [30, 31], as
well as by using the spectrum theory that poses by B.
Sandstede, A. Scheel et al [56, 60, 61, 69]. The far-field
breakup can be verified in field observation and is useful
to understand the population dynamics of oceanic ecolog-
ical systems. Such as that under certain conditions the
interplay between wake (or ocean) structures and bio-
logical growth leads to plankton blooms inside mesoscale
hydrodynamic vortices that act as incubators of primary
production. From Fig. 3 and corresponding the movies,
we see that spatial peridic bloom appear in the phyto-
plankton populations, and the details of spatial evolution
of the distribution of the phytoplankton population dur-
ing one bloom cycle, respectively.
In Ref. [70], the authors study the optimal control of
the model (1) from the spatiotemporal chaos to spiral
waves by the parameters for fish predation treated as a
multiplicative control variable. Spatial order emerges in a
range of spatial models of multispecies interactions. Un-
surprisingly, spatial models of multispecies systems often
0 0.5 1 1.5 2 2.5
Speed, c
Locus of Hopf bifurcation points
0.5 1 1.5 2
Speed, c
Locus of Hopf bifurcation points
FIG. 9: An illustration of the variations in parameter space of
the pattern ODEs (13). We plot the loci of Hopf bifurcation
points. (A) f − cp planes; (B) dh − cp planes. The parameter
values in (A) and (B) are the same as Fig. 3(A).
manifests very different behaviors from their mean-field
counterparts. Two important general features of spatial
models of multispecies systems are that they allow the
possibility of global persistence in spite of local extinc-
tions and so are usually more stable than their mean-field
equivalents, and have a tendency to self-organzie spa-
tially or regular spatiotemporal patterns [70, 71]. The
spatial structures produces nonrandom spatial patterns
such as spiral waves and spatiotemporal chaos at scales
much larger than the scale of interaction among individ-
uals level. These structures are not explicitly coded but
emerge from local interaction among individuals and lo-
cal diffusion.
As we know that plankton plays an important role in
the marine ecosystem and the climate, because of their
participation in the global carbon and nitrogen cycle at
the base of the food chain [72]. From the review [73], a
recently developed ecosystem model incorporates differ-
ent phytoplankton functional groups and their competi-
tion for light and multiple nutrients. Simulations of these
models at specific sites to explore future scenarios sug-
gest that global environmental change, including global-
warming-induced changes, will alter phytoplankton com-
munity structure and hence alter global biogeochemical
cycles [74]. The coupling of spatial ecosystem model to
global climate raises again a series of open questions on
the complexity of model and relevant spatial scales. So
the study of spatial model with large-scale is more impor-
tant in the ecological system. Basing on numerical simu-
lation on the spatial model, we can draft that the oceanic
ecological systems show permanent spiral waves and spa-
tiotemporal chaos in large-scale over a range of parame-
ter values dh, which indicates that periodically sustained
plankton blooms in the local area. As with all areas of
evolutionary biology, theoretical development advances
more quickly than does empiraical evidence. The most
powerful empirical approach is to conduct experiments
in which the spatial pattern can be measured directly,
but this is difficulties in the design. However, we can in-
directly measured these phenomenona by the simulation
and compared with the satellite pictures. For example,
the spatiotemporal chaos patterns agree with the per-
spective observation of the Fig. 3 in Ref. [73]. Also, some
satellite imageries [http://oceancolor.gsfc.nasa.gov] have
displayed spiral patterns that represent the phytoplank-
ton [the chlorophyll] biomass and thus demonstrated that
plankton patterns in the ocean occur on much broader
scales and therefore mechanisms thought diffusion should
be considered.
Acknowledgments
This work is supported by the National Natural Sci-
ence Foundation of China under Grant No. 10471040
and the Natural Science Foundation of Shan’xi Province
Grant No. 2006011009.
[1] R. E. Amritkar and Govindan Rangarajan. Spatially
synchronous extinction of species under external forcing.
Phys. Rev. Lett., 96(25):258102, 2006.
[2] Andrzej Pekalski and Michel Droz. Self-organized packs
selection in predator-prey ecosystems. Phys. Rev. E,
73(2):021913, 2006.
[3] Y.-Y. H. Sayama, M. A. M. de Aguiar, and M. Baranger.
Interplay between turing pattern formation and do-
main coarsening in spatially extended population models.
FORMA, 18:19, 2003.
[4] E. Gilad, J. von Hardenberg, A. Provenzale, M. Shachak,
and E. Meron. Ecosystem engineers: From pat-
tern formation to habitat creation. Phys. Rev. Lett.,
93(9):098105, 2004.
[5] Mark J Washenberger, Mauro Mobilia, and Uwe C
Täuber. Influence of local carrying capacity restrictions
on stochastic predator&ndash;prey models. J. Phys.:
Cond. Matt., 19(6), 2007.
[6] Mauro Mobilia, Ivan Georgiev, and Uwe Täuber. Phase
transitions and spatio-temporal fluctuations in stochastic
lattice lotkacvolterra models. J. Stat. Phys., 128(1):447–
483, 2007.
[7] Bernd Blasius, Amit Huppert, and Lewi Stone. Com-
plex dynamics and phase synchronization in spatially ex-
tended ecological systems. Nature, 399(6734):354–359,
1999.
[8] J. von Hardenberg, E. Meron, M. Shachak, and Y. Zarmi.
Diversity of vegetation patterns and desertification.
Phys. Rev. Lett., 87(19):198101, 2001.
[9] A. Provata and G. A. Tsekouras. Spontaneous formation
of dynamical patterns with fractal fronts in the cyclic
lattice lotka-volterra model. Phys. Rev. E, 67(5):056602,
2003.
[10] Alexander B. Medvinsky, Irene A. Tikhonova, Rubin R.
Aliev, Bai-Lian Li, Zhen-Shan Lin, and Horst Malchow.
Patchy environment as a factor of complex plankton dy-
namics. Phys. Rev. E, 64(2):021915, 2001.
[11] Alexander B. Medvinsky, Sergei V. Petrovskii, Irene A.
Tikhonova, Horst Malchow, and Bai-Lian Li. Spatiotem-
poral complexity of plankton and fish dynamics. SIAM
Review, 44:311–370, 2002.
[12] W. S. C. Gurney, A. R. Veitch, I. Cruickshank, and
G. McGeachin. circles and spirals: population persis-
tence in a spatially explicit predatorcprey model. Ecol-
ogy, 79(7):2516–2530, 1998.
[13] J. D. Murray. Mathematical biology. Interdisciplinary
applied mathematics. Springer, New York, 3rd edition,
2002.
[14] Sergei Petrovskii, Bai-Lian Li, and Horst Malchow. Tran-
sition to spatiotemporal chaos can resolve the paradox of
enrichment. Ecological Complexity, 1:37–47, 2004.
[15] Laurence Armi, Dave Hebert, Neil Oakey, James Price,
Philip L. Richardson, Thomas Rossby, and Barry Rud-
dick. The history and decay of a mediterranean salt lens.
Nature, 333(6174):649–651, 1988.
[16] Esa Ranta, Veijo Kaitala, and Per Lundberg. The Spa-
tial Dimension in Population Fluctuations. Science,
278(5343):1621–1623, 1997.
[17] L S Luckinbill. The effects of space and enrichment on a
predator-prey system. Ecology, 55:1142–1147, 1974.
[18] M Holyoak. Effects of nurient enrichment on predator-
prey metapopulation dynamics. J Anim. Ecol., 69:985–
997, 2000.
[19] Vincent A.A. Jansen. Regulation of predator-prey sys-
tems through spatial interactions: a possible solution to
the paradox of enrichment. Oikos, 74:384390, 1995.
[20] Vincent A.A. Jansen and Alun L. Lloyd. Local stabil-
ity analysis of spatially homogeneous solutions of multi-
patch systems. J. Math. Biol., 41:232252, 2000.
[21] Vincent A.A. Jansen. The dynamics of two diffusively
coupled predatorprey populations. Theor. Popul. Biol.,
59:119131, 2001.
[22] J. C. Allen, W. M. Schaffer, and D. Rosko. Chaos reduces
species extinction by amplifying local population noise.
Nature, 364:229232, 1993.
[23] M. C. Cross and P. C. Hohenberg. Pattern formation
outside of equilibrium. Rev. Mod. Phys., 65(3):851, 1993.
[24] Kyoung J. Lee, Raymond E. Goldstein, and Edward C.
Cox. Resetting wave forms in dictyostelium territories.
Phys. Rev. Lett., 87(6):068101, 2001.
[25] Satoshi Sawai, Peter A. Thomason, and Edward C. Cox.
An autoregulatory circuit for long-range self-organization
http://oceancolor.gsfc.nasa.gov
in dictyostelium cell populations. Nature, 433(7023):323–
326, 2005.
[26] Arthur T. Winfree. Varieties of spiral wave behavior:
An experimentalist’s approach to the theory of excitable
media. Chaos, 1(3):303–334, 1991.
[27] V. N. Biktashev, J. Brindley, A. V. Holden, and M. A.
Tsyganov. Pursuit-evasion predator-prey waves in two
spatial dimensions. Chaos, 14(4):988–994, 2004.
[28] M Garvie. Finite-difference schemes for reaction-diffusion
equation modeling predato-prey interactions in matlab.
Bull. Math. Biol., 69:931–956, 2007.
[29] Marcus R. Garvie and Catalin Trenchea. Optimal con-
trol of a nutrient-phytoplankton-zooplankton-fish sys-
tem. SIAM J. Contr. and Opti., 46(3):775–791, 2007.
[30] Markus Bär and Lutz Brusch. Breakup of spiral waves
caused by radial dynamics: Eckhaus and finite wavenum-
ber instabilities. New Journal of Physics, 6:5, 2004.
[31] Markus Bär and Michal Or-Guil. Alternative scenar-
ios of spiral breakup in a reaction-diffusion model with
excitable and oscillatory dynamics. Phys. Rev. Lett.,
82(6):1160–1163, 1999.
[32] Craig R Johnson and Maarten C Boerlijst. Selection at
the level of the community: the importance of spatial
structure. Trends Ecol. and Evol., 17:83–90, 2002.
[33] M Pascual. Diffusion-induced chaos in a spatial preda-
torprey system. Proc. R. Soc. Lond. B, 251:17, 1993.
[34] J. A. Sherratt, B. T. Eagan, and M. A. Lewis. Oscil-
lations and chaos behind predatorprey invasion: math-
ematical artifact or ecological reality? Phil. Trans. R.
Soc. Lond. B, 352:2138, 1997.
[35] S. V. Petrovskii and H. Malchow. Wave of chaos: new
mechanism of pattern formation in spatio-temporal pop-
ulation dynamics. Theor. Popul. Biol., 59:157174, 2001.
[36] J. A. Sherratt. Periodic travelling waves in cyclic preda-
torprey systems. Ecol. Lett., 4:3037, 2001.
[37] N. J. Savill and P. Hogeweg. Spatially induced speciation
prevents extinction: the evolution of dispersal distance in
oscillatory predator-prey models. Proc. R. Soc. Lond. B,
265(1390):25–32, 1998.
[38] E. R. Abraham. The generation of plankton patchiness
by turbulent stirring. Nature, 391:577–580, 1998.
[39] Carol L. Folt and Carolyn W. Burns. Biological drivers
of zooplankton patchiness. Nature, 14:300–305, 1999.
[40] M Scheffer. Should we expect strange attractors behind
plankton dynamics and if so, should we bother? J. Plank-
ton Res., 13:1291–1305, 1991.
[41] I. Hanski, P. Turchin, E. Korplmakl, and H. Henttonen.
Population oscillations of boreal rodents: regulation by
mustelid predators leads to chaos. Nature, 364:232235,
1993.
[42] S. Ellner and P. Turchin. Chaos in a noisy world: new
methods and evidence from time-series analysis. Am.
Nat., 145:343375, 1995.
[43] B. Dennis, R. A. Desharnais, J. M. Cushing, S. M. Hen-
son, and R. F. Costantino. Estimating chaos and com-
plex dynamics in an insect population. Ecol. Monogr.,
71:277303, 2001.
[44] M Scheffer. Fish and nutrients interplay determines algal
biomass: A minimal model. Oikos, 62:271–282, 1991.
[45] H. Malchow. Spatio-temporal pattern formation in non-
linear non-equilibrium plankton dynamics. Procc. R. Soc.
Lond. B, 251:103, 1993.
[46] M. Pascual. Diffusion-induced chaos in a spatial
predator-prey system. Procc. R. Soc. Lond. B, 251:1–7,
1993.
[47] Arnd Scheel. Radialy symmetric patterns of reaction-
diffusion systems. Mem. Amer. Math. Soc., 165:86, 2003.
[48] R A Satnoianu, M Menzinger, and P K Maini. Turing
instabilities in general system. J. Math. Biol., 41:493–
512, 2000.
[49] Quan-Xing Liu, Bai-Lian Li, and Zhen Jin. Resonant pat-
terns and frequency-locked induced by additive noise and
periodically forced in phytoplankton-zooplankton sys-
tem, 2007.
[50] Sven Erik Jørgensen and G. Bendoricchio. Fundamentals
of ecological modelling. Developments in environmental
modelling; 21. Elsevier, Amsterdam; New York, 3rd edi-
tion, 2001.
[51] Akira Okubo. Diffusion and ecological problems: mathe-
matical models. Biomathematics; v. 10. Springer-Verlag,
Berlin; New York, 1980.
[52] George Sugihara and Robert M. May. Nonlinear forecast-
ing as a way of distinguishing chaos from measurement
error in time series. Nature, 344(6268):734–741, 1990.
[53] A. M. Turing. The chemical basis of morphogenesis.
Philosophical Transactions of the Royal Society of Lon-
don. Series B, Biological Sciences, 237(641):37–72, 1952.
[55] Fagen Xie, Dongzhu Xie, and James N. Weiss. Inwardly
rotating spiral wave breakup in oscillatory reaction-
diffusion media. Phys. Rev. E, 74(2):026107, 2006.
[56] Björn Sandstede and Arnd Scheel. Absolute versus
convective instability of spiral waves. Phys. Rev. E,
62(6):7708–7714, 2000.
[57] S. M. Tobias and E. Knobloch. Breakup of spiral waves
into chemical turbulence. Phys. Rev. Lett., 80(21):4811–
4814, 1998.
[58] Q. Ouyang and J. M. Flesselles. Transition from spirals
to defect turbulence driven by a convective instability.
Nature, 379(6561):143–146, 1996.
[59] Qi Ouyang, H. L. Swinney, and G. Li. Transition from
spirals to defect-mediated turbulence driven by a doppler
instability. Phys. Rev. Lett., 84(5):1047–1050, 2000.
[60] Björn Sandstede and Arnd Scheel. Curvature effects
on spiral spectra: Generation of point eigenvalues near
branch points. Phys. Rev. E, 73:016217, 2006.
[61] Jens D.M. Rademacher, Björn Dandstede, and Arnd
Scheel. Computing absolute and essential spectra using
continuation. Physics D, 229:166–183, 2007.
[62] P. Wheeler and D. Barkley. Computation of Spiral Spec-
tra. SIAM J Appl. Dynam. Syst., 2006.
[63] Björn Sandstede and Arnd Scheel. Absolute and convec-
tive instabilities of waves on unbounded and large bound
domians. Physics D, 145:233–277, 2000.
[64] Dwight Barkley. Linear stability analysis of rotat-
ing spiral waves in excitable media. Phys. Rev. Lett.,
68(13):2090–2093, 1992.
[65] Dwight Barkley. Euclidean symmetry and the dynamics
of rotating spiral waves. Phys. Rev. Lett., 72(1):164–167,
1994.
[66] Igor Aranson, Lorenz Kramer, and Andreas Weber. Core
instability and spatiotemporal intermittency of spiral
waves in oscillatory media. Phys. Rev. Lett., 72(15):2316–
2319, 1994.
[67] Björn Sandstede and Arnd Scheel. Defects in oscillatory
media: Toward a classification. SIAM J. Appl. Dynam.
Syst., 3(1):1–68, 2004.
[68] A Dhooge, W Govaerts, Yu A Kuznetsov, W Mestrom,
A M Riet, and B Sautois. Matcont and Cl-Matcont: Con-
tinuation toolboxes in Matlab. Utrecht University, The
Netherlands, 2006.
[69] Paul Wheeler and Dwight Barkley. Computation of spiral
spectra. SIAM J. Appl. Dynam. Syst., 5:157–177, 2006.
[70] Maarten C. Boerlijst. The Geometry of Ecological Inter-
actions: Simplifying Spatial Complexity, chapter Spirals
and spots: Novel Evolutionary Phenomena through spa-
tial self-structuring, pages 171–182. Cambridge Univer-
sity Press, 2000.
[71] Ulf Dieckmann, Richard Law, and Johan A J Metz. The
Geometry of Ecological Interactions: Simplifying Spatial
Complexity. Cambridge University Press, 2000.
[72] J. Duinker and G. Wefer. Das co2-problem und die rolle
des ozeans. Naturwissenschaften, 81(6):237–242, 1994.
[73] M. Pascual. Computational ecology: From the complex
to the simple and back. Plos Comput. Biol., 1(2):101–
105, 2005.
[74] E. Litchman, C. A. Klausmeier, J. R. Miller, O. M.
Schofield, and P. G. Falkowski. Multi-nutrient, multi-
group model of present and future oceanic phytoplankton
communities. Biogeosciences, 3(4):585–606, 2006.
