The robustness of cyclic dominance under random fluctuations

Journal of Theoretical Biology 308 (2012) 79–87

Contents lists available at SciVerse ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

The robustness of cyclic dominance under random fluctuations Christoph K. Schmitt a,n, Christian Guill b, Barbara Drossel a a b

Institut f¨ ur Festk¨ orperphysik, TU Darmstadt, Hochschulstraße 6, 64289 Darmstadt, Germany Institute for Zoology and Anthropology, University of G¨ ottingen, Berliner Straße 28, 37073 G¨ ottingen, Germany

H I G H L I G H T S c c c c

We apply noise to different stages of the sockeye life cycle. Cyclic dominance is very robust under noise, with phase jumps occurring rarely. Noise can synchronize the oscillation of different sockeye populations. Noise can induce slow predator oscillations.

a r t i c l e i n f o

abstract

Article history: Received 16 March 2012 Received in revised form 25 May 2012 Accepted 25 May 2012 Available online 4 June 2012

We investigate the influence of random perturbations on a recently introduced three-species model that reproduces the empirically observed pattern of cyclic dominance in Fraser River sockeye salmon. Since the sockeye populations are subject to various types of fluctuations affecting their growth and survival, we investigate the robustness of the model under several types of noise. In particular, we evaluate the variation of population sizes around their values in the deterministic model, the frequency of phase shifts in the 4-year oscillation, the extent of synchronization between different sockeye populations, and the response to strong one-time perturbations. Our main conclusion is that cyclic dominance is very stable even under strong noise in this model. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Sockeye salmon Predator–prey cycles Oncorhynchus nerka Noise

1. Introduction The 4-year oscillations of the number of sockeye salmon (Oncorhynchus nerka) that return from the ocean to spawn in their native streams in the Fraser River basin in British Columbia (Canada) are a well-documented example of large-scale population oscillations (Ricker, 1950, 1997; Townsend, 1989). They belong to the class of single-generation cycles (Gurney and Nisbet, 1985; Murdoch et al., 2002). Every fourth year (which is the dominant generation time) the abundance of these fish is at very high levels, reaching several million fish in some spawning populations, but drops to numbers between several hundred and a few ten thousand individuals in the following years (Fig. 1). This characteristic pattern is often referred to as cyclic dominance. The oscillations were reported as early as the 19th century and are evident for instance in the extremely high catches by fisheries every fourth year (Rounsefell and Kelez, 1938). The phenomenon

n

Corresponding author. Tel.: þ49 6151 16 3106; fax: þ 49 6151 16 3681. E-mail addresses: [email protected] (C.K. Schmitt), [email protected] (C. Guill), [email protected] (B. Drossel). 0022-5193/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jtbi.2012.05.028

has important implications for the management and conservation of this both economically and ecologically important salmon species. Recently, Guill et al. (2011a) suggested a three-species model showing that cyclic dominance can be obtained as a population dynamics attractor that is due to predator–prey interactions occurring in the rearing lakes. The model produces cyclic dominance over a wide range of parameter values, but other dynamical patterns, such as fixed points, period-2 oscillations, quasiperiodic oscillations, or chaos are also found (Guill et al., 2011b). The model includes population dynamics of the sockeye fry, their predator (mainly rainbow trout), and their preferred food (zooplankton, mainly daphnia) in the lakes where the sockeye fry live for one season after hatching. Then, they migrate to the ocean, where they stay until they return to spawn. After spawning the salmon die, the important input of marine-derived nutrients to the freshwater ecosystems by their decomposing carcasses (Larkin and Slaney, 1997; Johnston et al., 2004; Hocking and Reynolds, 2011) also being accounted for. The view that the relevant processes generating the oscillations take place in the freshwater is supported by the fact that nowadays the oscillations of sockeye salmon originating from different lakes are out of phase. If the

C.K. Schmitt et al. / Journal of Theoretical Biology 308 (2012) 79–87

total escapement / millions

6 5 4 3 2 1 0 1950

1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 1950

1960

1970

1980 Year

1990

2000

total escapement / millions

total escapement / millions

total escapement / millions

80

1960

1970

1980 Year

1990

2000

4 3.5 3 2.5 2 1.5 1 0.5 0 1950

1960

1970

1980 Year

1990

2000

1960

1970

1980 Year

1990

2000

1.2 1 0.8 0.6 0.4 0.2 0 1950

Fig. 1. Four examples for cyclic dominance. (a) The salmon escapement (number of female and male salmon that escape the fisheries) of the late Shuswap Lake run. (b) The same data for Quesnel Lake, (c) for the late Stuart run and (d) for Chilko lake (Pestal et al., 2008).

cause of the oscillations was in the ocean, this cause would affect different spawning populations in the same way, leading to a synchronization of the oscillations. While the predator–prey model (Guill et al., 2011a) captures the core mechanism that generates cyclic dominance, due to its deterministic nature it misses important effects that shape the time series of natural spawning populations. The salmon populations are subject to ubiquitous fluctuations in their natural environment, which may either increase or decrease the magnitude of the oscillations, or may even cause phase jumps (a shift of the dominant line from one brood line to another) (Myers et al., 1998). In particular human activities such as fisheries or temporal blocking of migration routes cause severe perturbations. For several years, the catch quotas for adult salmon were as high as 80% of the total migrating population (Ricker, 1950). Despite these disturbances and the large random fluctuations that are superimposed to the dynamics, several spawning populations show remarkably stable oscillations over more than 70 years. The sockeye escapement data shown in Fig. 1 illustrate the stability of the dominant line: even though it may temporarily be lower than another line (as for instance in 1960 and 1964 in the Chilko lake data), it often becomes again the strongest line afterward. Escapement data of several lakes show instances where the subdominant line becomes temporarily stronger than the dominant line (in the examples shown here this can be seen in 2002 in the Quesnel lake data). Other striking features of these data are the huge difference between the dominant and weak years (with a factor of more than 100 not being unusual) and variations by more than a factor of two between subsequent returns of the same brood line (i.e., after 4 years). (The complete breakdown of the cyclic dominance pattern in Chilko lake after 1990 is due to lake fertilization.) In order to gain a better understanding of stochastic effects in this unique ecological system, we therefore expand the predator– prey model for sockeye salmon in this paper by analyzing the stability of cyclic dominance under stochastic perturbations. We apply noise at different stages of the life cycle of the sockeye salmon and find that phase jumps occur rarely in the model, just

as observed in the empirical data. Even though they are rare, phase jumps become more frequent as the level of noise increases. Additional insight is gained by applying a large onetime perturbation to the otherwise deterministic model. Both types of investigations show that the time needed for a phase jump to become established is long compared to the generation time of the sockeye salmon and extends from several decades up to 100 years. Since different spawning populations may encounter similar stochastic effects (for instance due to regional weather conditions affecting several lakes in the same way, or to increased river temperature during migration), we investigate also the capability of noise to synchronize these populations. By starting with several model populations that have their population maximum in different years and by subjecting all these populations to the same stochastic perturbations, we find that over time their oscillations can synchronize. This may represent the state of the natural system before the sockeye populations in the Fraser River were nearly driven extinct by human activities in the early 20th century (Ricker, 1950).

2. Model 2.1. The original deterministic model The most important mechanism for creating cyclic dominance is a negative feedback of a population at one moment in time on this population at a later moment in time. This is a generic property of predator–prey models: larger than average prey numbers lead to predator growth, which in turn increases prey consumption and reduces the prey population at a later moment in time. The predator–prey model introduced in Guill et al. (2011a) was designed to include only the essential mechanisms required for the occurrence of cyclic dominance, not to quantitatively predict the population dynamics of all species in the corresponding lake ecosystems. The model uses standard continuous population dynamics equations for sockeye fry, sn(t), their predator, pn(t), and their zooplankton food, zn(t), during the

C.K. Schmitt et al. / Journal of Theoretical Biology 308 (2012) 79–87

d zn ðtÞ  sn ðtÞ sn ðtÞ  pn ðtÞ sn ðtÞ ¼ l  asz aps ds  sn ðtÞ, dt 1 þ cs  sn ðtÞ þ zn ðtÞ 1 þ cp  pn ðtÞ þsn ðtÞ   d zn ðtÞ zn ðtÞ  sn ðtÞ zn ðtÞ ¼ r  zn ðtÞ  1 , asz dt Kn 1þ cs  sn ðtÞ þzn ðtÞ d sn ðtÞ  pn ðtÞ p ðtÞ ¼ l  aps dp  pn ðtÞ: dt n 1 þ cp  pn ðtÞ þsn ðtÞ

ð1Þ

sn, zn and pn are biomass densities. They are scaled such that the half-saturation density of the sockeye-plankton predation term is 1. The parameter choice is based on the bioenergetics approach by Yodzis and Innes (1992), in the form used by Brose et al. (2006) and by Heckmann et al. (2012). Per capita rates of respiration, consumption, and mortality scale with body mass to the power 0.75 (Brown et al., 2004). By estimating that the predator body mass is approximately 1500 times larger than the fry body mass, we obtain the following parameter values (Guill et al., unpublished), which are used throughout this paper: l ¼ 0:85 is the food to predator biomass conversion rate of carnivores. The maximal per unit biomass ingestion rates are asz ¼ 12 for salmon fry, and aps ¼ 2 for rainbow trout. The parameters ds ¼3 and dp ¼0.5 denote biomass loss from respiration and mortality. Using the same allometric scaling laws it can be shown that the time scale of species dynamics scales with body mass to the power of 0.25 (Kartascheff et al., 2010), so the dynamics of the rainbow trout are about six times slower than those of the sockeye fry. The feeding terms are based on a Beddington functional response (Beddington, 1975). They include saturation at high prey densities. In contrast to the more widely used Holling type II functional response, the Beddington functional response also includes a predator interference term, motivated by the observation that fish affect each others’ foraging behavior (Walters and Christensen, 2007). By applying again allometric scaling with body mass and choosing the interference parameters such that they have a noticeable effect on the dynamics but do not become dominant, we obtain cs ¼0.8 and cp ¼0.25. Logistic growth is used for the zooplankton, with maximum growth rate r ¼15 and carrying capacity Kn, which depends on the number of spawning sockeye and is given by   sn ð0Þ : ð2Þ Kn ¼ K0 þ k k0 þ sn ð0Þ The three parameters determining the carrying capacity are the base carrying capacity K 0 ¼ 20, the maximal added capacity due to the nutrient input by sockeye carcasses k ¼ 20, and the half-saturation density k0 ¼ 10. Measurements show that in some rearing lakes in years with many spawning salmon the plankton carrying capacity is nearly doubled (Hume et al., 2004), so K0 and k should be of the same order of magnitude, and k0 should be somewhat smaller than the maximum value of sn seen in the model.

35 30 values of sn (T)

growth season from spring (t ¼0) to fall (t ¼ T ¼ 1) in year n. This is combined with a rule for calculating the three populations at the start of each season based on the populations at the end of previous seasons. The sockeye population at the end of a season determines the number of spawning adults (and thus the initial fry biomass) 3 years later. A small fraction of sockeye stay in the ocean for one more year, causing a mixing between brood lines. The spawning sockeye represent a nutrient input that affects the lake’s zooplankton carrying capacity. This means that more food is available to the sockeye fry in strong sockeye years. The differential equations governing the continuous part of the model are

81

25 20 15 10 5 0

0

5

10

15 K0

20

25

30

Fig. 2. Bifurcation diagram of the three-species model. In black, the values of sn(T) on the attractor are plotted versus the base plankton carrying capacity K0. At low K0 the attractor is a fixed point, so there is only one value of sn(T). Around K0 ¼ 7.5 there is a Neimark–Sacker bifurcation, and salmon biomass starts to oscillate quasiperiodically. In this regime infinitely many different values of sn(T) are possible. Around K 0 ¼ 12 phase locking occurs, and the oscillation period becomes fixed at four years. This means that only four different values of sn(T) occur on the attractor. In gray 100 biomass values occurring with noise of type A (s ¼ 0:3, see below) are shown for each value of K0.

From empirical data, we concluded that effects of previous years on the zooplankton population size can be neglected (Levy and Wood, 1992), and we therefore initialized the zooplankton biomass in the spring of year n with the value Kn. The plankton carrying capacity varies between lakes, and depending on the value of K0 the model produces different dynamical attractors (Fig. 2). For low carrying capacity (which would reflect the situation in ultraoligotrophic lakes), the dynamics goes to a fixed point, while it shows a very pronounced pattern of cyclic dominance at larger K0. For the value of K0 used in our simulations, the system is well within the regime of cyclic dominance (Guill et al., unpublished). Varying other parameters leads to similar bifurcation diagrams (Guill et al., 2011a, unpublished). The matching conditions that determine the spring biomass densities from the biomass densities at the end of previous seasons are given by sn þ 1 ð0Þ ¼ gðð1EÞsn3 ðTÞ þ Esn4 ðTÞÞ, zn þ 1 ð0Þ ¼ K n þ 1 , pn þ 1 ð0Þ ¼ pn ðTÞ,

ð3Þ

with E ¼ 0:1 denoting the fraction of adult sockeye spawning at age 5. Predator biomass does not change between fall and spring in this model (but this property is not essential for the results). The parameter g ¼ 0:1 summarizes all factors influencing smolt to fry biomass conversion, including smolt survival, ocean survival, spawning success and egg to fry survival. It is important to note that no equilibrium is reached in the continuous part of the dynamics: at the start of each season the sockeye biomass is set to a small value, and the zooplankton biomass to a high value. The season ends before an equilibrium is reached. At a fixed point of the dynamics, 1=g is equal to the factor of salmon biomass growth during each season, and its value could therefore be estimated from empirical data. Obviously, setting this parameter to a constant value is a gross simplification, since there may occur large fluctuations in all factors that influence g. This model reproduces the periodic oscillation of spawning sockeye, but due to its deterministic nature it does not show the fluctuations visible in the empirical data (compare Figs. 1 and 3a). The oscillation with period 4 was shown to be caused by a socalled ‘‘strong resonance’’ near a Neimark–Sacker bifurcation (Guill et al., 2011a) (Fig. 2). It is an attractor of the model over a wide range of parameter values (Guill et al., unpublished).

82

C.K. Schmitt et al. / Journal of Theoretical Biology 308 (2012) 79–87

20

25 20 snA(T)

sn(T)

15 10 5

10 5

0

0 0

10

20 30 Time (years)

40

50

25

20

20

15

15

A sn (T)

25

A

sn (T)

15

10 5

0

10

20 30 Time (years)

40

50

0

10

20 30 Time (years)

40

50

10 5

0

0 0

10

20 30 Time (years)

40

50

Fig. 3. Simulated salmon abundance including different amounts of noise. Noise is modeled by multiplying the spawning sockeye biomass with a lognormally distributed random variable with a width s. Shown here are data obtained with noise of type A (see below), data for other noise types are similar. The values of s used are (a) 0 (no noise), (b) 0.2, (c) 0.5 and (d) 1.5.

2.2. Model with noise In order to analyze the effects of fluctuating conditions in the ocean and during migration on cyclic dominance, noise was introduced into the discrete part of the model. Four types of noise were examined. In our computer simulations, we used only one kind of noise at a time in order to see more clearly the different effects caused by different types of noise. We modeled noise by multiplying g or E in year n with a lognormally distributed random variable xn. The probability density for xn is thus given by pðxn Þ ¼

1 ðln xn Þ2 pffiffiffiffiffiffi exp , xn 4 0: 2s2 xn s 2p

ð4Þ

The logarithm of xn is normally distributed with mean 0 and standard deviation s. Using multiplicative noise with a lognormal distribution means that doubling and halving the value of g or E occurs with the same probability. By varying s, the strength of the noise can be modified. Strong currents or warm water during migration to the spawning grounds reduce spawning success of adult sockeye (Farrell et al., 2008), which is contained in the parameter g. This random perturbation of spawning success uniformly affects all adults spawning in a given year, both of age 4 and 5, and we call this type of noise ‘‘type A’’ in this paper. Factors that affect egg to fry survival lead also to type A noise. To model the effect of this type of noise, Eq. (3) is replaced by sAn þ 1 ð0Þ ¼ xn gðð1EÞsAn3 ðTÞ þ EsAn4 ðTÞÞ:

ð5Þ

We did not perturb the continuous part of the model or the initial predator and plankton biomass. After migrating from the rearing lake to the ocean, sockeye smolts live near the river mouth for some time to adjust to salt water. During this time, they are particularly susceptible to food scarceness or diseases. Such changing conditions during and after smolt migration are modeled by noise that uniformly affects all smolts leaving the lake (type B), no matter whether they will return to spawn at age four or five. As for type A (and type C and D), only a change to equation (3) is required to incorporate this into the model sBn þ 1 ð0Þ ¼ gðð1EÞxn3 sBn3 ðTÞ þ Exn4 sBn4 ðTÞÞ:

ð6Þ

Sockeye salmon live in the ocean for about 3 years, during each of which conditions can be different. All salmon which are in the ocean at the same time are affected by fluctuating conditions in the ocean. Noise in the ocean (type C) thus affects all brood lines that are in the ocean in the same way. Eq. (3) is changed to sCn þ 1 ð0Þ ¼ gðð1EÞxn1 xn2 xn3 sCn3 ðTÞ þ Exn1 xn2 xn3 xn4 sCn4 ðTÞÞ: ð7Þ Finally the generation length of sockeye salmon may vary. More or fewer salmon can spawn at age 5, which is modeled by varying the parameter E. This generation length noise (type D) causes smaller population fluctuations than the other three types. Eq. (3) is changed to assign a different generation length to each sockeye brood: D D sD n þ 1 ð0Þ ¼ gðð1xn3 EÞsn3 ðTÞ þ Exn4 sn4 ðTÞÞ:

ð8Þ

A similar distinction between different types of noise is also made by other authors modeling salmon dynamics. In Worden et al. (2010), noise similar to our types B, C, and D was investigated. For larger s, the brood lines show larger fluctuations, as can be seen in Fig. 3b–d. These data also show that the dominant brood line may temporarily become weaker than another line, and that the relative strength of the other lines can also show large variations.

3. Results 3.1. Effects of continuous noise Introducing noise into the system has several effects: (i) a slow periodic oscillation of the predator biomass is induced. (ii) The biomasses of the four salmon lines fluctuate around their values on the deterministic attractor. (iii) Phase jumps (i.e., shifts of the dominant line to another year) are induced. The magnitude of each of these effects depends on noise strength s and the type of noise. In the following, we discuss these three effects one after the other. (i) Slow predator oscillations: All types of noise cause the predator biomass to oscillate with a period of about 25–35 years. The amplitude of this slow oscillation is smaller than that of the

20

20

15

15 pn (T)

pn (T)

C.K. Schmitt et al. / Journal of Theoretical Biology 308 (2012) 79–87

10 5

83

10 5

0

0 0

10

20

30 40 Time (years)

50

60

70

0

10

20

30 40 Time (years)

50

60

70

Fig. 4. Oscillating predator abundance. Slow predator oscillations in the model without (left) and with noise (right). In the absence of noise, the predator oscillation is a transient phenomenon. In the presence of noise (here type A with s ¼ 0:5) the oscillation is permanently excited.

30 25 values of snA(T)

4-year oscillation. This oscillation occurs also in the deterministic model, but is only a transient phenomenon because it is dampened. Both cases are depicted in Fig. 4. Mathematically, this oscillation is due to a pair of complex-conjugate eigenvalues characterizing the dynamics close to the bifurcation (Guill et al., 2011a). Salmon biomass is only affected weakly by this predator oscillation. Excitation of dampened oscillations by noise is a generic phenomenon in dynamical systems (Gang et al., 1993?). The reason is that the fourier spectrum of the noise contains also the frequency of the dampened oscillation, leading to a resonance. (ii) Biomass fluctuations of the four salmon lines: We evaluated the fluctuations of the salmon biomass around its average value for each brood line. Fig. 5 shows the increasing spread of biomass values with increasing spawning success noise (type A). The subdominant brood line fluctuates more strongly than the dominant one at equal s (also for other noise types). The amount of fluctuations can be quantified by measuring the coefficient of variation (the standard deviation divided by the mean) of each brood line. Because the average of a brood line can change due to a phase jump, we used moving averages over 80 years (20 per brood line, see Fig. 6) to determine the coefficient of variation. Finally, the coefficients of variation of the dominant and subdominant brood lines were averaged to quantify the overall variation of the two strongest brood lines (the two weak brood lines are ignored because of their negligible quantitative contribution to the total salmon population). The fluctuations grow larger with increasing noise, eventually suppressing the periodic behavior. Fig. 7 shows the increase of the fluctuation strength with the noise strength s, for all four types of noise. At equal s, the different types of noise cause different amounts of fluctuations. Doubling the fraction of sockeyes that spawn at age 5 has a much smaller effect on the dominant and subdominant brood line than doubling the number of spawning sockeyes or smolts. Noise types A and B lead to very similar effects. The effect of oceanic noise (type C) is largest, because each brood line experiences this effect several times. To determine at what level of noise cyclic dominance breaks down we use the following criterion: if the smallest of the moving averages over brood lines is larger than one-third of the largest moving average in at least 20% of 4-year periods, cyclic dominance is no longer present. This occurs around s ¼ 1:3 for spawning success noise (A), s ¼ 1:4 for smolt survival noise (B), s ¼ 1:5 for noisy generation length (D), and already at s ¼ 0:8 for oceanic noise (C). Plots are cut off at these points because our methods of analyzing data with moving averages over the separate brood lines can not be properly applied to time series without cyclic dominance. (iii) Phase jumps: Sufficiently strong noise can change the phase of the salmon oscillation, causing a change of the dominant broodline. Two time series of sockeye biomass during a phase

20 15 10 5 0 0

0.1

0.2

0.3 σ

0.4

0.5

0.6

Fig. 5. Spreading of biomass values with increasing noise. In black, the values of sn(T) on the deterministic attractor are plotted (for K 0 ¼ 20). In grey biomass values occurring with spawning success noise (type A) of varying strength are shown.

jump are shown in Fig. 6. The phase jump mechanism is the same for all noise strengths and types: the biomass of the brood line that will become dominant slowly rises to its new level, while the originally dominant one declines to the level of a weak or subdominant line. Any brood line can become the dominant one. It usually takes between 40 and 100 years for the dominant brood line to drop to the biomass level of a weak line, the drop to a subdominant line is significantly faster. Fig. 6b shows the frequency of phase jumps as a function of the noise strength for all four noise types. For small amounts of noise, there are no phase jumps. The dependence of the number of phase jumps on the type of noise is similar to that of the fluctuations: given the same noise strength s, oceanic noise (type C) causes the largest number of phase jumps, while perturbing the generation time (type D) leads to fewer phase jumps than other types of noise. Noise in spawning success (type A) and in smolt to adult survival (type B) produces similar numbers of phase jumps, which lie between the two other types of noise. The four examined types of noise produce different numbers of phase jumps not only at equal s but also at equal amounts of fluctuations, albeit to a lesser extent (Fig. 6c). This shows that the occurrence of phase jumps depends on the correlation between perturbations. Now the perturbations in generation length cause the largest number of phase jumps and noise in the ocean the smallest number of phase jumps. This means that moving biomass from one year to the next is more likely to cause a phase change than simply increasing spawner biomass. Changes in oceanic conditions cause the fewest phase jumps because the populations in subsequent years are similarly affected by them. The number of phase jumps seems to rise linearly with fluctuation strength once it crosses a certain threshold. The threshold depends on the noise type, but the slope of the linear part does not.

84

C.K. Schmitt et al. / Journal of Theoretical Biology 308 (2012) 79–87

25

snA(T)

20 15 10 5 0

0

20

40

60

80

100

60

80

100

Time (years) 25

s n A(T)

20 15 10 5 0

0

20

40 Time (years)

Fig. 6. Phase jumps and moving averages of sockeye biomass. The four smooth lines represent moving averages over 80 years of the four brood lines (represented by four different symbols). Spawning success noise (type A) with s ¼ 0:6 was used in (a), and with s ¼ 0:5 in (b). In these data sets phase jumps occur where the brood line represented by the filled squares changes from dominant to weak, and the first weak line (symbol *) becomes dominant. In (a) this happens in a single phase jump, in (b) the same result is accomplished with two backward jumps in quick succession.

In principle, three different types of phase jumps can occur. These are forward phase changes (the subdominant line becomes dominant and the dominant line becomes weak, phase þ1), backward phase changes (the dominant line becomes subdominant and the second weak line becomes dominant, phase þ3, example in Fig. 6b), and phase jumps of two years (the first weak line becomes dominant and the dominant line becomes weak, phase þ2, Fig. 6a). We evaluated the frequency with which these three types of phase jumps occur. With strong noise, phase jumps of 1, 2, or 3 years are about equally likely, while weak noise causes mostly backward phase jumps (Fig. 8). The strong preference for backward jumps at low noise strength is found for all noise types, but disappears for parameter values considerably closer to the quasiperiodic regime (i.e., for smaller K0).

examples for how the system responds to a one-time perturbation of the dominant and the subdominant brood line. Changing survival of spawning adults or smolt survival (again by a factor of up to 10) for any single brood line produced no lasting phase change. This remarkable stability of the deterministic system with respect to strong one-time perturbations depends on the simulation parameters. With parameters closer to the quasiperiodic regime (i.e., for smaller base plankton carrying capacity K0, see Fig. 2), lasting phase jumps could also be caused by other perturbations, such as very strong perturbations of the dominant brood line.

3.3. Synchronization of multiple spawning populations 3.2. Phase jumps after a strong one-time perturbation In addition to random continuous noise, sometimes strong or even catastrophic one-time events occur. For example the 1914 blockade of Hell’s Gate (a narrowing of the Fraser river canyon) nearly wiped out the Fraser sockeye and damaged salmon runs for decades. A less severe event occurred in 2009/2010: salmon returns in 2009 were unexpectedly low, while returns in 2010 were unexpectedly high (without an increase in the number of five year old salmon). To better understand what can happen after a strong perturbation of the system, single perturbations of the deterministic system from its attractor were examined. Without continuous noise the system returns to the period four attractor after the perturbation, but the phase can change. With the chosen parameters, such a phase jump occurs only rarely. In fact, we found phase jumps only for noise type D, i.e., by delaying a brood line, causing most of the salmon that would usually spawn in a given year to spawn one year later instead. Changes up to multiplication and division of the parameter E by ten during the year of the perturbation were examined for every brood line, but with the chosen parameters the system almost always returns to the same attractor phase. Only delaying the spawning of the subdominant brood line caused a phase jump, and only if the fraction of salmon that spawn at the age of five was set to at least 78% (E Z0:78). The two examples in Fig. 9 show

All four types of noise used in our study take place in the ocean and during migration. Therefore, noise should affect salmon populations from different spawning regions similarly. We investigated whether this can lead to the synchronization of different populations. First, we simulated 100 salmon populations by using exactly the same model parameters for each population. Only the initial population sizes were chosen differently. We applied to all 100 populations exactly the same noise, but otherwise the populations were completely independent from each other. When noise was so weak that no phase jumps can occur, all four possible phases of the oscillation occurred even after a long time. Populations with the same phase quickly became completely synchronized. When noise was strong enough that phase jumps did occur, eventually a full synchronization of all populations was reached (Fig. 10). This is because more and more populations that have the same phase become synchronized and have exactly the same dynamics. Once n populations have become synchronized, they stay synchronized forever, and they perform their phase jumps together. After such a phase jump, they may become synchronized with populations that have already their new phase. In this way, n increases with time. Eventually all populations have the same phase and oscillate synchronously.

C.K. Schmitt et al. / Journal of Theoretical Biology 308 (2012) 79–87

85

80 phase +1 phase +2 phase +3

1.2 percentage of phase jumps

1 fluctuation strength

70

noise type: spawning success (A) smolt survival (B) oceanic conditions (C) generation length (D)

0.8 0.6 0.4

40 30 20

0 0.4

0.6

0.8

1

1.2

1.4

σ

0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

σ

noise type: spawning success (A) smolt survival (B) oceanic conditions (C) generation length (D)

12000 phase jumps in 106 years

50

10

0.2

10000 8000 6000 4000

0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

σ

noise type: spawning success (A) smolt survival (B) oceanic conditions (C) generation length (D)

12000 10000

Fig. 8. Directions of phase change for different noise strength s. The direction of phase jumps changes with noise strength s, backward jumps (phase þ 3) are much more likely at low noise strength. The distribution shown is for spawning success noise (type A), it is very similar for other noise types.

their rearing lakes have different size and different nutrient content. We therefore performed additional computer simulations, where either noise or parameters were made different between different populations. Fig. 10b shows a simulation where in addition to identical noise of type A with s ¼ 0:5, each population experiences an additional, independent noise contribution of type A with s ¼ 0:1. This means that noise was similar for all populations, but not identical. The figure shows that for most of the time more than 25% of the populations have the same phase, and during some periods more than 90% have the same phase, but full synchronization is not achieved. Fig. 10c shows a simulation where the base plankton carrying capacity K0 was different for different populations, ranging from K 0 ¼ 15 to K 0 ¼ 25. Noise was identical for all 100 populations. Almost always more than 50% of all populations are in the same phase. Before the Hell’s Gate blockage, the Fraser sockeye populations were to a large extent synchronized (Ricker, 1950). This shows that after very long times partial synchronization of salmon runs does happen in nature, just as in the simulations shown in Fig. 10b and c.

2000

phase jumps in 106 years

60

8000 6000

4. Discussion

4000

We have investigated the effect of random fluctuations (‘‘noise’’) on cyclic dominance of pacific sockeye salmon, using the three-species model recently introduced by Guill et al. (2011a). We implemented four different types of noise, corresponding to the following different instances at which population dynamics can be affected by changing environmental conditions: migration to the ocean, survival in the ocean, growth in the ocean, and migration from the ocean. The main finding of this paper is that the salmon oscillation in the three-species model is very stable under noise. Even strong noise has little effect on the phenomenon of cyclic dominance, only causing fluctuations around the mean values of the four lines and infrequent phase jumps. Phase jumps take several decades to become fully established. Only noise with s considerably larger than 1 can destroy the cyclic dominance for most types of noise (see Discussion at the end of Section 3.1(ii)). After a strong onetime perturbation, the system usually slowly returns to the attractor, retaining its phase. However, sometimes such a perturbation can induce a phase jump, which may become visible only several decades after the perturbation. Different spawning populations that are subject to the same noise can become synchronized with respect to the dominant year. We also found that noise

2000 0 0

0.2

0.4

0.6

0.8

1

fluctuation strength Fig. 7. Relation between noise strength, fluctuation in population size, and number of phase jumps. Top: fluctuation strength vs noise strength s. Middle: number of phase jumps in 106 seasons vs noise strength s. Stronger noise leads to stronger fluctuations and more phase jumps. Bottom: number of phase jumps in 106 seasons vs fluctuation strength. Plots are cut off when cyclic dominance is fully suppressed by noise.

A complete synchronization of all salmon populations within the same river system is not observed in nature, and of course our computer simulations were unrealistic in several respects: first, the noise experienced by different populations is not completely identical because they do not migrate exactly at the same time and do not encounter exactly the same conditions in their spawning streams. Second, the values of population dynamics parameters are not identical for different populations, because

C.K. Schmitt et al. / Journal of Theoretical Biology 308 (2012) 79–87

25

20

20

15 sn (T)

sn (T)

86

15 10

10 5

5 0

0 0

10

20

30 Time (years)

40

50

60

0

20

40

60 80 Time (years)

100

120

140

100 80 60 40 20 0 0

400

800

1200

1600

2000

Synchronized Pops.

Synchronized Pops.

Synchronized Pops.

Fig. 9. Response of the system to a one-time perturbation of generation length. In (a), E was set to 0.95 (from 0.1) in year 9, a spawning year of the dominant brood line. This forces a temporary phase change of the salmon oscillation that slowly relaxes into the previous state. In (b) the spawning of the subdominant brood line was delayed instead, E was set to 0.80 in year 10. This causes the phase of the oscillation to shift in the opposite direction much later.

100 80 60 40 20 0 0

400

800

Time (years)

1200

Time (years)

1600

2000

100 80 60 40 20 0 0

400

800

1200

1600

2000

Time (years)

Fig. 10. Synchronization of salmon oscillations by noise. We show the maximum number of populations (out of 100) with the same phase, with an initial number of 25 populations in each phase. In (a) all salmon populations have the same parameters and are affected by identical noise of type A with s ¼ 0:5. This leads to complete synchronization. In (b) an additional noise of type A with s ¼ 0:1 is applied to each population independently. (c) uses again identical noise, but the populations have different base plankton carrying capacities K0 in the range between 15 and 25 instead of the default value of 20. In (b) and (c) there is no permanent, complete synchronization, but a strong bias towards oscillation with the same phase.

induces a slow oscillation of the predator species. Such an oscillation occurs in the deterministic model only as a transient, dampened oscillation. These observations can be understood when viewing the model as a dynamical system: in the absence of noise, the system has four attractors, because cyclic dominance can occur with four different phases, depending on which year is the dominant one. Each attractor has a basin of attraction, which is given by the set of all combinations of the initial population sizes for which the dynamics approach the considered attractor. Apparently, these basins are relatively ‘‘deep’’ and large, so that it is difficult to push the system by noise from one basin to another. When noise is small, it causes the population sizes to fluctuate around the values that they have on the attractor, but cannot push the system over the attractor boundary. Only when noise is large enough, it can cause the system to move to a different basin. When noise is applied only once, a phase jump can only occur if the perturbation is so large that it places the system in a different basin of attraction. When noise is applied continuously, as was done throughout most of our investigation, the system can be pushed to a different basin of attraction by an appropriate combination of many smaller fluctuations. Even after the boundary to a different attractor basin has been crossed, the system may be pushed back by the following fluctuations, as long it has not yet moved sufficiently far away from the boundary. This picture is supported by the slowness of the phase change. It suggests dynamical emergence instead of a series of perturbations directly changing biomass values to the new phase. Also the distribution of phase jumps in time is not uniform. Phase jumps seem to be more likely a short time after another phase jump. This indicates that the system is more susceptible to further phase jumps after a phase jump than during regular dynamics, which supports the idea that the dynamic state of the system is far from the stable attractors during phase jumps. The finding that phase jumps are so difficult to achieve and take such a long time in the model is due to the four brood lines

being coupled via the predator. Even if one brood line becomes strongly perturbed, the predator experiences the usual growth conditions in the other three years, and in this way the strong perturbation becomes attenuated with time. In order to achieve a phase jump, all four brood lines and the predator population must change their phase, and this can only occur if the perturbation moves all five population densities to a point that lies within the basin of attraction of the oscillation with the other phase. Let us compare our findings with the empirical data shown in Fig. 1: these data suggest that phase jumps are extremely rare in nature, as also observed in Myers et al. (1998). In spite of the rareness of phase jumps, the population size of each brood line fluctuates considerably, and the dominant brood line can become from time to time weaker than another brood line. The empirical data also shows that the strongest and weakest line within a fouryear period may vary by a factor of more than 100. All these features are reproduced by our model at realistic levels of noise, due to the fact that the ‘‘strong resonance’’ responsible for cyclic dominance in our model is an extremely robust attractor. Unfortunately, there is very little data about the rainbow trout population, so the low-frequency predator oscillation could not be verified. Because our model uses uncorrelated noise, long-term changes in sockeye population dynamics, such as the decline of the Stuart lake population during the last 20 years, are not captured by the model. Earlier models for salmon oscillations (Myers et al., 1998; Worden et al., 2010), which do not include the coupling to a predator and in which the oscillation can be sustained only in the presence of noise, do not match the features of the empirical data as well as our model. If cyclic dominance is not due to the coupling of the salmon populations to another species, the changes caused by noise persist much more easily from one generation to the next, as is the case in simulations where cyclic dominance is obtained via stochastic forcing of the Ricker model (Myers et al., 1998). In such simulations phase jumps are much more frequent, and the

C.K. Schmitt et al. / Journal of Theoretical Biology 308 (2012) 79–87

dominant–subdominant–weak–weak pattern frequently observed in nature occurs only episodically, if at all. Also, the difference between the strongest and weakest line varies more strongly and is typically much smaller than in nature. The model by Worden et al. (2010) was not designed to study cyclic dominance of sockeye salmon, but to examine the effects of environmental noise on chinook and coho salmon at several stages of their life history. Depending on the type of salmon, it has a 3- or a 4-year generation time. The type of investigation is similar to the one performed in Guill et al. (2011a), with the essential difference that there is no coupling to other species. The simulation data shown in this paper display oscillations corresponding to the dominant generation time, which are excited by noise. Just as in Myers et al. (1998), these oscillations are much weaker than those displayed by our model. Interestingly, Worden et al. (2010) also found that the stage in the life history at which noise is applied significantly changes the system’s reaction. We conclude that for cyclic dominance to be as pronounced and robust as it is in nature, the population dynamics must be on an attractor with a large basin of attraction. The three-species model by Guill et al. (2011a) has this feature, while other models that generate salmon oscillations do not have this property and therefore lead to less regular, less pronounced, and less stable oscillations. In this respect, the model by Guill et al. (2011a) resembles another three-species model, which describes the cyclic outbreaks of the spruce budworm (Ludwig et al., 1978). This model contains also only a minimal set of ingredients that are required for producing the oscillation, and it has a very stable attractor, leading to strong resilience of the system (Holling, 1998). Certainly, more models of this type are needed in order to better understand what stabilizes ecosystems.

Acknowledgments This work was supported by the DFG under contract number Dr300/9-1. We thank Eddy Carmack for useful discussions, and Carl Walters for pointing out to us realistic parameter values. References Beddington, R., 1975. Mutual interference between parasites or predators and its effects on searching efficiency. J. Anim. Ecol. 44, 331–340. Brose, U., Williams, R.J., Martinez, N.D., 2006. Allometric scaling enhances stability in complex food webs. Ecol. Lett. 9, 1228–1236. Brown, J.H., Gillooly, J.F., Allen, A.P., Savage, V.M., West, G.B., 2004. Toward a metabolic theory of ecology. Ecology 85, 1771–1789.

87

Farrell, A.P., Hinch, S.G., Cooke, S.J., Patterson, D.A., Crossin, G.T., Lapointe, M., Mathes, M.T., 2008. Pacific salmon in hot water: applying aerobic scope models and biotelemetry to predict the success of spawning migrations. Physiolog. Biochem. Zool. 81, 697–709. Gang, H., Ditzinger, T., Ning, C.Z., Haken, H., 1993. Stochastic resonance without external periodic force. Phys. Rev. Lett. 71, 807–810. Guill, C., Drossel, B., Just, W., Carmack, E., 2011a. A three-species model explaining cyclic dominance of Pacific salmon. J. Theoret. Biol. 276, 16–21. Guill, C., Reichardt, B., Drossel, B., Just, W., 2011b. Coexisting patterns of population oscillations: the degenerate Neimark–Sacker bifurcation as a generic mechanism. Phys. Rev. E 83, 021910. Guill, C., Drossel, B., Carmack, E. Cyclic Dominance of Sockeye Salmon can be Obtained from a Predator–Prey Model, unpublished. Gurney, W.S.C., Nisbet, R.M., 1985. Fluctuation periodicity, generation separation, and the expression of larval competition. Theoret. Populat. Biol. 28, 150–180. Hume, J., Shortreed, K., Whitehouse, T., Sockeye fry, smolt, and nursery lake monitoring of Quesnel and Shuswap lakes in 2004, 2005. Available at /http://www.unbc.ca/qrrc/historical_research.htmlS. Heckmann, L., Drossel, B., Brose, U., Guill, C., 2012. Interactive effects of body-size structure and adaptive foraging on food-web stability. Ecol. Lett. 15, 243–250. Hocking, M.D., Reynolds, J.D., 2011. Impacts of salmon on riparian plant diversity. Science 331, 1609–1612. Holling, C.S., 1998. Temperate forest insect outbreaks tropical deforestation and migratory birds. Mem. Entomolog. Soc. Can. 120, 21–32. Johnston, N.T., MacIsaac, E.A., Tschaplinski, P.J., Hall, K.J., 2004. Effects of the abundance of spawning sockeye salmon Oncorhynchus nerka on nutrients and algal biomass in forested streams. Can. J. Fisher. Aquat. Sci. 61, 384–403. Kartascheff, B., Heckmann, L., Drossel, B., Guill, C., 2010. Why allometric scaling enhances stability in food web models. Theoret. Ecol. 3, 195–208. Larkin, G.A., Slaney, P.A., 1997. Implications of trends in marine-derived nutrient influx to south coastal British Columbia salmonoid production. Fisheries 22, 16–24. Levy, D.A., Wood, C.C., 1992. Review of proposed mechanisms for sockeye salmon population cycles in the Fraser River. Bull. Math. Biol. 54, 241–261. Ludwig, D., Jones, D.D., Holling, C.S., 1978. Qualitative analysis of insect outbreak systems spruce budworm and forest. J. Anim. Ecol. 47, 315–332. Murdoch, W.W., Kendall, B.E., Nisbet, R.M., Briggs, C.J., McCauley, E., Boiser, R., 2002. Single-species models for many-species food webs. Nature 417, 541–543. Myers, R.A., Mertz, G., Bridson, J.M., Bradford, M.J., 1998. Simple dynamics underlie sockeye salmon Oncorhynchus nerka cycles. Can. J. Fisher. Aquat. Sci. 55, 2355–2364. Pestal, G., Ryall, P., Cass, A., Collaborative Development of Escapement Strategies for Fraser River Sockeye: Summary Report 2003–2008. Appendix 8. Technical Report. Canadian Department of Fisheries and Oceans, 2008. Ricker, W.E., 1950. Cycle dominance among the Fraser sockeye. Ecology 31, 6–26. Ricker, W.E., 1997. Cycles of abundance among Fraser River sockeye salmon (Oncorhynchus nerka). Can. J. Fisher. Aquat. Sci. 54, 950–968. Rounsefell, G.A., Kelez, G.B., 1938. The salmon and salmon fisheries of Swiftsure Bank Puget Sound, and the Fraser River. Bull. Bur. Fisher. 49, 692–823. Townsend, C.R., 1989. Population cycles in freshwater fish. J. Fish Biol. 35 (Supplement A), 125–131. Walters, C., Christensen, V., 2007. Adding realism to foraging arena predictions of trophic flow rates in Ecosim ecosystem models: shared foraging arenas and bout feeding. Ecolog. Model. 209, 342–350. Worden, L., Botsford, L.W., Hastings, A., Holland, M.D., 2010. Frequency responses of age-structured populations: pacific salmon as an example. Theoret. Populat. Biol. 78, 239–249. Yodzis, P., Innes, S., 1992. Body size and consumer-resource dynamics. Am. Natural. 139, 1151–1175.