Rev. Acad. Canar. Cienc., XIV (Núms. 1-2), 121-137 (2002)
Hopf bifurcations and slow-fast cycles in a
model of plankton dynamics
Isabel Fernández, José M. López, José M. Pacheco, César Rodríguez
Universidad de Las Palmas de Gran Canaria,
Departamento de Matemáticas,
Campus de Tafira, 35017 Las Palmas de Gran Canaria,
SPAIN
Tlf: +34 928 458818 ; Fax: +34 928 458811
e-mail: pacheco@dma.ulpgc.es
Abstract
The existence of Hopf bifurcations and slow-fast cycles in the dynamics
of a model for the interaction between phyto- and zooplankton is considered.
A sensible and easily interpretable dimensionless version of the
model is presented, followed by a numerical bifurcation analysis in a twodimensional
parameter space. The biological meaning of the parameters
and qualitative features in phase space are stressed.
Keywords:plankton modela, bifurcation analysis, slow-fast cycles
121
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
1 Introduction
In the open sea, plankton populations live in an ever-changing environ.ment
where equilibrium situations are the exception rather than the rule. Quite
often, sudden and dramatic variations in biomass happen due to the interna!
mechanisms that rule the behaviour of the planktonic components.
Therefore, it is natural to look for models of plankton dynamics able to
account for the onset of oscillations of various classes; among them, the slow-fast
dynamics related to excitable variables is frequently observed in natural
environ.ments. There is no need of formulating complicated models: Even fairly
easy ones can account for very rich dynamics.
In this paper, as in Truscott (1995), plankton is divided into only two classes,
phytoplankton and zooplankton, and they stand in a prey-predator relation-ship.
See also Steele and Henderson (1992). p(t) and z(t) shall represent the
spatially averaged phyto- and zooplankton biomasses, repectively. Both classes
are related to each other through grazing; by which biomass is transferred from
phytoplankton to zooplankton. In the absence of grazing, the time evolution
of phytoplankton is modeled by a logistic law with growth rate r and carrying
capacity k:
p' = rp(l - t)·
Grazing is introduced as a Holling III interaction term:
this means that for large p the consumption rate of phytoplankton per unit of
grazer biomass tends to a saturation value Rma:x· a is the so-called semisat-
122
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
uration constant, which is a measure of how fast the grazing term becomes
asymptotically constant: The closer a is to O, the faster is Rmax approached.
On the other hand, the dynamics of z is a linear one in the absence of phyto-plankton:
z1 = -μz, and the positive grazing effect is obtained by adding the
grazing term in the p dynamics times sorne effectivity coefficient 'Y, a parameter
that can embody the delay effect of the biomass transfer from phytoplankton
into zooplankton. The above considerations are summed up in the ODE model:
'- p p2
p - rp(l - k) - R.naxz a 2 + p2
' p2
z = -μz + 'YR.naxz 2 ...2 • a +y
The aims of this paper are:
a) To present a different dimensionless version of the ODE model where
parameters are easily interpreted from a biological viewpoint.
b) To develop a bifurcation analysis leading to the understanding of the
plankton dynamics through Hopf bifurcations and then to slow-fast cycles. Hopf
bifurcations are one of the most employed tools in describing the onset of oscil-lations
in many systems, ranging from mechanical to biological to social ones.
See Fernández and Pacheco (2001); Pacheco et al. (1997).
e) To obtain further insight of the biological interpretations of the model.
2 Qualitative analysis of the model
2.1 Adimensionalization
Obtaining dimensionless forms for the equations is a primary task before em-barking
on further analyses in any model study. As a rule, equations can be
123
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
written in simpler forms, the number of parameters is reduced, and the remain-ing
ones are found to be easily interpretable combinations of the original ones.
The procedure is not unique and usually units must be carefully chosen. A
most interesting discussion can be found in Fowler (1997), the standard tech-nique
is explained in Edelstein-Keshet (1988), and many interesting examples
are described in Murray (1989).
In this case a nice dimensionless version of the model can be obtained by
scaling the variables
t = t*t, p = p*p, z = z*z
with the new units:
- 1 r
t = -, p = k, Z = D k.
r ~"max
This choice is different from the one proposed in 'l'ruscott (1995), and amounts
to using the linear relaxation time and the carrying capacity of phytoplankton
as new time and phytoplankton biomass units. The zooplankton biomass unit is
simply the fixed fraction of k given by the dimensionless quotient ;,,ax relating
the linear growth rate of phytoplankton and the maximum consumption rate
per unit grazer. By plugging these values into the equations and dropping the
asterisks of the new variables, the model is written in the simpler form:
p' = p(l - p) - z__i__
v2 +p2
z' = /j [ v2 : p2 - w] z.
The original six constants are reduced to only three non-dimensional ones:
Q μ Rmax
V= -k' 8 = -, and /j = --7,
r r
124
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
a μ R.nax
v = -, 8 = -, and f3 =---y,
k r r
a
where v = k expresses the semisaturation constant a as a fraction of the carry-ing
capacity k, 8 = !:!:. is the quotient of the linear growth rates of both species,
r
and /3 = R.nax 'Y is the analogue of 'Y -note that the factor is the inverse fraction
r
of the one used to define the zooplankton unit. Actually, f3 plays no role in the
analyses to follow; instead, the dimensionless combination w = % is introduced
for the sake of simplicity, notational convenience, and interpretability. It is
interesting to observe that w measures how large the linear decay rate μ of zoo-plankton
is with respect to the maximum rate -rR.nax of zooplankton biomass
creation out of phytoplankton. This remark will often be used in the sequel.
2.2 Singular points and their stability
For any initial condition (p0 , zo) in the first orthant of the phase plane, the
corresponding orbit never leaves it, for it follows from the model construction
that both the p-axis and the z-axis are trajectories of the system. The analysis
is therefore restricted to this area.
The manifold z' = O is the union of the p-axis and the line
p=vJ1 -ww '
while the p' = O manifold is formed by the z-axis and the curve, asymptotic to
the z-axis when p-+ O, defined by the equation
(1 - p)(v2 + p2)
z= .
p
The above equations show immediately that two conditions must be fulfilled
for the system to have a singular point interior to the first orthant:
125
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
w<l andp<l
the other only two singular points being {O, O) and {1, O) for any choice of the
parameter values. Note also the equivalence
p<l{=:}v<J~-1.
It is straightforward to show that {O, O) is a saddle point whose stable and
unstable manifolds are the z-axis and the p-axis respectively . The point {1, O)
is also a saddle point with the p-axis as stable manifold -stated simply, this
means that the restricted dynamics is given by a logistic-, while the unstable
manifold is locally given by the curve z = {1 - p)(v2 + p2)/p.
Once the two conditions w < 1 and p < 1 are met, the singular point
in the interior of the first orthant is the unique intersection point of the curves
(Figure 1):
p=vv1-~w
{1 - p)(v2 + p2)
z= .
p
The relationship between v and w governs the stability of this singular point, and
in case there should existan stability shift for sorne pair(s) (w, v), a bifurcation
problem must be considered. This is the aim of the next section.
2.3 A bifurcation analysis
2.3.1 The shape of the p' = O manifold
From a geometric viewpoint the stability of the singular point depends on the
shape of the curved component of the p' = O manifold, the curve
126
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
(1 - p)(112 + p2)
Z= .
p
112 + p2
Its slope ---¡¡.-- + 2(1- p) depends (Figure 2) on the parameter 11. For large
11 the slope is always negative and the curve descends monotonically from +oo
to O along the interval (O, 1]. This is easily proved: Writing
112 + p2 112
--¡¡.--+2(1-p) = 1-2p- p2 <o
2
in the form 2p+ ~ > 1, it follows trivially that it suffices to take 11 large enough
for the inequality to hold. To quantify how large 11 must be, consider that the
minimum m of the expression 2p + 11: occurs when Pm satisfies 2 - 2: = O,
P Pin
2
or equivalently Pm = 11213. Therefore m > 1 if 211213 + ; 13 = 311213 > 1, or
11 > llcrit = ¡¡; = 0.19245 ...
For 11 < 11 crit, two extrema -a minimum and a maximum- appear at the
values Pmin and Pmax1 with O < Pmin < Pmax < l. Therefore the slope of the
curve is positive in the interval (pmin, Pmax) , its graph shows a characteristic
hump form, and the linear stability analysis of the point (p*, z*) depends on
whether p* E (pmin, Pmax) or not. To see when this is the case, write:
.~ Pmin <p =11y~ <Pmax
to obtain the following estimates for w:
that must be satisfied (remember that 11 < llcrit) for the singular point (p*, z*)
to be in the "uphill" part of the hump, to the left of the maximum. When this is
the case, an annular closed region around the singular point can be determined
127
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
in phase space such that any orbit entering it will remain there forever, so the
Poincaré-Bendixson theorem is applied and a limit cycle exists. See Arrowsmith
and Place (1994); Verhulst (1990).
2.3.2 Numerical analyses in the (w, v) plane
The model dynamics depends essentially on the two parameters w and v, and
the bifurcation analysis amounts to study those sets of pairs (w, v) for which
the qualitative aspect of the phase portrait is the same. In this model no simple
expressions for Pmin and Pmax can be obtained for general 11, so the following step
will be the numerical analysis of sorne geometrical shapes in the (w, v) plane.
Remember that the conditions for the existence of a singular point in the
first orthant amount to
w < 1 and v < V~ -1
so the set of feasible pairs in the (w, v) plane is the open set n limited by the
positive w-axis, the positive 11-axis and the graph of v = J ~ -1 (Figure 3).
By establishing a grid over n and computing the eigenvalues ..\(w, v) of the
jacobian matrix at the grid points, the transition from two negative real eigen-values
to negative real part complex ones, and from these to positive real part
complex ones can be tracked across two parabola-like curves that can be de-termined
by fitting adequate trigonometric polynomials to the observed points.
This procedure can be easily programmed in Mathematica or Maple.
Therefore, n is union of three open regions and the two parabola-like curves
separating them1, fN-S and fs-H· fN-S runs from the origin to (1,0) with
1 Here N means "node", S, "spiral", and H, "Hopf".
128
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
a single ma.ximum at (0.4, 0.8); r s-H runs from (0.5, O) to (1, O) and its single
ma.ximum is at (0.75,0.19). The three open regions are (Figure 4):
nN, above rs-H: For (w, 11) in this region, (p*' z*) is a stable node (two neg-ative
real eigenvalues) whose attraction basin is the interior of the first orthant.
ns, limited by both rN-S, rs-H 'and the w-axis: For (w, 11) in this region,
(p*, z*) is a stable spiral point (negative real part of the complex eigenvalues)
and, again, the attraction basin is the interior of the first orthant. If (p*, z*) E
rN-S, it behaves as a degenerate node.
nH, below rs-H: For (w,11) in this region, (p*,z*) is a unstable spiral
point surrounded by a stable limit cycle. Then, rs-H is the bifurcation set
of a transcritical Hopf bifurcation: When crossing it, the stable spiral point
(p*, z*) splits into an unstable spiral point anda stable limit cycle surrounding
it. The attraction basin is the interior of the first orthant as well. Numerical
computations show that the Hopf bifurcation condition (see Arrowsmith and
Place (1994))
[ 8(Re.X(w,11)]
art > 0
rs-n
is also met.
3 Biological interpretation
To obtain a sensible explanation in biological terms of the above discussion the
definition and interpretation of the parameters w and 11 must be remembered:
ó μ a
w = - = --, 11 = -
fJ "YRmax k
129
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
w expresses the relationship between the linear decay rate μ of zooplankton and
the maximum rate 1Rma.x of zooplankton biomass creation out of phytoplankton.
v indicates how fast the maximum grazing rate is approached: Actually,
it is a measure of how much phytoplankton is needed for grazers to reach a
satiated state. Small values indicate rapid satiation, while larger values are the
signal of a harder struggle for existence.
Given that for the existence of (p*, z*) it is required that w < 1, i.e.
μ < 1Rma.x, small values of w mean that the nonlinear "birth" rate /Rma.x
is much larger than the linear decay rate μ. Therefore, if v is small, i.e. the
ma.ximum grazing rate Rma.x is rapidly approached, the net balance is positive
for zooplankton growth. It must be noted that for large v the interval of feasible
w values becomes small quite rapidly: For instance, when v = 1 the interval is
(O, 1/2). Should any w > 1/2 be chosen, the model dynamics shows that the
populations (p, z) -+ (1, O), that is, zooplankton disappears as phytoplankton
settles comfortably at its carrying capacity.
Summing up, the biology underlying the model is characterized by the conflict
between "v, or how much phytoplankton is needed for grazers to reach a
satiated state" and "w, or how efficient is zooplankton in the transformation of
phytoplankton biomass into its own biomass". This is reflected in the following,
purely qualitative table:
v\w small lar ge
small Stable populations Cyclic behaviour
lar ge Stable populations Zooplankton extinction
130
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
It must be observed that the analysis in 2.3.2 shows that stable populations
may be approached either directly (stable nodes) or after decaying oscillations
(stable spiral points).
4 Slow-fast cycles
Roughly speaking, a slow-fast cycle is a cycle along which the travelling speed
of one state variable is dramatically altered on sorne parts of the cycle, and this
variable is called a "fast" variable. The other state variable is the "slow" one.
The graph of the time evolution of the fast variable shows a typical oscillatory
pattern: A wave train of steep humps with more slowly decreasing lees (Figure
5).
Inspection of the relative position of the p' = O and z' = O manifolds reveals
a slow-fast cycle situation if one of them -the one corresponding to the fast
variable- is an S-shaped curve and the other one crosses it transversally through
sorne point in the central part of the S. This is exactly the case of the model for
(w, v) E nH: To see it, just rotate the phase plane by an angle of rr/2 around
the origin to "discover" the S-shapedness of the p' = O manifold and how the
straight line z' = O goes across it. In this case p is the fast variable.
A slow-fast cycle will typically show (p, z) travelling counterclockwise along
p' = O (the curve z = (1 - p)(v2 + p2 )/p) up to CPmax, Zmax), then jumping
horizontally to the descending part of the same curve, then descending to the
point CPmin, Zmin), where a new horizontal track is traversed, and so on (Figure
6). See Verhulst (1990).
According to the Hopf bifurcation theorem, the amplitude of the limit cycle
created when an stable spiral points splits into an unstable one and a stable
131
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
cycle grows as the pair of bifurcation parameters ( w, v) crosses transversally the
bifurcation set rs-H into the region nH. Therefore, any sensible analysis will
track the parameter values (w, v) E D.H in such a way that the gradient of the
amplitude is maximum. The numerical experiments show that this is the case
for the right comer of D.H; this means that v must be small -grazers are soon
satiated- and w rather large -the nonlinear rate '"Yllrna:x is smaller than the linear
decay rate μ.
5 Conclusions and views
The model considered in this paper stresses the role of various relationships
when interpreted from a biological viewpoint and accounts for severa! interesting
features of plankton dynamics.
v --or its dimensional counterpart a- controls how fast predators, i. e. zooplankton,
are satiated: Small values of the parameter correspond with fast
satiation. From the Oceanography viewpoint, a small v can correspond to a
stage where phytoplankton biomass grows quickly: a so-called "algal bloom".
The other interesting parameter is w, it is a measure of the relative importance
of the linear decay rate of zooplankton and the maximum nonlinear rate of
zooplankton biomass creation out of phytoplankton. Both parameters span a
two-dimensional space where a feasible subset n of parameter pairs is obtained.
Moreover, there exists a bifurcation set -the curve fs-H - separating stable
spiral behaviour from unstable spiral behaviour and a stable limit cycle.
Therefore the Hopf bifurcation is a characteristic feature of this model meaning
that for small values of v and rather large values of w plankton can develop
stable oscillating patterns out of any initial distribution of phyto- and zoo-
132
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
plankton. H vis very small and w ~ 1, the stable cycle shows a typical slow-fast
behaviour. Informally, this means that when predators are rapidly satiated and
their efficiency in transforming prey into predator biomass is not large enough
to balance the linear decay, sudden and dramatic changes can occur to both
components of the planktonic system. These results agree with the general
conclusions of Truscott (1995), exception made of the Hopf bifurcation stage.
An interesting point for further study is the introduction of periodic behaviour
in sorne of the model parameters. Periodic upwellings controlled by
variations in the depth of the thermocline -as in the EN phase of ENSO- can
be represented by a modulation like v = v(t) with v(t) = v(t + T), where
the "only" difficulty is the determination of the period T. See Fernández and
Pacheco (2000) and references therein. In addition to that, upwellings usually
imply the appearance of predators externa! to the planktonic systems -fish- that
can highly modify the dynamics and due to their mobility make it necessary, as
in Malchow (1994), to consider spatial distribution effects as well.
References
[1] Arrowsmith, D. and Place, C.,1994. An introduction to dynamical systems,
Cambridge UP, Cambridge, 423 pp.
[2] Edelstein-Keshet, L., 1988. Mathematical models in biology. Random
House-Birkhauser, New York, 586 pp.
[3] Fernández, l. and Pacheco, J., 2000. Bases para la predicción de ENSO.
In: R.García and E. Hernández (Editors), El Niño: Climatología, efectos y
predicción, Universidad Complutense-Mapfre, Madrid, pp. 93-132.
133
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
[4] Fernández. I. and Pacheco, J., 2001. Sorne mathematical aspects in the
modelling of environmental quality. In: M. Matthies (Editor), Integrative
approaches to natural and social dynamics, Springer Verlag, Berlin, pp.
235-248.
[5] Fowler, A.,1997. Mathematical models in the applied sciences. Cambridge
Univ. Press, Cambridge, 402 pp.
[6] Malchow, H., 1994. Non-equilibrium structures in plankton dynamics, Ecol.
Modelling, 75: 123-134.
[7] Murray, J.,1989. Mathematical biology. Springer Verlag, New York, 767 pp.
[8] Pacheco, J., Rodríguez, C. and Fernández, I.,1997. Hopf bifurcations in a
predator-prey model with social predator behaviour, Ecol. Modelling, 105:
83-87.
[9] Steele, J. and Henderson, W., 1992. The role of predation in plankton
models, J. Plankton Res., 14(1): 157-172.
[10] Truscott, J. E.,1995. Environmental forcing of simple plankton models, J.
Plankton Res., 17(12): 2207-2232.
[11] Verhulst, F.,1990. Nonlinear differential equations and dynamical systems.
Springer Verlag, Berlin, 277 pp.
134
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
FIGURES
0.3
0.2
o.
0,2 0.4 0.6 0.8 1 i;hyto
Figure 1: The manifolds p' = O and z' = O. The curve has the equatipn
Z = (1-p)(v'+p'), with 11 = 0.1
p
0.8
o.a 1 ~
Figure 2: The curve in Fig. 1 and the graph of ita derivative. Note the
characteristic hump due to the changing sign of the derivative.
135
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
y
3
2.5
.....
1.5
1 Q
0.5
w 0.2 o.a 1
Figure 3: The feasible set !1 for the existence of a singular point in the first
orthant. The limiting curve has the equation v = J t - 1
y
1.5
1
w
0.2 0.4 0.6 o.a 1
Figure 4: The feasible region in Fig. 3, divided into three components.
The stability type of the singular point (p'", z'") is the same for (w, v) in each
region. The curve rs-H is the bifurcation set for Hopf bifurcations.
136
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
Phase Plana
0,6
0,48
0,36
oo
N
0,24
0,12
o
o 0,2 0,4 0,6 0,8 1
Phyto
Figure 5: Phase plane for v = 0.05 and w = 0.95. Shown are the atable
cycle and the unstable spiral point (1), with coordinates (0.22, 0.18)
Phytoplankton vs. time
0,8
0,6
o
>.
.i= a..
0,4
0,2
o
o 40 80 120 160 200
time
Figure 6:The slow-fast cycle of phytoplankton for v = 0.05 and w = 0.95.
Note the typical shape of the wavetrain.
137
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017