Effect of demographic noise in a phytoplankton-zooplankton model of bloom dynamics

0.33 MB·English

by  Piero Olla

Demographic noise as a triggering mechanism for algal blooms
Piero Olla
ISAC-CNR and INFN, Sez. Cagliari, I–09042 Monserrato, Italy.
(Dated: July 27, 2012)
Population models, such as those for plankton dynamics, are often based on a mean-field ap-
proximation of individual behaviors. A weakly stable mean-field configuration, however, can be
destabilized by demographic noise. In certain cases, such destabilization persists even in the ther-
modynamiclimit. Thepossibilerelevanceofsuchaneffecttotheonsetofalgalbloomsisdiscussed.
PACSnumbers: 87.23.Cc,87.10.Mn,02.50.Ey,05.40.-a
2 I. INTRODUCTION zooplankton control and grows to the carrying capacity
1
of the medium, before also the zooplankton population
0
grows appreciably, and is able to bring that of the phy-
2 Phytoplankton provide the basis of the food chain in
l the ocean. They are constituted for the large part by toplankton back to its pre-bloom level.
u At a sufficiently coarse grained scale, the system dy-
microscopic algae and their abundance is such that they
J namicscanbedescribedbyphytoplanktonandzooplank-
6 are responsible for roughly one half of the total photo- tonconcentrationfieldsP¯(x,t)andZ¯(x,t)[17,18]. This
synthesis going on in the planet [1].
2 constitutes a mean field approximation for the stochas-
Although phytoplankton are present in virtually all of tic individualdynamics; at differentlevels of complexity,

] the world hydrosphere (more precisely, the top 100m suchadescriptionmayincludethe nutrientandeventhe
E
layer of the water column, called the euphotic layer, detritus concentration fields N¯(x,t) and D¯(x,t) (NPZ
P
wherethereisenoughlightforphotosynthesis),theirdis-
[19, 20] and NPZD [21] models). The tininess of both
.
o tributioninspaceandtimeisfarfromuniform[2]. Space phyto-andzooplanktonindividuals,aswellasthesizeof
i patchiness is observed at scales that go from that of the
b typicalareasofinterestforthestudyoftheconcentration
individualtoseveralhundredskilometers[3,4]. Thespa-
- dynamics (at leastseveralmeters), suggeststhat a mean
q tial patterns, revealed e.g. by remote sensing, include
fielddescriptionisindeedthemostappropriate. Thefact
[ filaments, fronts and more irregular shapes, suggesting
thatbloomeventsmaybetriggeredbyfluctuationsinthe
that turbulent transport by sea currents may play an
1 system, however, should suggest caution.
v important role [5–9]. The point that we want to examine in the present pa-
4 Variabilityintime, onthe otherhand,ischaracterized per is precisely whether microscopic fluctuations, disre-
5
by sporadic bloom events in which the plankton density gardedinameanfieldapproach,canresurfaceatmacro-
2
in extended regions can grow by up to three orders of scopicscaleandcontributetotriggerbloomevents. This
6
. magnitude in few days [10]. These are seasonal events isanexampleofdemographicfluctuationinducedbreak-
7
that may or may not recur annually, and are superim- up of mean field descriptions of reaction-diffusion sys-
0
posedtotheweakerday-nightandseasonalcycles. Their tems, a type of phenomenon that have received a great
2
1 importance is utmost for several reasons. Depending on deal of attention in recent years [22–26]
: thespecies(e.g. diatoms),duetotheirabundance,bloom In the present paper we shall consider a simple PZ
v
events can have consequences even at the level of the model due to Truscott and Brindley (TB model), which
i
X globalbiogeochemicalcycle [11]. Bloomevents involving has the nice characteristics that the plankton behaves
r other species (e.g. dinoflagellates) can have harmful ef- likeanexcitablemedium[18]. Theoriginalmodeliszero-
a fects on water quality and fishery [12, 13]. Being able dimensional,butcanbeeasilygeneralizedtoincludespa-
to predict the onset of algal blooms is clearly an issue of tial effects such as advection and diffusion [27]. In the
great practical importance. present analysis a spatially homogeneous domain will be
Several mechanisms of both biological and physical considered with no advection, but with diffusive terms
nature contribute to blooms, not all of them probably accounting for small scale motions in the water column.
known yet. From the point of view of biology, the domi- The game to play will be to determine the demographic
nant controlling factor for phytoplankton growth is usu- fluctuations implicitly neglected in the TB model, and
allyconsideredtobegrazingbyzooplankton[14](though see if they candestabilize the systemevenat the levelof
nutrient transport, light, and the presence of a thermo- spatial averages.
cline surely play an important role [15]). This is the so This paper is organized as follows. In Sec. II, the
called mismatch issue [16]: since the life cycle of phyto- main properties of the TB model are reviewed. In Sec.
planktonistypicallyoneorderofmagnitudeshorterthan III, the expression for the stochastic contribution to de-
thatofzooplankton,apositivefluctuationinphytoplank- mography are derived. In Sections IV and V, the effect
ton productivity is not compensated by a simultaneous of demographic noise on the dynamics of the TB model,
increase in the zooplankton population. This produces without and with seasonal forcing, is analyzed. Section
the bloom event: the phytoplankton population escapes VI is devoted to the conclusions. Some technical details
2
on the master equation treatment of demographic fluc-
3.5
tuations, in spatially extended domains, are provided in
the Appendix for reference. 3
2.5
II. THE TRUSCOTT-BRINDLEY MODEL
¯
Z/α 2
We review here the main properties of the TB model. 1.5
Forthemomentwedisregardthespatialstructureofthe
fields and focus on the original zero-dimensional version 1
•
of the model, which is described by the equations:
0.5
P¯˙ = r P¯ 1− P¯ −R P¯2Z¯ , 0 2 4 6 P¯8/α10 12 14 16
0 (cid:16) K(cid:17) mP¯2+α2
Z¯˙ = −µZ¯+γRmP¯P2¯+2Z¯α2. (1) (FfIaGt.do1t: tTorlaojwecetrolreieftscaoprnperoraocfhipnigctuthree)fistxaerdtinpgoifnrtom(Pdf,iffZefr)-
ent initial conditions; going from outer to inner trajectory:
The main characteristics of the dynamics are the follow-
(P¯,Z¯)=(0.5,0.5); (0.5.0.6); (0.6,0.66); (0.6,0.68). The pa-
ing:
rametersarethoseofEq. (2). Thedottedlineindicatechange
• Logistic reproductive behavior of the phytoplank- of signature of the Jacobian matrix for Eq. (1), and could
ton, with the carrying capacity K determining the roughly be identified as thepoint where thebloom begins.
maximumconcentrationthatthe mediumcansup-
port at steady state.
We see that the dynamics is characterized by two inde-
• Grazing by zooplankton characterized by a so
pendent small parameters: γ and ǫ. The dynamics de-
called Holling-III kind of behavior [28]. Zooplank-
scribed by Eq. (1) has a fixed point at the scale of the
ton are able to graze on phytoplankton with opti-
malrateR Z¯,onlyiftheconcentrationofthe sec- half-saturation constant α:
m
ond is above the level fixed by the half-saturation q
P = α ,
concentration α. Below this threshold, the grazing f r1−q
rateisquadraticinP¯: R (P¯/α)2Z¯,reflectingboth
m rˆ P αrˆ
areducedgrazingabilityofthezooplanktoninadi- Z = 0 (1− f)(α2+P2)≃ 0 . (4)
lute environment, and the actual reduced amount f Pf K f q(1−q)
of food available. p
For small γ,ǫ ≪ 1 and q < 1/2, this fixed point is glob-
ally attracting (a Holling-III functional form for grazing
• Acarryingcapacityofthe medium supposedmuch
appears to be crucial for stability). However, if the ini-
largerthanthehalfsaturationconcentration,K ≫
α, meaning that at high values of P¯, the ability tialzooplanktonconcentrationistoolow,beforereaching
thefixedpoint,the systemwillmakeanexcursiontothe
of zooplankton to control phytoplankton growth is
highP¯ ∼K range,whichcouldbeinterpretedasabloom
limited.
event. The situation is illustrated in Fig. 1. As shown
• A zooplankton reproductive dynamics supposed in figure, the onset of bloom could roughly be identified
slowerthanthatofthephytoplankton: µ/Rm ≪1, in the P¯Z¯ plane by the line where the largesteigenvalue
while r0/Rm ∼ 1. A conversion efficiency γ as- of the Jacobian of Eq. (1) crosses to positive, and phase
sumed consistently small. points start to separate exponentially.
The above picture of bloomtriggering by zooplankton
The defaultvaluesofthe constantsthatareutilizedin
depletioncanbeimprovedincludingtheeffectofseasonal
the zero-dimensional case are [18]:
forcing. Model equations (1) can accommodate this ef-
r =0.3/day, R =0.7/day, µ=0.012/day, fect by letting the phytoplankton productivity r become
0 m
K =108mgC/m3, α=5.7mgC/m3, (2) dependentonthe temperature, assuggestedin[29]. The
parameterizationthat we adopt is the same as in [29]: a
γ =0.05,
Van’t Hoff kind of dependence for r [30]:
where the units “mgC” stand for carbon milligrams in r →r(T)=r 2υ(T−T0) (5)
dry weight. Fromhere we canextractthree independent 0 0
dimensionless groups andasinusoidaldependence ontime ofthe temperature:
r µ
rˆ0 = 0 ≃0.43, q = ≃0.34, T(t)=T0+∆T sin(Ωt+φ). (6)
R γR
m m
α withΩ=2π/(365days)toallowforanannualcycle,and
ǫ= ≃0.053. (3) φ/(2π) = 0.59, to have that setting t = 0 on January
K
3
γ
0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
12 18
11 16
∆
10 12T 16
q
Z¯/α o∆T(C)γ 6789 0.1 0.2 0q.3 0.4 48 (C)o 1124ǫ∆T(Co
5 10
1 )
• 4
3 8
2
6
1
0.1 1 10 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
¯ ǫ
P/α
FIG. 2: Bloom and no-bloom cycles (heavy and thin lines FIG. 3: Dependenceon the critical forcing amplitude ∆Tcrit
respectively) in the seasonally forced TB model. The phase ontheparametersγ (thinline),ǫ(heavyline) andq (insert).
pointscirclecounterclockwisealongthecycle. Noticethatthe In the three cases ∆Tγ,ǫ,q indicate the value of ∆Tcrit ob-
intersections of the two cycles corresponds to phese points tained varying γ,ǫ,q respectively and keeping the remaining
that in the two cycles are associated to different (although parameters fixed to thevalues in Eq. (3).
close)instantsoftime. Thedottedverticallinegivestheposi-
tionofthefixedpointsthatwouldcorrespondtothedifferent
values taken by rˆduring the year. The fat dot still identifies
the fixed point corresponding to rˆ= rˆ0. The parameters in ical phytoplankter. A reasonable estimate (with large
thegraph are those in Eqs. (5-6). variations) for the mass of a copepod could be, for in-
stance:
m ∼20µgC, (7)
Z
1st, causes the first temperature minimum to occur on
March 1st and the first maximum on August 29th. For correspondingto a typicalindividual size in the millime-
the sake of definitenes, as in ..., we set ∆T = 6oC and ter range [6, 31]. In a no-bloom regime P¯,Z¯ ∼α with α
υ =0.1oC−1. given as in Eq. (2), Eq. (7) would lead to a numerical
Addingaseasonalforcing,turnsouttomodifythe dy- density of the order of one zooplankton individual per
namics in important way, with the single fixed point in liter. Given the much smaller size of microscopic algae,
the autonomous case leaving way to two stable limit cy- phytoplankton could in turn be treated as a continuum
cles [29]. The situation is illustrated in Fig. 2: a small at those scales. We thus expect that the demographic
amplitude no-bloom cycle coexists with a large ampli- fluctuations in the systembe drivenby the zooplankton.
tude bloom cycle, each one characterized by a well de- Locally,demographicfluctuationswillbetheresultofa
fined(time-dependent)basinofattraction. Asillustrated competitionbetween stochasticityatthe individual level
in Fig. 3, the small cycle will remain stable only if of the birth-death process, and mixing by spatial trans-
the seasonal temperature excursion ∆T is below a crit- port. We shall consider a two-dimensional situation, in
ical threshold ∆T (rˆ ,q,ǫ,γ), that, for the values of which the vertical structure of the water column is not
crit 0
the parameters quoted in Eqs. (2) is ∆T ≃ 6.1oC. resolved, and parameterize horizontal mixing through a
crit
Similarly, it can be shown that, if ∆T is too small, only diffusivity κ, supposed equal for both phyto- and zoo-
the no-bloom cycle will survive. The no-bloom destabi- plankton. This is probably the only viable strategy to
lizationthresholdis verycloseto the value ofthe forcing describe a very complex situation, in which microscopic
∆T =6oCconsideredin[29]. The analysisinthatpaper swimming, stirring by larger organisms, turbulence gen-
showed in fact that addition of a fluctuating component erated by perturbations at the water surface, all play an
to the forcing, lead to random switch from year to year important role [32]. Likewise, we shall neglect all effects
between the two regimes. As it is clear from Fig. 2, the oflargescaleadvection,includingthe forcinginducedby
destabilizationislikelytotakeplaceneartheintersection the formation of fronts, where the two plankton popula-
betweenthe bloomandno-bloomtrajectories,where the tions may get out of balance [5].
separationbetweenthephasepointsofthetwocycles(at Including the possibility of a 2D spatial structure, the
equal times) is smaller. original model equations (1) can be written in the form
P¯˙ = B (P¯,Z¯)−D (P¯,Z¯)+κ∇2P¯,
P P
III. DEMOGRAPHIC FLUCTUATIONS Z¯˙ = B (P¯,Z¯)−D (P¯,Z¯)+κ∇2Z¯, (8)
Z Z
We consider a situation in which the typical size of a whereB andD givethelocalbirthanddeathrates
PZ PZ
zooplanktonindividualismuchlargerthanthatofatyp- in the population, and P¯ =P¯(x,t) and Z¯ =Z¯(x,t) give
4
the plankton concentration averaged over the height of domain):
the water column. Expressing time in units R−1 and
m
concentrations in units α (and lengths in terms of some ∂ρ = Π ρ+Ω1/2 P¯˙ ∂ρ +Z¯˙ ∂ρ
reference scale), the parameters entering Eq. (8) can be ∂t Xi (cid:26) i Z h i∂φi i∂ζii
written in the form
+ Ω exp −Ω−1/2 ∂ −1 B
P C ∂φ P,i
nh (cid:16) i(cid:17) i
BP =rˆP¯; DP =ǫrˆP¯2+P¯2Z¯/(1+P¯2); −1/2 ∂
+ exp Ω −1 D ρ
B =γP¯2Z¯/(1+P¯2); D =γqZ¯. (9) C ∂φ P,i
Z Z h (cid:16) i(cid:17) i o
+ Ω exp −Ω−1/2 ∂ −1 B
Z Z ∂ζ Z,i
From now on, unless otherwise stated, all relations will nh (cid:16) i(cid:17) i
be expressed in dimensionless form. + exp Ω−1/2 ∂ −1 D ρ , (12)
Rather than working in a field theoretical setting [33], h (cid:16) Z ∂ζi(cid:17) i Z,io (cid:27)
we prefer to apply the standard system size expansion where Ω1/2 = Ω−1/2Ω , and the additional term Π ρ
approachofvanKampendirectlytothemasterequation C Z P i
accounts for the exchange of plankton between adja-
for the system [34]. Let us therefore partition the do-
cent volumes produced by diffusion (see Appendix). Ex-
maininvolumes ofhorizontalsize ∆x, eachonecontain-
panding to lowest order in Ω−1 , both the exponen-
ing instantaneously N ,N individuals of the P and Z P,Z
P Z
tials and the reaction rates in Eq. (12) [recall that
groups. We can introduce instantaneous concentrations
P,Z =Ω−P,1ZNP,Z, coarsegrainedathorizontalscale ∆x, BP = BP(P¯ +ΩZ−1/2φ,Z¯+ΩZ−1/2ζ),...], we obtain the
with Fokker-Planckequation:
∂ρ ∂
+ a φ +a ζ ρ
ΩP,Z =h(∆x)2/mP,Z (10) ∂t Xi n∂φi(cid:16) PP,i i PZ,i i(cid:17)
∂
+ a φ +a ζ ρ
parameterizingthesizeofthepopulationinthevolumes, ∂ζ ZP,i i ZZ,i i
i(cid:16) (cid:17) o
andhindicatingtheheightofthewatercolumn. Indicate 1 ∂ ∂
with P˜ = P −P¯ and Z˜ = Z −Z¯ the fluctuating part of = Πiρ+ Ξijρ, (13)
2 ∂ζ ∂ζ
the coarse-grained field; in the present situation of fluc- Xi Xij i j
tuations driven by zooplankton discreteness, we expect
P˜/P¯,Z˜/Z¯ = O(N−1/2). Let us define normalized fluc- whereΞij =[BZ(P¯i,Z¯i)+DZ(P¯i,Z¯i)]δij,andtheaPP,...
Z give the entries of the Jacobian matrix for Eq. (1)
tuation fields φ,ζ with this scaling contribution factored
(a ≡ ∂P¯˙/∂P¯,...). Notice that the only noise corre-
out: PP
lator present in the equation, Ξ , is the one associated
ij
with the variable ζ, i.e. with the zooplankton fluctua-
N = Ω P :=Ω (P¯+Ω−1/2φ), tions.
P P P Z
N = Ω Z :=Ω (Z¯+Ω−1/2ζ). Going back to the original variables P˜,Z˜, and tak-
Z Z Z Z ing the continuous limit, we see that the Fokker-Planck
equation (13) is equivalent to the system of Langevin
equations
Thestandardapproachistohypothesizethatthebirth
arenfldecdteabtirhthraatnesdBdePa,ZtharnadteDsPa,tZthaet tinhdeivpiodpuuallatleiovnel,sccaolne-, ∂∂Pt˜ +aPPP˜+aPZZ˜ =κ∇2P˜
ditionedtotheinstantaneousvaluesoftheconcentration
∂Z˜
fields P(x,t) and Z(x,t). At the population level, the +a P˜+a Z˜ =κ∇2Z˜+ξ, (14)
∂t ZP ZZ
transition probabilities dP in an interval dt will be:
where
NP →NP +1: dP =ΩPBP(P,Z)dt, hξ(x,t)ξ(0,0)i=Ξ(x,t)δ(x)δ(t),
NP →NP −1: dP =ΩPDP(P,Z)dt, Ξ(x,t)=mˆ [B (P¯,Z¯)+D (P¯,Z¯)], (15)
Z Z Z
N →N +1: dP =Ω B (P,Z)dt,
Z Z Z Z
NZ →NZ −1: dP =ΩZDZ(P,Z)dt. (11) Za¯n(dx,wt)e(hmaovreepduettamiˆlsZin=thmeZA/php,enP¯di≡x).P¯(x,t) and Z¯ ≡
Equations(13-15)provideamesoscopicdescriptionfor
A standard procedure [34] leads then, from Eq. (11), to the system dynamics that is valid as long as the fluctu-
the master equation for the probability density function ations at scale ∆x can be considered small. It is imme-
(PDF)ρ({φ ,ζ },t)(theindexilabelsthevolumesinthe diately clear that this condition cannot be satisfied in a
i i
5
bloom regime. In fact, as illustrated in Fig. 1, bloom At small scales, the fluctuations are smeared out by dif-
is associated with change of signature in the Jacobian fusion, with the asymptotic law CPPk ≃ a2PPΞ/(κk2)3,
matrix ofEq. (1)andwith exponentialgrowthoffluctu- CPZk ≃ aPZΞ/(κk2)2, CPZk ≃ Ξ/(κk2). The transition
ationsinEq. (14). WeshalldealwiththisregimeinSec. occurs at κk2 ∼ a a /a , which sets the crossover
ZP PZ PP
V; let us concentrate for the moment the case of small length,fromEqs. (3)and(17),backtodimensionalunits:
fluctuations.
λ = κ/µ. (21)
c
p
Thisisthetypicaldistancetravelledbyazooplankterina
IV. DYNAMICS NEAR THE FIXED POINT
lifetimeandcorrespondstothecharacteristicwavelength
of the chemical waves supported by the system in the
Let us consider first the case of a TB model without mean field. Notice that, contrary to the case considered
seasonal forcing, and study the fluctuations around the in [27], for the choice of parameters utilized, no Turing
fixed point given by Eq. (4). The analysis is similar instability is present.
to the one carried on in [26] on another PZ model (the Thecorrelationspectrumthathasbeenobtained,char-
Levin-Segel model [35]), which focused on the destabi- acterizedbyaplateauatkλ <1,andadecayatkλ >1,
c c
lization of Turing patterns. The evolution equations for corresponds to fluctuations with a correlation scale λ .
c
thecorrelationfunctionsCPP(x,t)=hP˜(x,t)P˜(0,t)i,... Thefluctuationamplitudecanbeestimatedapproximat-
can be obtained fromEqs. (14-15),and take the form in ing the decay at kλ > 1 with a step function and ap-
c
Fourier space: proximatingthesolutionforkλ <1withEq. (20). This
c
gives for the fluctuation amplitude
1
2C˙PPk+(aPP +κk2)CPPk+aPZCPZk =0 CPP(0)=Fx−=10[CPPk]∼λ−c2CPP0,
C˙PZk+aZPCPPk+(aPP +aZZ +2κk2)CPZk with CPP0 as given in Eq. (20), and similar expressions
+aPZCZZk =0 for CPZ(0) and CZZ(0). From Eqs. (20) we get for the
1 ratio of the fluctuation amplitude to the mean (back to
2C˙ZZk+aZPCPZk+(aZZ +κk2)CZZk =Ξ. (16) dimensional units):
At the fixed point we find immediately, using Eq. (9) CPP(0) CPZ(0) CZZ(0) mZ
∼ ∼ ∼ , (22)
and working to lowest order in ǫ and γ: P¯2 P¯Z¯ Z¯2 αhλ2
c
a ≃−rˆ(1−2q), a =−q, that is the ratio of the zooplankter mass and the typical
PP PZ
total plankton mass in a water column of height h and
a ≃2γrˆ(1−q), a =0,
ZP ZZ horizontal extension λ .
c
Ξ≃2γmˆ q/(1−q). (17)
Z
p
Equation (16) becomes at steady state: V. DESTABILIZATION OF THE NO-BLOOM
REGIME
(aPP +κk2)CPPk+aPZCPZk =0,
aZPCPPk+(aPP +2κk2)CPZk+aPZCZZk =0, Let us pass to consider a seasonally forced situation
aZPCPZk+κk2CZZk =Ξ, andaskwhether demographicnoiseis ableto destabilize
the global bloom and no-bloom cycles.
that has solution
As already discussed, the linearized theory of Sec. IV
CPPk ≃ ∆−1a2PZΞ, cannot be utilized in the present case, as the trajecto-
CPZk ≃ ∆−1aPZ(−aPP +κk2)Ξ, rTiehsuse,vonluvemfeorricmalossotluoftitohneiorftitmhee imnaasnteurnesqtaubalteiorneg(i1o2n).
CZZk ≃ ∆−1(−aPP +κk2)2Ξ, (18) becomes necessary, either by direct means, or by Mon-
tecarlo techniques. In alternative, one may try to ex-
with
tendtheLangevinequation(14)beyondthelinearregime
∆=a a a +a2 κk2−2a κ2k4+κ3k6. (19) by replacing the Jacobian matrix a with the full RHS
PP PZ ZP PP PP ij
(right hand side) of the mean field TB model (8). It
FromEqs. (18-19),weseethatthereisalongwavelength
is convenient to express lengths in units λ , so that the
c
range, dominated by demography; using Eq. (17):
forced equation can be written in the form
CPPk ≃ rˆ02(1−qqmˆ)(Z1−2q)r1−q q, PZ˙˙ == BBPZ((PP,,ZZ))−−DDZP((PP,,ZZ))++λλqq∇∇22ZP,+ξ. (23)
mˆ q
CPZk ≃ Z , [Notice the dependence of the terms on RHS of Eq. (23)
rˆ (1−q)r1−q
0 on the fluctuating fields P and Z; similarly, the depen-
mˆZ(1−2q) q dence on P¯ and Z¯ in Eq. (15) is replaced by one on P
CZZk ≃ q(1−q) r1−q. (20) and Z].
6
10−2
♦♦♦♦
♦♦
♦♦
10−3 ♦♦♦♦♦♦♦
mˆZ ♦♦♦
♦
10−4 ♦
♦
10−5 ♦
5 5.2 5.4 5.6 5.8 6
∆T(oC)
FIG. 4: No-bloom cycle destabilization threshold in func-
tion of ∆T, from numerical simulation in a periodic domain
200∆X ×200∆X with ∆X = λ /10. The no-bloom cycle
c
is considered destabilized if the system crosses the threshold
P¯ = 10 before t = 10years. All parameters except ∆T and
mˆ set to the values in Eqs. (2) and (7). Initial conditions
Z FIG. 5: Snapshots of the phytoplankton concentration field
set equal to(Pf,Zf) uniformly in the domain. from a simulation with mˆ = 10−4; ∆T = 6oC and other
Z
parameters set as in Eqs. (2) and (7). Size of the do-
main 200∆X ×200∆X; periodic boundary conditions, with
∆X = 0.3λ . Initial conditions set equal to (P ,Z ) uni-
c f f
Inthis way,theroleofeffectivediffusivity isplayedby formly in the domain. A: day 700 (December 1st: no-bloom
theproductλq,whiletheparametermˆZ,which,through condition); B: day 890 (May 9th: pre-bloom condition); C:
Eq. (15), determines the amplitude of the noise ξ, takes day 924 (June13th: bloom peak); D: day 970 (July 29: con-
the form, in terms of dimensional parameters: centration minimum after bloom).
m µ
Z
mˆ = . (24)
Z
αhκ
This quantity will play the role of a control parameter
for the theory. [We notice by the way that mˆ coincides
Z
with the amplitude ratio in Eq. (22)].
The approximation leading to Eq. (23) is clearly un-
controlled, as no account is taken of the modifications
producedinthenoiseterm. Nevertheless,exceptforvery
largevaluesofmˆ ,MontecarlosimulationwithEq. (11),
Z
and direct solution of Eq. (23), have given very close re-
sults. This appears to be due the fact that the noise
plays a crucial role in destabilizing the cycles, only in
the narrow regions around the “intersections” (see Fig.
2). In the unstable regions, the dynamics is dominated
by amplification of the separation between phase points
produced by the deterministic part of the equation; the
noise plays in the unstable regions a minor role. Thus,
the fact that inthese regionsthe noise in Eq. (23)is not
well approximated, does not produce dangerous effects.
As expected, a sufficiently high level of noise destabi- FIG. 6: Sameas Fig. 5 in thecase of thezooplankton field.
lizesthesmallcycleandleadstolockingthesysteminthe
large bloom cycle. As illustrated in Fig. 4, the thresh-
old in mˆ becomes lower as the critical forcing ∆T is
Z crit
approached. with a persistence time between change of directions of
For default values of the parameters, such as those in the order of one second [6].
Eqs. (2)and(7),with∆T =6oC,thethresholdwouldbe AsillustratedinFig. 5,evenwhenthesystemislocked
atmˆ ≃4·10−5, which, fora depth h≃5m,wouldcor- globally on a bloom cycle, its dynamics is characterized
Z
respond to a diffusivity κ ≃ 0.2m2/day ≡ 0.024cm2/s byspatialfluctuations,withregionsofsize∼λ inwhich,
c
and a correlation length λ ≃ 4.1m. In comparison, during bloom events, P remains well below its typical
c
κ ∼ 0.01cm2/s would be the diffusivity that would be bloom value.
produced by microscopic swimming at speed ∼ 1mm/s This picture is confirmed by looking at the evolution
7
16 destabilizationcouldbejustafinitesizeeffect,andwould
14
12 disappear in the thermodynamic limit. One suggestion
¯Z10
8 that this is not the case comes from the invasive char-
, 6
¯P 4 acter of the destabilization phenomenon. The sequence
2
0 B−C inFig. 5givesahintoftheprocess: regionsofsize
0 1 2 3 4 5 6 7 8 9 10 ∼λ ,characterizedbyhighvaluesofP,throughdiffusive
2 c
Z1.5 coupling,destabilizenearbyregionswithlowvaluesofP,
ˆcP 1 and push them in the bloom cycle. This picture is cor-
,Z0.5 roborated by the behavior of the system in the noiseless
P, 0 mˆZ =0case. AsshowninFig. 8,intheabsenceofnoise,
ˆσ
-0.5 the phenomenon persists: phytoplankton transfer across
0 1 2 3 4 5 6 7 8 9 10 a ∼ λ distance, from a bloom to a no-bloom region,
t (years) c
destabilizes the second one. Thus, if a bloom bubble is
sufficientlylarge(onthescaleofλ ),itwillsurvivediffu-
c
FIG. 7: Top figure: evolution of the mean concentra- sivetransferandgraduallyexpandattheexpensesofthe
tion fields P¯ (thin line) and Z¯ (heavy line). Bottom fig- surroundings. It is possible to see that the situation is
ure: evolution of the normalized RMS fluctuations σˆP = confirmed in the case of a system with just two homoge-
hP˜2i1/2/P¯ (thin line), σˆ = hZ˜2i1/2/Z¯ (heavy line), cˆ =
Z PZ neous compartmens coupled diffusively, one evolving on
sign(hP˜Z˜i)|hP˜Z˜i/(P¯Z¯)|1/2 (dotted line). Same choice of pa-
a bloom cycle, the other on a no-bloom cycle: the no-
rameters as in Figs. 5 and 6.
bloom cycle will always be destabilized and the system
will lock on a global bloom cycle.
This suggests the following picture:
• Demographic noise continuously generates local
fluctuations in P andZ (especially in Z) that may
push some regions of space in a bloom regime.
• If the noise amplitude is sufficiently high, some of
these regions will be large enough for not being
destroyed at once by diffusivity (see passage from
year 1 to year 2 in Fig. 8).
• At this point diffusion, through phytoplankton
transfer, quickly destabilizes the surrounding no-
bloom regions. Noise is not expected to be impor-
tant any more in this phase.
FIG. 8: Evolution due to the effect of diffusion, in the ab- VI. CONCLUSIONS
sence of noise, of domains characterized by bloom dynamics
(yellow in figure; a bloom event at a given pixel is identified
One can imagine a hierarchy in the description of a
by crossing during the year of the threshold P = 12). The
plankton population, that goes, increasing the degree
initial condition at t=0 was a random distribution in space
in refinement, from zero-dimensional models, to models
ofvaluesP,Z inthebloomandno-bloombasinofattraction.
Initial fraction of bloom points: 0.235 (for larger fractions, based on the dynamics of concentration fields, such as
thesystem locks immediately on a bloom-cycle; theopposite P and Z, to individual based models (IBM), taking full
for lower fractions). Domain characteristics and parameters account of demographic stochasticity. In some ways, PZ
as in Figs. 5-6. models(andalike)maybeseenasthemeanfieldapprox-
imation of some IBM. The correspondence between the
two levels of descriptionis established imposing that the
individual reproduction probabilities in the IBM (condi-
of the RMS fluctuations, as depicted in Fig. 7. The tional to the instantaneous local values of the concen-
stronger fluctuation level in the phytoplankton concen- tration fields), have the same form as the corresponding
tration field, at the onset of the bloom events and soon mean rates at the population level.
aftertheirdisappearance,thatcanbeseenincaseB and By definition, a mean field approximation breaks-up
D in Fig. 5, is paralleled by the double peaks in σˆP in when fluctuations are important. Usually, this means
Fig. 7. large amplitude fluctuations. If the mean-field dynamics
The above picture of global destabilization of the no- is close to an instability, however, this is not necessarily
bloom cycle is based on numerical evidence from sim- thecase,andglobalchangesinthesystemcanbeinduced
ulations in a finite domain. One may question whether by small fluctuations.
8
Wehaveprovidedanexampleofsuchasituationinthe equation for ρN can be written in the form
case of a PZ model (the TB model [18]), with destabi-
lization by demographic noise resulting in switching the ∂ρN = exp Ω−1/2 ∂ − ∂
∂t ∂ζ ∂ζ
systemgloballyfromano-bloomtoabloomseasonalcy- Xi nkX=±1 n (cid:16) i+k i(cid:17)o
cle. [Theparameterdeterminingthestrengthofthenoise
is the the dimensionless constant mˆ of Eq. (24)]. From × Wi+k−2Wi ρN. (A1)
Z o
the point of view of the PZ model, this corresponds to a
To derive an equation for ρ≡ρ , we exploit the relation
ζ
decrease of the instability threshold in the seasonaltem-
perature forcing with respect to the mean-field case (see ΩK/2∂ρN = ∂ρ −Ω1/2 Z¯˙ ∂ρ.
Fig. 4). Forshallowwaterconditions(fewmetersdepth), ∂t ∂t i∂ζ
Xi i
and mixing in the water column produced by diffusion
(diffusivity in the range 0.1m2/day), demographic noise Substituting into Eq. (A1) and expanding to O(Ω−1),
appears sufficient to cause switching between regimes. we obtain:
Ttehmepeeffreactturiseoflfutchteuastaimonesmcaognnsiidtuerdeedaisnth[2a9t]o(ffltuhcetugalotiboanl ∂ρ = Ω1/2 Z¯˙i∂ρ
∂t ∂ζ
amplitude∼1oCwithcorrelationtimeequalto35days). Xi i
1 ∂ ∂
As in other systems [25, 26], the modifications of the
+ Ω 1+ −
mean-field dynamics, produced by demographic noise, Xi nnh Ω1/2(cid:16)∂ζi−1 ∂ζi(cid:17)
are expected to be maintained in the case of an infi-
1 ∂ ∂ 2
nitesystems(thermodynamiclimit). Thereasonispartly + − wi−1−wi
trivial: thereproductionratesatthepopulationlevelare 2Ω(cid:16)∂ζi−1 ∂ζi(cid:17) i o
1 ∂ ∂ 1
nonlinear functions of the concentration fields P and Z,
+ 1+ − +
andfluctuationsleadnecessarilytotheirrenormalization. nh Ω1/2(cid:16)∂ζi+1 ∂ζi(cid:17) 2Ω
Suchapicture,however,isincomplete,asthedestabiliza- ∂ ∂ 2
× − w −w ρ, (A2)
tion process, unless the demographic noise is very large, ∂ζ ∂ζ i+1 i
(cid:16) i+1 i(cid:17) i oo
destabilizes the system only locally. [The characteristic
where we have introduced the rate density w =
scale of the fluctuations coincides with that of the Tur- i
Ω−1W = Z κˆ . Equation (A2) could be further sim-
ing patterns of the system, and has nothing to do with i i i
plified exploiting the relation
stochastic demography]. Only later, by an invasive pro-
cess, which appears to be insensitive to system size, the
[wi+1+wi−1−2wi]ρ=0.
destabilized regions inglobate the inactive surroundings,
Xi
leading to a global bloom cycle.
At this point we write w = w¯ +Ω−1/2κˆ ζ , w¯ ≡ κˆ Z¯ ,
i i i i i i i
and expand Eq. (A2) in powers of Ω. We find, to
O(Ω−1/2):
Acknowledgments
[Z¯˙i−(w¯i+1+w¯i−1−2w¯i)]∂ρ =0, (A3)
∂ζ
I wish to thank M. Gatto, R. Casagrandi, A. Lugli`e i
and B. Padedda for interesting and helpful discussion. which gives, taking the continuous limit ∆x → 0, the
This research was funded in part by Regione Autonoma diffusion equation
della Sardegna.
∂Z¯(x,t) ∂2(κ(x)Z¯(x,t))
= . (A4)
∂t ∂x2
Going to next order, we find the equation for the fluctu-
Appendix A: Master equation treatment of diffusion
ations
∂ρ
We consider for simplicity a one-dimensional domain = Π ρ, (A5)
i
∂t
and a single species, say Z. Discretize the domain and Xi
indicate with N the number of individuals in slot i. We
i where
can define the PDF’s:
∂ ∂
Πi = κˆi−1 ζi−1+κˆi+1 ζi+1
ρN(N) ≡ ρN({ΩZ¯i+Ω1/2ζi,i=1,...K}) ∂ζi−1 ∂ζi+1
∂
= Ω−K/2ρζ(ζ). − (κˆi−1ζi−1+κˆi+1ζi+1)
∂ζ
i
1 ∂ ∂ 2
Suppose that individuals are transferreddiffusively from + κˆi−1Z¯i−1 −
slot i to slot i ± 1 with a rate Wi→i+1 = Wi→i−1 ≡ 2h (cid:16)∂ζi−1 ∂ζi(cid:17)
W = N κˆ , where κˆ = κ /(∆x)2, with κ ≡ κ(x ) the ∂ ∂ 2
i i i i i i i + κˆ Z¯ − . (A6)
diffusivity and ∆x the width of the slot. The master i+1 i+1 ∂ζ ∂ζ
(cid:16) i+1 i(cid:17) i
9
From Eqs. (A5-A6) we obtain the equation for the cor- and the remaining terms give the diffusion equation for
relation C (t)=hζ (t)ζ (t)i: C(x,y;t)≡hζ(x,t)ζ(y,t)i:
ij i j
C˙jk = −2(κˆj +κˆk)Cjk+κˆj−1Cj−1,k+κˆk−1Cj,k−1 ∂C(x,y;t) ∂2(κ(x)C(x,y;t))
=
1 ∂t ∂x2
+ κˆ C +κˆ C + π (A7)
j+1 j+1,k k+1 j,k+1 jk
2 ∂2(κ(y)C(x,y;t))
+ . (A8)
where ∂y2
{w¯j−1+2w¯j +w¯j+1}, k =j, Thecorrespondingdynamicsforthe fluctuationfieldζ is
πjk =−{w¯j−1+w¯j}, k =j−1, diffusive as well:
0, k <j−1,
∂ζ(x,t) ∂2(κ(x)ζ(x,t))
 = , (A9)
and π =π . Taking the continuous limit we find that ∂t ∂x2
jk kj
the source term π is infinitesimal:
ij
which leads to the diffusion terms to the RHS of the
π →−4∆xZ¯(x)κ(x)δ′′(x−y)→0, Langevin equations (14).
ij
[1] C.B. Field, M.J. Beherenfeld, J.T. Randerson, and P.G. 48, 591 (1990)
Falkowski, Science 281, 237 (1998) [21] A.M. Edwards, J. Plankton Res. 23, 386 (2001)
[2] P.J.S. Franks, Limnol. Oceanogr. 42, 1297 (1997) [22] W.R. Young, A.J. Roberts and G. Stuhne, Nature 412,
[3] N. Blackburn, T. Fenchel, and J. Mitchell, Science 282, 328 (2001).
2254 (1998) [23] A.J. McKane and T.J. Newman, Phys. Rev. Lett. 94,
[4] A.P.Martin,Philos.Trans.R.Soc.LondonA365,2873 218102 (2005)
(2005) [24] C.R. Doering, K.V. Sargsyan and L.M. Sander, Multi-
[5] R. Reigada, R.M. Hillary, M.A. Bees, J.M. Sancho and scale Modelling and Simulation 3, 283 (2005)
F. Sagues, Proc. R. Soc. Lond.B 270, 875 (2003) [25] T. Butler and D. Reynolds, Phys. Rev. E 79, 032901
[6] A.M. Metcalfe, T.J. Pedley and T.F. Thingstad, J. Ma- (2009)
rineSys. 49, 105 (2004) [26] T. Butler and N. Goldenfeld, Phys. Rev. E 84, 011112
[7] M. L`evy,Lect. Notes Phys.774, 219 (2008) (2011)
[8] A.Bracco,S.ClaytonandC.Pasquero,J.Geophys.Res. [27] I. Siekmann and H. Malchow, Math. Model. Nat. Phe-
114, C02001 (2009) nom. 3, 114 (2008)
[9] W.J. McKiver and Z. Neufeld, Phys. Rev. E 83, 016303 [28] C.S.Holling,Can.Entomol.91,293(1959);C.S.Holling,
(2011) Can. Entomol. 91, 385 (1959)
[10] P.J.S. Franks, Limnol. Oceanogr. 42, 1273 (1997) [29] J.A.Freund,S.Mieruch,B.Scholze,K.WiltshireandU.
[11] U. Siegenhalter and J.L. Sarmiento, Nature 356, 589 Feudel, Ecol. Complex. 3, 129 (2006)
(1993) [30] J.A. Berges, D.E. Varela and P.J. Harrison, Mar. Ecol.
[12] T.J. Smayda,Limnol. Oceanogr. 42, 1137 (1997) Prog. Ser. 225, 139 (2002)
[13] C.S.Yentsch,B.E. Lapointe,N.Poulton andD.A.Phin- [31] E.J. Gonz´alez, T. Matsumura-Tundisi and J.G. Tundisi,
ney,Harmful Algae, 7, 817 (2008) Braz. J. Biol. 68, 69 (2008)
[14] M. Scheffer, S. Rinaldi, Y.A. Kutznetsov and E.H. Van [32] To be precise, in order for the two species to have iden-
Nes, Oikos80, 519 (1997) tical diffusion properties, we could assume that passive
[15] T.J. Smayda,Limnol. Oceanogr. 42, 1132 (1997) transport be dominant over swimming.
[16] D.H.Cushing, Adv.Mar. Biol. 26, 249 (1990) [33] M.Doi,J.Phys.A9,1465(1976);A.S.Mikhailov,Phys.
[17] G.T. Evans and J.S. Parslow, Biol. Oceanogr. 3, 327 Lett.85,214(1981); N.Goldenfeld,J.Phys.A17,2807
(1985) (1985); L. Peliti, PJ. Physique 46, 1469 (1985); H.K
[18] J.E. Truscott and J. Brindley, Bull. Math. Biol. 56, 981 Janssen and U.C. Tauber, Ann.Phys. 315, 147 (2005)
(1994) [34] N.G. Van Kampen, Stochastic processes in physics and
[19] J.H. Steele and E.W. Henderson, Am. Naturalist 117, chemistry (Elsevier, New York,1992)
676 (1981) [35] S.A. Levin and L.A.Segel, Nature259, 659 (1976)
[20] M.Fasham,H.DucklowandS.McKelvie,J.MarineRes.