Mixed-mode oscillations and chaos in aprey-predator system with dormancy of
predators
Masataka Kuwamura∗and Hayato Chiba†

Abstract
It is shown that the dormancy of predators induces mixed-mode
oscillations and chaos in the population dynamics of a prey-predator
system under certain conditions. The mixed-mode oscillations andchaos are shown to bifurcate from a coexisting equilibrium by meansof the theory of fast-slow systems. These results may help to find experimental conditions under which one can demonstrate chaotic population dynamics in a simple phytoplankton-zooplankton(-resting eggs)
community in a microcosm with a short duration.
keywords: mixed-mode oscillations, chaos, dormancy, prey-predator,
fast-slow systemrunning head: mixed-mode oscillations and chaoscorresponding author: Masataka KUWAMURA∗

Graduate School of Human Development and Environment, Kobe University, Kobe
657-8501, Japan, e-mail: kuwamura@main.h.kobe-u.ac.jp†
Graduate School of Informatics, Kyoto University, Kyoto 606-8501, Japan. Current
address: Faculty of Mathematics, Kyushu University, Fukuoka 819-0395, Japan, e-mail:chiba@math.kyushu-u.ac.jp1

It is well known that zooplankton mainly produce subitaneous
eggs in comfortable environments; however, they produce fertilized eggs (resting eggs, dormant state) to survive periods of harshliving conditions. In this study, we investigate the eﬀects of dormancy on the population dynamics of phytoplankton-zooplanktonin freshwater ecosystems through a simple prey-predator modelthat considers the case in which predators enter dormancy in theevent of a shortage of prey. The result is that the dormancy of
predators can induce mixed-mode oscillations and chaos in the population dynamics of a prey-predator system under certain conditions. This suggests that the population dynamics of zooplankton
species with a long life and a high foraging ability may exhibit
complex behavior such as mixed-mode oscillations and chaos in anenvironment where food deﬁciency occurs for an extended duration. In this case, the period of mixed-mode oscillations (chaos)
has the same order as the average dormancy period. These pieces
of information may help to ﬁnd experimental conditions underwhich one can demonstrate chaotic population dynamics in a simple phytoplankton-zooplankton(-resting eggs) community in a microcosm with a short duration. In addition, a mechanism of theonset of the mixed-mode oscillations and chaos is revealed by meansof the theory of fast-slow systems.

2

1

Introduction

Many short-lived organisms have ways of coping with variable environmental conditions that signiﬁcantly impact individual growth, reproduction, andsurvivorship. In particular, zooplankton mainly produce subitaneous eggs
in comfortable environments; however, they produce fertilized eggs (resting
eggs, dormant state) to survive periods of harsh living conditions or catastrophic events. The ability to produce resting eggs is most probably an earlyfeature in zooplankton evolution. In addition, the production of resting eggs
gives rise to a biotic reservoir, analogous to plant seed banks that preserve
genetic material over long periods. Therefore, dormancy has been a central
subject in the study of zooplankton, as seen in the review paper by Gyllströmand Hansson [15].
In this study, we focus on the eﬀect of dormancy on the population dynamics of phytoplankton-zooplankton in freshwater ecosystems. With regard
to this issue, McCauley et al. [27] found that dormancy has a crucial eﬀect
on the population dynamics of Daphnia and its algal prey; that is, the amplitude of the prey-predator cycles of Daphnia and its algal prey in microcosms
increases when a portion of females producing resting eggs is replaced byasexually reproducing gravid females. This suggests that the dormancy of
predators can stabilize the population dynamics of prey-predator systems.
In order to explain these experimental results, a simple prey-predator model
was proposed in [22] according to the natural strategy of a predator that
produces resting eggs in the event of a shortage of prey. The model is athree-component autonomous ordinary diﬀerential equation as follows:

p
dp


= r(1 − )p − f (p)z


dt
K



dz
(1.1)
= k1 µ(p)f (p)z + αw − d1 z

dt





 dw = k2 (1 − µ(p))f (p)z − αw − d2 w,
dt
where p and z denote the population densities of prey and predators, respectively, and w denotes the population density of predators with a dormancy
state (resting eggs). r and K denote the intrinsic growth rate and the carrying capacity of prey, respectively. The function f (p) represents the Hollingtype II functional response deﬁned by(1.2)

f (p) =

bp,
c+p

where b and c denote the maximum foraging rate and the half saturation
constant, respectively. k1 and k2 denote the increasing rates of predatorsin the active and dormant states, respectively. A switching function µ(p)3deﬁned as the sigmoid function(1.3)µ(p) =

p − η 
1
1
tanh
+1
 −2(p−η)  =2σ
1 + expσ

controls the induction of dormancy, where η and σ denote the switchinglevel and the sharpness of the switching eﬀect, respectively. This function
implies that predators produce more resting eggs than subitaneous eggs when
the prey density decreases below a certain level η. d1 and d2 denote themortality rates of the active and dormant predators, respectively. α denotes
the hatching rate, i.e., resting eggs have a dormancy period with 1/α on
average. It should be noted that (1.1) is reduced to a classical prey-predator
system known as the MacArthur-Rosenzweig model [33]:

(1.4)

 dp
p


 dt = r(1 − K )p − f (p)zdz

= k1 f (p)z − d1 z,dtwhen α = 0 and µ(p) ≡ 1. It is well known that as the carrying capacity K
increases, the population dynamics of (1.4) is destabilized. In other words, a
coexisting equilibrium of (1.4) becomes unstable and a stable periodic orbit(prey-predator limit cycle) appears through the super-critical Hopf bifurcation.
In this study, we show that (1.1) exhibits complex dynamics by means of
the theory of fast-slow systems. It is found that the dormancy of predators
can induce mixed-mode oscillations and chaos in the population dynamics
of the prey-predator system (1.1) under certain conditions. This suggeststhat the population dynamics of zooplankton species with a long life anda high foraging ability may exhibit complex behavior such as mixed-mode
oscillations and chaos in an environment where food deﬁciency occurs foran extended duration. In this case, the period of mixed-mode oscillations(chaos) has the same order as the average dormancy period. These pieces
of information may help to ﬁnd experimental conditions under which one
can demonstrate chaotic population dynamics in a simple phytoplanktonzooplankton community in a microcosm with a short duration (day). On theother hand, in a complex plankton community consisting of many species ina large area with a long duration (year), the existence of chaotic population
dynamics has been reported in [5].
The remainder of this paper is organized as follows. In the next section,
we show numerical results of (1.1); this set of equations exhibits mixed-mode
oscillations and chaos for certain parameter values. In section 3, we introduce
a fast-slow system related to (1.1) to understand the mechanism of the onsetof mixed-mode oscillations and chaos in (1.1) in the framework of the bifurcation theory [23] and the geometric singular perturbation theory [19, 20].
4

Table 1: Deﬁnitions of variables and parameters and default parameter values
used in simulations.
Deﬁnitions
Variablespzw
Parameters
r
K
b
c
d1
d2
k1
k2
α
σ
η

ValuesUnits
mg L−1
mg L−1
mg L−1phytoplankton density
zooplankton density
resting-egg density
intrinsic growth rate of phytoplakton
carrying capacity of phytoplakton
maximum foraging rate
half-saturation constantmortality rate of zooplanktonmortality rate of resting eggs
increasing rate of zooplankton
increasing rate of resting eggs
hatching rate
sharpness of switching eﬀectlevel of switching eﬀect

0.5
day−1
variablemg L−1
2.0
day−1
2.0
mg L−10.2
day−10.0001day−1
0.6 dimensionless
0.12 dimensionless
0.02
day−1
0.1 dimensionless
1.0
mg L−1

It is shown that a mixed-mode oscillation bifurcates from a coexisting equilibrium in the fast-slow system and it induces chaotic dynamics according to
the theory described in [7]. Section 4 presents the concluding remarks.

2

Numerical results: mixed-mode oscillations
and chaosIn this section, we numerically investigate the bifurcation structure of (1.1)
with respect to K, and show that (1.1) exhibits mixed-mode oscillations andchaos under certain conditions.

2.1

Parametrization

Based on our previous work [22], we choose default parameter values for thenumerical simulations in this paper. As seen in [13, 14, 34, 36], these values,
listed in Table 1, are not inconsistent with experimental results. They can beconsidered as a reference to study the qualitative properties of the population
dynamics of phytoplankton-zooplankton communities under the model (1.1),
although they do not correspond to one identical species.
5A few remarks on some parameters are in order: k1 and k2 are dimensionless constants satisfying 0 < k1 < 1 and 0 < k2 < 1, respectively. They
are deﬁned by the eﬃciency at which a predator converts food into reproduction energy divided by the costs of producing subitaneous and resting
eggs, respectively. Because the ratio of the costs for producing subitaneousand resting eggs is estimated to be 1 : 5 by [15], we set k2 = k1 /5. Because the mortality rate of resting eggs d2 is suggested to be very smallin [17], we choose d2 = 0.0001. This implies that the decrease in resting eggs is essentially caused by the hatching of resting eggs, i.e., d2  α.
Hence, the death of resting eggs does not aﬀect the dynamics of (1.1) because−αw − d2 w = −α(1 + d2 /α)w ≈ −αw holds in the right-hand side of thethird equation of (1.1).
Throughout this paper, we use the parameter values and units listed inTable 1, and do not mention these again unless any ambiguity occurs.

2.2

Bifurcation analysis

In this subsection, we study the bifurcation structure of (1.1) with respect
to K, and consider conditions under which (1.1) exhibits mixed-mode oscillations and chaos. We observe how the bifurcation structure changes if
the maximum foraging rate b and the mortality rate of active predators d1
vary. In addition, we numerically estimate the ratio between the period of
mixed-mode oscillations (chaos) and the average dormancy period given by
1/α. The parameters b, d1 , and α have the same dimension (day−1 ) and theyplay important roles in directly controlling the dynamics of active and/or
dormant predators.
First, we consider the non-negative equilibria of (1.1). It is easy to verify
that (1.1) has two equilibria (K, 0, 0) and (p∗ , z ∗ , w∗ ) except for the trivial
one (0, 0, 0), where p∗ , z ∗ , and w∗ are deﬁned byk1 µ(p∗ )f (p∗ ) +

(2.1)

k2 α
(1 − µ(p∗ ))f (p∗ ) = d1 ,α + d2

and
(2.2)

z∗ =

r(1 − p∗ /K)p∗,f (p∗ )w∗ =k2
(1 − µ(p∗ ))f (p∗ )z ∗ .
α + d2It has been shown in [22] that (p∗ , z ∗ , w∗ ) transcritically bifurcates from(K, 0, 0) as K increases from a suﬃciently small positive value. That is,(K, 0, 0) is stable for 0 < K < K tc whereas it is unstable for K > K tc , where
K tc is a transcritical bifurcation point. Such transcritical bifurcation is also
observed in the MacArthur-Rosenzweig model (1.4). Moreover, (p∗ , z ∗ , w∗ )is stable for K tc < K < K H whereas it is unstable for K > K H , where K H
is the Hopf bifurcation point. In this case, we call (p∗ , z ∗ , w∗ ) a coexisting
equilibrium of (1.1).
6

By using AUTO [11], a software package used for studying the bifurcation
structure of diﬀerential equations, we can capture the bifurcation structureof (1.1) with respect to K under the parameter values listed in Table 1.
In this case, the coexisting equilibrium (p∗ , z ∗ , w∗ ) bifurcating from (K, 0, 0)
becomes unstable as K increases and a periodic orbit appears through the
sub-critical Hopf bifurcation at K H ≈ 5.76, the bifurcation diagram for whichis shown in Fig.1 (a). Moreover, almost every solution of (1.1) converges toa stable periodic orbit on the bifurcating branch for K > K H (see [22] fordetails).
Next, we focus on the maximum foraging rate b, and study how thebifurcation structure with respect to K and the dynamics of (1.1) change
when we increase b from b = 2.0 and other parameter values are ﬁxed as listed
in Table 1. Fig.1 (b) shows a numerical result of the bifurcation diagram of(1.1) with respect to K for b = 7.0. In this case, a stable periodic orbit
bifurcates from the coexisting equilibrium through the Hopf bifurcation at
K H ≈ 2.98. However, the periodic orbit becomes unstable through the period
doubling bifurcation at K D ≈ 4.05. In fact, AUTO detects the multiplier−1 of the linearization of the Poincaré map around the bifurcating periodic
orbit, and numerical simulations show that almost every solution converges
to the periodic solution shown in Fig.2.
As K increases considerably, numerical simulations show that the attractor of (1.1) for b = 7.0 exhibits more complex dynamics. For example,
when K = 5.0, an attracting mixed-mode oscillatory solution appears asshown in Fig.3, and almost every solution of (1.1) converges to it. As K
increases further, the mixed-mode oscillation undergoes a succession of bifurcations and a chaotic attractor appears. Fig.4 shows a chaotic solution
for K = 6.0. In fact, the maximum Lyapunov exponent is numerically obtained as λ ≈ 0.0042 > 0. Similar numerical results can be obtained when
we decrease the mortality rate of zooplankton d1 from d1 = 0.2 and otherparameter values are ﬁxed as listed in Table 1. In fact, when d1 = 0.05, we
can obtain a mixed-mode oscillation for K = 2.7 and chaos for K = 3.0.
These mixed-mode oscillations and chaos cannot be observed for suﬃciently small α because the dynamics of (1.1) is almost decoupled to thoseof the (p, z)- and w-components. Similarly, we cannot observe mixed-mode
oscillations and chaos for suﬃciently large α because w decays exponentiallyfast; hence, the asymptotic dynamics of (1.1) can be reduced to the dynamics of the (p, z)-component. Thus, (1.1) has mixed-mode oscillations andchaos only for intermediate (slightly small) values of α under conditions thatthe maximum foraging rate b is large and/or the mortality rate of active
predators d1 is small (as compared to the reference values listed in Table 1).
From numerical simulations for various parameter values, the period ofa mixed-mode oscillation T is roughly estimated as 1 < T /s < 10, wheres = 1/α is the average dormancy period; that is, T has the same order as
the average dormancy period.
7

size

size

(b)

(a)

Fig. 1: The bifurcation diagrams of (1.1) with respect to K for (a) b = 2.0
and (b) b = 7.0. Other parameter values are as listed in Table 1. Here,the horizontal and vertical axes indicate K and the size (in the L2 -sense)of solutions of (1.1), respectively. The solid line indicates asymptoticallystable solutions, whereas the dashed line indicates unstable ones. The whiteand black squares indicate the transcritical and the Hopf bifurcation points,
respectively. The circle indicates the period doubling bifurcation point.

Fig. 2: An attracting solution of (1.1) for K = 4.05 when b = 7.0. Other
parameter values are as listed in Table 1. The solution orbit in the phasespace (p, z, w) is represented.

8

Fig. 3: An attracting mixed-mode oscillation of (1.1) for K = 5.0 whenb = 7.0. Other parameter values are as listed in Table 1. The graphs of p, z,
and w in t (on the left) and the solution orbit in the phase space (p, z, w)(on the right) are represented.

It should be noted that the prey density satisﬁes p < η in Figs.3 and 4,
where η is the level of switching of the dormancy. Numerical simulations forvarious parameter values show that p(t) ≤ η holds for a rather wide rangeof values of t ( p(t) > η holds for some t under certain parameter values).
This suggests that mixed-mode oscillations and chaos of the population of
zooplankton can be observed in an environment where food deﬁciency occurs
for an extended duration.
Figs.3 and 4 show that the derivatives dw/dt along the orbits are smaller
than dp/dt and dz/dt on average. We expect that the onset of mixed-mode
oscillations and chaos may be well explained by using the fast-slow system(3.1) described below. Indeed, the shapes of the orbits shown in Figs.3 and
4 are typical for a three-dimensional fast-slow system, as reported in [7].

3

An approach to mixed-mode oscillations andchaos via a fast-slow system

In this section, we investigate a possible mathematical mechanism to generate
the mixed-mode oscillations and chaos of (1.1). Our approach is based onstandard arguments of the bifurcation theory and the geometric singular
perturbation theory with the aid of numerical computations [9, 12, 19, 20,
9Fig. 4: A chaotic attractor of (1.1) for K = 6.0 when b = 7.0. Other
parameter values are as listed in Table 1. The graphs of p, z, and w in t (onthe left) and the solution orbit in the phase space (p, z, w) (on the right) arerepresented. The maximum Lyapunov exponent is numerically obtained as
λ ≈ 0.0042 > 0.

10

23, 24, 29].
In Sec.3.1, we introduce a fast-slow system related to (1.1) and investigate
its bifurcation diagram with respect to K. In this system, a small parameter
ε is introduced to apply a singular perturbation method. The mathematicalmethod has a limitation in that the value of ε must be suﬃciently small,
although ε = 1.0 in the original system (1.1). However, it is often thecase that the dynamic structure found by the singular perturbation method
qualitatively persists even for an ε value that is not so small. This appears to
be true in our case as well as demonstrated by the numerical simulations in
the previous section. We further conﬁrm this point by numerically showingthat solution orbits (Fig.4) on the chaotic attractor for (1.1) move along
the critical manifold (introduced in Sec.3.2) of the fast-slow system (Fig.9).
Although the fast-slow system (3.1) does not have a direct physical relevance,
it is a useful mathematical procedure for investigating the dynamical behaviorof the original system (1.1). In Sec.3.3, the bifurcation diagram of the fastslow system with respect to b is investigated and it is shown that a mixedmode oscillation bifurcates from a coexisting equilibrium. In Sec.3.4, thebifurcation diagram of the fast-slow system with respect to ε is investigatedand a chaotic attractor is shown to bifurcate from the mixed-mode oscillation.
Throughout this section, we perform numerical computations by using theparameter values related to the chaotic attractor shown in Fig.4. However,
the results in this section are applicable to understand the chaos of (1.1) for
other parameter values.

3.1

Fast-slow system

Motivated by the argument presented at the end of Section 2, we introduce
a small parameter ε > 0 and deﬁne the fast-slow system related to (1.1) to
be

p
dp


= r(1 − )p − f (p)z


dt
K



dz
(3.1)
= k1 µ(p)f (p)z + αw − d1 z

dt





 dw = ε(k2 (1 − µ(p))f (p)z − αw − d2 w).
dt
Our purpose is to study the bifurcation structure of (3.1) with respect to K
and b for small ε.
First, we set b = 7.0, ε = 0.2, and other parameter values as listed inTable 1, and investigate the bifurcation diagram of (3.1) with respect to K
by using AUTO to compare it with that of (1.1) obtained as shown in Fig.1(b). The result is shown in Fig.5. The left-hand side of Fig.5 shows thatthe bifurcation diagram of (3.1) is similar to that of (1.1) (see Fig.1 (b))
when K is smaller than the Hopf bifurcation point (K < K H ). However,
11

size

1.70e+00

1.60e+00

2

11.50e+00

3
1.40e+00

4
1.30e+00
0.00e+00

2.00e-01

4.00e-01

6.00e-01

8.00e-01

1.00e+00

Fig. 5: The bifurcation diagram of (3.1) with respect to K (on the left) andthe phase portrait (on the right) for ε = 0.2 and b = 7.0, where b = 7.0 is
the same as that in Fig.4. Other parameter values are as listed in Table 1.
Here, the horizontal and vertical axes indicate K and the size (in the L2 sense) of solutions of (3.1), respectively. The solid line in the left-hand side
indicates asymptotically stable solutions, whereas the dashed line indicatesunstable ones. The white and black squares indicate the transcritical andthe Hopf bifurcation points, respectively. Periodic orbits corresponding to
the numbered points in the bifurcation diagram are respectively drawn in
the projected phase plain (p, w) on the right-hand side.

it is remarkable that as K increases, a periodic orbit grows and becomes amixed-mode oscillation along the stable branch bifurcating from a coexisting
equilibrium, as shown in the right-hand side of Fig.5. This suggests thatthe fast-slow system (3.1) is suitable for studying the bifurcation structureof mixed-mode oscillations of (1.1). In the following subsections, we clarify
the mechanism of the onset of mixed-mode oscillations and chaos using thetheory described in [7].

3.2

Critical manifold

It is well known [7, 19, 20, 29] that the dynamics of the fast-slow system (3.1)
can be explained by the critical manifold of (3.1) and the dynamics on it asfollows.

12

saddle-focus

transcritical

w

2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

unstable focus N1
stable focus N2
stable node N3
saddle
N4

Hopf

0.25
0.2
0.15
0.1
0.05
0

z

-0.2

0

0.2

0.4

p

0.6

0.8

1

Fig. 6: The critical manifold of (3.1) for b = 7.0 and K = 6.0. These values
are the same as those in Fig.4, and other parameter values are as listed inTable 1.

Setting ε = 0 in (3.1), we have the unperturbed system of (3.1):

p
dp


= r(1 − )p − f (p)z


dt
K



dz
(3.2)
= k1 µ(p)f (p)z + αw − d1 z

dt





 dw = 0 .
dt
The critical manifold of (3.1) is deﬁned by the set of ﬁxed points of (3.2) as(3.3)
p
M = { (p, z, w) ∈ R3 | r(1 − )p − f (p)z = 0, k1 µ(p)f (p)z + αw − d1 z = 0 }.
K
We consider M in the region p ≥ 0, z ≥ 0, and w ≥ 0 because p, z,
and w denote population densities. In this region, the critical manifold M
for b = 7.0 and K = 6.0 is numerically calculated and the result is shown in
Fig.6. It is found that M consists of the line M1 = {p = 0, αw−d1 z = 0} andthe parabola-like curve M2 . When regarding w as a bifurcation parameter,
the red points on M2 (branch N1 ) in Fig.6 indicate unstable focuses; thatis, the linearized matrix of the right-hand side of (3.2) corresponding toeach red point on M2 has two complex eigenvalues with positive real parts.
The green branch on M2 (branch N2 ) consists of stable focuses and it isconnected to the red branch at the Hopf bifurcation point. On the other
hand, the blue and purple branches on M2 (branches N3 and N4 , respectively)
consist of stable nodes and saddle points, respectively. The purple branchconnects to the green branch at the fold point of M2 , which is the saddle-focus
bifurcation point. Similarly, the purple and blue branches on M1 consist of13

saddle points and stable nodes, respectively. The intersection of M1 and M2
is the transcritical bifurcation point. It should be noted that the criticalmanifold of (3.1) has the same structure as that shown in Fig.6 for a rather
wide range of values of b and K.
The dynamics of w of (3.1) on the critical manifold M is given by theslow system
(3.4)

dw
= ε(k2 (1 − µ(p))f (p)z − αw − d2 w)
dt(p,z,w)∈M

that is, the dynamics of (3.1) on M is reduced to the dynamics deﬁned by
the one-dimensional equation (3.4).

3.3

Bifurcation from a coexisting equilibrium to a mixedmode oscillation

In this subsection, we investigate the behavior of a coexisting equilibriumon the critical manifold M and study the bifurcation structure of (3.1) withrespect to b in order to reveal the onset of mixed-mode oscillations. It will beshown that a mixed-mode oscillation bifurcates from a coexisting equilibrium
through the Hopf bifurcation at the fold point of M (cf. [19]). Throughout
this subsection, we set ε = 0.2, K = 6.0, and we increase b from b = 2.0 to
b = 7.0, where b = 2.0 is the reference value listed in Table 1 and b = 7.0 is
the same as that in Fig.4.
Setting b = 2.0, we have an unstable equilibrium (coexisting equilibrium)
near the red branch of M2 (branch N1 ) and a stable periodic orbit (preypredator cycle) of (3.1) in a region p > 0, z > 0, and w > 0 (see Fig.7).
As we increase b from b = 2.0, the coexisting equilibrium moves upward
along M2 and into the green branch on M2 (the branch N2 ). This implies
that the coexisting equilibrium becomes stable, and an unstable periodic
orbit appears through the Hopf bifurcation at b ≈ 2.16 (See Fig.8).
As b increases considerably, the stable coexisting equilibrium continues tomove upward and approaches the fold point of M2 . On the other hand, thestable periodic orbit (prey-predator cycle) collides with the unstable periodic
orbit and disappears through the saddle-node bifurcation at b ≈ 3.14 (seeFig.8).
As we further increase b, the coexisting equilibrium becomes unstable
near the fold point of M2 , and a stable periodic orbit appears through the
Hopf bifurcation at b ≈ 6.46. The bifurcating periodic orbit grows and turns
into a mixed-mode oscillation as b increases (see Figs.8 and 9 (a)).
Thus, we see that the mixed-mode oscillations of (3.1) bifurcate from a
coexisting equilibrium when we increase b from the reference value b = 2.0
listed in Table 1. Moreover, it is found that (3.1) exhibits mixed-mode oscillations when we choose parameter values such that a coexisting equilibrium
is near the fold point of the critical manifold. For example, taking large
14

w

unstable focus N1
stable focus N2
stable node N3
saddle
N4

2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

2.5
2
1.5
1
0.5
0

z

-0.5

0

0.5

1

1.5

p

2

2.5

3

3.5

Fig. 7: The critical manifold, a stable periodic orbit, and an unstable equilibrium of (3.1) for K = 6.0, ε = 0.2, and b = 2.0, where K = 6.0 is thesame as that in Fig.4 and b = 2.0 is the reference value of b listed in Table 1.
Other parameter values are as listed in Table 1. The black curve and crosspoint denote a periodic orbit and an equilibrium, respectively.

size
1.90e+00

6.5
1.80e+00

6.6
1.70e+00

1.60e+00

6.7

1.50e+00

1.40e+00

7.0
1.30e+00
0.00e+00

2.00e-01

4.00e-01

6.00e-01

8.00e-01

1.00e+00

Fig. 8: The bifurcation diagram of (3.1) with respect to b (on the left) andphase portrait (on the right) for K = 6.0 and ε = 0.2, where K = 6.0 is
the same as that in Fig.4. Other parameter values are as listed in Table
1. Here, the horizontal and vertical axes indicate b and the size (in the L2 sense) of solutions of (3.1), respectively. The solid line in the left-hand side
indicates asymptotically stable solutions, whereas the dashed line indicatesunstable ones. The black square indicates the Hopf bifurcation point. Periodic orbits on the stable branch bifurcating from a coexisting equilibrium at
b = 6.5, 6.6, 6.7, and 7.0 are respectively drawn in the projected phase plain(p, w) on the right, where b = 7.0 is the same as that in Fig.4.

15

unstable focus N1
stable focus N2
stable node N3
saddleN4

1.9

unstable focus N1
stable focus N2
stable node N3
saddle
N4

1.9

1.8

1.8

1.7

1.7

1.6

1.6

w 1.5

w 1.5

1.4

1.4

1.3

1.3
1.2

1.2
0.25
0.2
0.15
0.1
0.05
0

0.25
0.2
0.15
0.1
0.05

z

-0.2

0

0.2

0.4

p

0.6

0.8

1

z

-0.2

0

0.2

0.4

p

0.6

0.8

1

0

Fig. 9: The critical manifold and an attractor of (3.1) when b = 7.0 and
K = 6.0; (a) a mixed-mode oscillation for ε = 0.2 and (b) a chaotic attractor
for ε = 1.0. Other parameter values are as listed in Table 1.

Fig. 10: A schematic view of the Poincaré section Σ.

values of b (as compared to the reference value listed in Table 1) is one ofconditions under which a coexisting equilibrium is near the fold point of thecritical manifold.

3.4

Bifurcation from the mixed-mode oscillation to chaos

According to the parameter values related to the chaotic attractor shownin Fig.4, we set b = 7.0, K = 6.0, ε = 0.2, and other parameter values aslisted in Table 1. Then, we obtain the critical manifold and a mixed-mode
oscillation of (3.1), as shown in Fig.9 (a). Next, when we increase ε from 0.2
to 1.0, the mixed-mode oscillation becomes chaotic (Fig.9 (b)) through thecascade bifurcation, as seen in Section 2. Based on the theory described in[7], the mechanism of the onset of mixed-mode oscillations and chaos can be
explained as follows.
Let us deﬁne a Poincaré section Σ that transversally intersects an attrac16

(a)

(b)

R
Fig. 11: Positional relationship of the rectangle R and Π(R).

tor (the mixed-mode oscillation or the chaotic attractor) of (3.1) denoted by
A, as shown in Fig.10. Take a small open set U ⊂ Σ including Σ ∩ A andconsider a solution of (3.1) starting from y ∈ U . Since the green branch on
M2 (branch N2 ) consists of stable focuses of (3.2), the orbit approaches thisbranch. Since the right-hand side of the slow dynamics (3.4) is positive on
the green branch on M2 , the orbit moves upward rotating around M2 andreaches the fold point of M2 . Then, the orbit leaves the fold point and approaches the blue branch on M1 (branch N3 ) consisting of stable nodes of(3.2). Since the right-hand side of (3.4) is negative on M1 , the orbit moves
downward along M1 and reaches the transcritical point. Then, the orbit
leaves the transcritical point and returns to Σ. In this manner, we can deﬁnea Poincaré map Π : U → Σ.
In order to investigate the properties of the Poincaré map Π, take arectangle R on U ⊂ Σ and consider how it behaves when it runs along
the ﬂow of (3.1). The rectangle R is folded into a ring-shaped domain R
when moving upward around the green branch on M2 , and the radius of R isdecreasing. If ε is suﬃciently small, then R stays near the green branch of M2
for a long duration because the dynamics of w deﬁned by (3.4) is very slow.
In this case, R suﬃciently shrinks when passing through the vicinity of M2 ,and hence, Π(R) ⊂ R holds (Fig.11(a)). This implies that Π is a contraction
map, and thus, (3.1) has a stable periodic orbit (Fig.9 (a)). In contrast, when
ε is not so small, R cannot shrink suﬃciently. Therefore, Π(R) transversallyintersects R and a horseshoe is formed, as shown in Fig.11(b). Therefore,
the dynamics of (3.1) becomes chaotic through the cascade bifurcation as εincreases (see [7] for the proof).
In this section, we introduced the small parameter ε to apply the theory
described in [7] and revealed the mechanism of the onset of mixed-mode
oscillations and chaos of the fast-slow system. Nevertheless, we emphasizethat the mechanism of the onset of mixed-mode oscillations and chaos forthe original system (1.1) (i.e. ε = 1.0) is similar to that of the fast-slow
17system. In fact, Figs.3 and 4 show that dw/dt is much smaller than dp/dt
and dz/dt along the orbits. This implies that (1.1) can be regarded as the
fast-slow system without introducing the small parameter ε. Furthermore,
dw/dt along the chaotic attractor of (1.1) (see Fig.4) is larger than that along
the mixed-mode oscillation of (1.1) (see Fig.3) on average. Thus, the mixedmode oscillation of (1.1) becomes chaotic by the same mechanism discussed
above.

4

Concluding remarks

In this study, we show that the dormancy of predators can induce mixedmode oscillations and chaos in the population dynamics of a prey-predator
system under certain conditions. Noting that the time derivative of the density of dormant predators is smaller than those of prey and active predators,
we introduced the fast-slow system to show that a mixed-mode oscillation
and a chaotic attractor bifurcate from a coexisting equilibrium. This approach enables us to understand the mechanism of the onset of mixed-mode
oscillations and chaos in a prey-predator system with dormancy of predators.
To the best of our knowledge, mixed-mode oscillations and chaos have notyet been experimentally demonstrated in a simple phytoplankton-zooplankton(resting eggs) community in a microcosm with a short duration. We expectthat such complex behavior will be experimentally demonstrated in the future
based on our theory. In fact, as reported in [26], the period of a prey-predator
cycle in Daphnia-algal systems was 21.4 days on average under certain experimental conditions. On the other hand, the period of a prey-predator cyclein our model (1.1) is approximately 20.9 (day) under the parameter values
listed in Table 1 and K ≈ 5.17 (this value corresponds to the limiting pointon the bifurcating branch of periodic orbits in the bifurcation diagram of (1.1)
shown in Fig. 1(a)). This suggests that the experimental prey-predator cyclein Daphnia-algal systems corresponds to the theoretical one in our model.
Thus, it is expected that the Daphnia-algal prey-predator cycle will become
unstable and the systems will exhibit mixed-mode oscillations and chaos if
the experimental conditions are varied because in our model, mixed-mode
oscillations and chaos bifurcate from the prey-predator cycle. Our model
may serve as a reference for empirical researchers to ﬁnd a good example ofchaos in plankton communities, and it may contribute toward understanding the ecological implications of dormancy for short-lived organisms suchas zooplankton, though it should be emphasized that our theoretical prediction concerns the qualitative properties of population dynamics in a simple
phytoplankton-zooplankton(-resting eggs) community.
Complex dynamics such as mixed-mode oscillations and chaos has also
been reported in many literatures concerning food webs consisting of three ormore species (prey-predator is a basic building block of a food web). Roughly18

speaking, they are concerned with chaotic dynamics in (extensions of) a
three-species model of the food-prey-predator type known as the HastingsPowell model [18]. The reader should consult classical and recent papers[9, 12, 18, 25, 30, 31] and the references therein. [9, 12, 24] introduced fastslow systems to investigate bifurcation structures. However, the mechanismof generating mixed-mode oscillations and chaos in our system is diﬀerent
from those in [9, 12, 24]. For example in [9, 12], the dimension of the fast
direction is one whereas it is two in our system. Our analysis of the fastslow system is based on the theory described in [7], in which the existence of
mixed-mode oscillations and chaos is proved for a three-dimensional fast-slow
system with a two-dimensional fast direction.
Our model requires some extensions in order to reﬂect more realistic conditions. For example, to consider intraspecies interaction (density-dependent
inhibition) among active predators, we propose the system

p
dp


= r(1 − )p − f (p)z


dt
K



dz
(4.1)
= k1 µ(p)f (p)z + αw − d1 z − γz 2

dt





 dw = k2 (1 − µ(p))f (p)z − αw − d2 w,
dt
where γ is a positive constant. The eﬀects of intraspecies interaction among
predators for the MacArthur-Rosenzweig model and the Hastings-Powellmodel have been studied by [4, 28] and [10, 30, 37], respectively (an extension of the MacArthur-Rosenzweig model with intraspecies interaction
among predators is called the Bazykin model). From [10, 30, 37], it appears
that (4.1) have chaotic attractors for a narrow range of parameter values due
to a stabilizing eﬀect of intraspecies interaction. Combining singular perturbation methods and scaling arguments, the diﬀerence between (1.1) and(4.1) should be clariﬁed by detailed bifurcation analysis with the aid of numerical computations. It is one of the important problems to be consideredin further studies.
It may also be proﬁtable to seek other possibilities to generalize the model(1.1) and explore their bifurcation structures. This provides many pieces ofuseful information for understanding how the dormancy of predators aﬀects
the population dynamics of prey-predator systems.
Dormancy has begun attracting considerable interest among freshwater
ecologists, and recent empirical studies have provided fundamental and profound knowledge about the dormancy of zooplankton [2, 6, 15, 16, 27, 32].
Moreover, dormancy appears to have evolved independently among a wide
variety of living organisms [1, 3, 8, 21, 35]. It is an important adaptive
response to selective pressures that are common to a large number of organisms. Therefore, theoretical approaches such as those presented in this paper
are indispensable for understanding the ecological implications of dormancy.
19

Acknowledgments The authors express their sincere gratitude to Dr.
Takefumi Nakazawa, Professors Toshiyuki Ogawa and Kunimochi Sakamoto
for their useful advice and comments. Moreover, the authors would like toappreciate the referees for their useful suggestions and comments, which have
improved the original manuscript. The ﬁrst author was supported in partby the Grant-in-Aid for Scientiﬁc Research (C) No.18540120 JSPS, and thesecond author by JSPS Fellowships for Young Scientists.

References[1] Alekseev, V.R., De Stasio, B.T., Gilbert, J.J., Diapause in Aquatic Invertebrates: Theory and Human Use, Kluwer Academic (2007)[2] Alekseev, V., Lampert, W., Maternal control of resting-egg productionin Daphnia, Nature 414 (2001), 899-901.
[3] Baskin, C.C., Baskin, J.M., Seeds - Ecology, Biogeography, and Evolution of Dormancy and Germination, Academic Press (1998)[4] Bazykin, A.D., Nonlinear Dynamics of Interacting Populations, World
Scientiﬁc (1998)[5] Beninca, E., Huisman, J., Heerkloss, R., Johnk, K.D., Branco, P., Van
Nes, E.H., Scheﬀer, M., Ellner, S.P., Chaos in a long-term experimentwith a plankton community, Nature 451 (2008) 822-825.
[6] Carvalho, G.R., Hughes, R.N., Eﬀect of food availability, female culturedensity and photoperiod on ephippia production in Daphnia magna
Strauss (Crustacea: Cladocera). Freshw. Biol. 13 (1983), 37-46.
[7] Chiba, H., Fast-slow systems with Bogdanov-Takens type fold points,
submitted.
[8] Danks, H.V., Insect dormancy : an ecological perspective, Entomological Society of Canada (1987)[9] Deng, B., Food chain chaos with canard exposition, Chaos 14 (2004)
1083-1092[10] Deng, B., Equilibriumizing all food chain chaos through reproductive
eﬃciency, Chaos 16 (2006) 043125[11] Doedel, E.J., Champneys, A.R., Fairgrieve, T.F., Kuznetsov, Y.A.,
Oldeman, B.E., Paﬀenroth, R.C., Sandstede, B., Wang, X., AUTO 2000:
Continuation and bifurcation software for ordinary diﬀerential equations(with HomCont) (2000).
20

[12] Duarte, J., Januário, C., Martins, N., Topological invariants in the studyof a chaotic food chain system, Chaos 18 (2008) 023109.
[13] Genkai-Kato, M., Yamamura, N., Unpalatable prey resolves the paradoxof enrichment. Proc. R. Soc. Lond. B. 266 (1999), 1215-1219.
[14] Genkai-Kato, M., Yamamura, N., Proﬁtability of prey determines the
response of population abundances to enrichment. Proc. R. Soc. Lond.
B. 267 (2000), 2397-2401.
[15] Gyllström, M., Hansson, L.-A., Dormancy in freshwater zooplankton:
Induction, termination and the importance of benthic-pelagic coupling.
Aquat. Sci. 66 (2004), 274-295.
[16] Hairston, N.G.Jr., Hansen, A.M., Schaﬀner, W.R., The eﬀect of diapause emergence on the seasonal dynamics of a zooplankton assemblage,
Freshw. Biol. 45 (2000), 133-145.
[17] Hairston, N.G.Jr., Van Brunt, R.A., Kearns, C.M., Age and survivorshipof diapausing eggs in a sediment egg bank, Ecology 76 (1995), 1706-1711.
[18] Hastings, A., Powell, T., Chaos in a three-species food chain, Ecology
72 (1991), 896-903.
[19] Hoppensteadt, F.C., Izhikevich, E.M., Weakly connected neural networks, Springer-Verlag, New York (1997)[20] Jones, C.K.R.T, Geometric singular perturbation theory, Dynamicalsystems (Montecatini Terme, 1994), pp.44-118, Lecture Note in Math.
1609, Springer-Verlag, Berlin (1995)[21] Kaprelyants, A., Gottschal, J.C., Kell, D., Dormancy in non-sporulating
bacteria, FEMS Microbiol. Rev. 10 (1993), 271-285.
[22] Kuwamura, M., Nakazawa, T., Ogawa, T., A minimum model of preypredator system with dormancy of predators and the paradox of enrichment, Journal of Mathematical Biology 58 (2009), 459-479.
[23] Kuznetsov, Y.A., Elements of Applied Bifurcation Theory (3rd edition),
Springer-Verlag (2004).
[24] Kuznetsov, Y.A., Rinaldi, S., Remarks on food chain dynamics, Math.
Biosci. 133 (1996) 1-33.
[25] McCann, K., Yodzis, P., Bifurcation structure of a three-species foodchain model, Theor. Popl. Biol. 48 (1995), pp.93-125.

21

[26] McCauley, E., Nelson, W.A., Nisbet, R.M., Small-amplitude cycles
emerge from stage-structured interactions in Daphnia-algal systems. Nature 455 (2008), pp.1240-1243.
[27] McCauley, E., Nisbet, R.M., Murdoch, W.W., de Roos, A.M., Gurney,
W.S.C., Large-amplitude cycles of Daphnia and its algal prey in enrichedenvironments. Nature 402 (1999), pp.653-656.
[28] McGehee, E.A., Schutt, N., Vasquez, D.A., Peacock-López, E., Bifurcations, and temporal and spatial patterns of a modiﬁed Lotka-Voltera
model, International Journal of Bifurcation and Chaos 18 (2008), 22232248.
[29] Mischenko, E., Rozov, N., Diﬀerential equations with small parameters
and relaxation oscillations, Plenum Press New York (1980)[30] Peet, A.B., Deutsch, P.A., Peacock-López, E., Complex dynamics ina three-level trophic system with intraspecies interaction, Journal of
Theoretical Biology 232 (2005), 491-503.
[31] Rai, V., Upadhyay, R.K., Evolving to the edge of chaos: chance ornecessity? Chaos, Soliton and Fractals 30 (2006), 1074-1087.
[32] Ricci, C., Dormancy patterns in rorifers. Hydrobiologia 446 (2001), 111.
[33] Rosenzweig, M.L., MacArthur, R.H., Graphical representation and stability conditions of predator-prey interactions. Am. Nat. 47 (1963), 209223.
[34] Scheﬀer, M., de Boer, R.J., Implications of spatial heterogenety for theparadox of enrichment, Ecology 76 (1995), 2270-2277.
[35] Stewart, T.H.M., Wheelock, E.F., Cellular Immune Mechanisms andTumor Dormancy, CRC Press, (1992)[36] Vos, M., Kooi, B.W., DeAngelis, D.L., Mooij, W.M., Inducible defences
and the paradox of enrichment, Oikos 105 (2004), 471-480.
[37] Xu, C., Li, Z., Inﬂuence of intraspeciﬁc density dependence on a threespecies food chain with and without external stochastic disturbances,
Ecol. Model. 155 (2002), 71-83.

22