# Irregular timing of population cycles

Jonathan E Rubin, David J Hearn, Priscilla E Greenwood, Todd L. Parsons, Karen C Abbott

<figure><img src="image_1.png" /></figure>

## To cite this version:

Jonathan E Rubin, David J Hearn, Priscilla E Greenwood, Todd L. Parsons, Karen C Abbott. Irregular timing of population cycles. Oikos, In press, pp.e09290. 10.1111/oik.09290. hal-03445622v2

HAL Id: hal-03445622

https://cnrs.hal.science/hal-03445622v2

Submitted on 28 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Research article

# Irregular population cycles driven by environmental stochasticity and saddle crawlbys

Jonathan E. Rubin, David J. D. Earn, Priscilla E. Greenwood, Todd L. Parsons and Karen C. Abbott

J. E. Rubin (https://orcid.org/0000-0002-1513-1551) ✉ (jonrubin@pitt.edu), Dept of Mathematics, Univ. of Pittsburgh, Pittsburgh, PA, USA. – D. J. D. Earn (https://orcid.org/0000-0002-7562-1341), Dept of Mathematics & Statistics, McMaster Univ., Hamilton, ON, Canada. – P. E. Greenwood, Dept of Mathematics, Univ. of British Columbia, Vancouver, BC, Canada. – T. L. Parsons, Laboratoire de Probabilités, Statistique et Modélisation (UMR 8001), CNRS & Sorbonne Univ., Paris, France. – K. C. Abbott, Dept of Biology, Case Western Reserve Univ., Cleveland, OH, USA.

Oikos

2022: e09290

doi: 10.1111/oik.09290

Subject Editor: Meike Wittmann

Editor-in-Chief: Dries Bonte

Accepted 15 September 2022

Despite considerable study of population cycles, the striking variability of cycle periods in many cyclic populations has received relatively little attention. Mathematical models of cyclic population dynamics have historically exhibited much greater regularity in cycle periods than many real populations, even when accounting for environmental stochasticity. We contend, however, that the recent focus on understanding the impact of long, transient but recurrent epochs within population oscillations points the way to a previously unrecognized means by which environmental stochasticity can create cycle period variation. Specifically, consumer-resource cycles that bring the populations near a saddle point (a combination of population sizes toward which the populations tend, before eventually transitioning to substantially different levels) may be subject to a slow passage effect that has been dubbed a ‘saddle crawlby’. In this study, we illustrate how stochasticity that generates variability in how close predator and prey populations come to saddles can result in substantial variability in the durations of crawlbys and, as a result, in the periods of population cycles. Our work suggests a new mechanistic hypothesis to explain an important factor in the irregular timing of population cycles and provides a basis for understanding when environmental stochasticity is, and is not, expected to generate cyclic dynamics with variability across periods.

Keywords: periodic orbits, predator–prey dynamics, Rosenzweig–MacArthur, saddle points

## Introduction

Recurrent population cycles – eruptive increases and crashes in population size – are widespread (Kendall et al. 1999, Barraquand et al. 2017). Consumer-resource and host–enemy interactions are famously known to promote cycles (Turchin and Hanski 2001, Myers 2018), but models of these interactions predict a far more regular cycle period than is observed in many real cyclic populations (Dwyer et al. 2004). Ecological data are famously noisy, so the observation that cycle periods appear irregular may seem unsurprising and perhaps even incidental. However, our inability to

<figure><img src="image_3.png" /></figure>
explain the irregularity could reflect a fundamental gap in our understanding of the mechanisms driving cycles. This possibility is troubling because we rely on our understanding of population cycles in pest management, conservation and other important contexts (Barraquand et al. 2017). Proposed mechanisms acting on a population may be accepted or rejected based on whether they predict realistic cycle period variability (Dwyer et al. 2004). More broadly, the variability in a quantity can be as informative as its mean, and repeated calls have been made for improving our ability to extract understanding from patterns in variance (Benedetti-Cecchi 2003, Violle et al. 2012, Holyoak and Wetzel 2020, Shoemaker et al. 2020). The development of general theory for the origins of cycle period variation thus would deepen our understanding of the factors influencing fluctuating populations and improve our ability to manage and protect them.

While it is tempting to implicate environmental stochasticity as a source of cycle period irregularity, the intrinsic period of ecological cycles appears to be surprisingly robust to noise. Past theoretical studies have shown that although white noise readily generates variance in cycle amplitudes, it may have little effect on the regularity of the cycle period. Specifically, Fig. 1a in Dwyer et al. (2004) shows the very regular outbreaks predicted by a stochastic difference equation model for gypsy moth dynamics, and Fig. 9 in Greenman and Pasour (2011) similarly shows virtually no effect of noise on the interepidemic period in a noisy two-strain susceptible-infected model of dengue. As such, other mechanisms, such as chaos, stochastic transitions between alternative attractors, eco-evolutionary feedbacks and resonance between intrinsic cycles and seasonal forcing, have been invoked to explain the irregular timing of some population cycles (Dwyer et al. 2004, Nguyen and Rohani 2007, Ives et al. 2008, Greenman and Pasour 2011, Benincà et al. 2015, Bengfort et al. 2017). These alternative ideas are convincing and no doubt play an important role in some ecological systems. Still, the ubiquity of environmental stochasticity makes it appealing as a parsimonious and perhaps more general explanation for other, unexplained instances of cycle period variability. Though we know from counter-examples (i.e. the aforementioned cases shown in Dwyer et al. 2004, Greenman and Pasour 2011) that environmental stochasticity does not necessarily cause significant variance in cycle period, we do not know that it cannot cause such variance, so we contend that it deserves a closer look.

In this paper, we reconsider the possibility that the irregular timing of population cycles can originate from the interaction between density-dependent consumer-resource interactions and certain forms of parametric noise. In addition to being prone to limit cycles, consumer-resource models generally share another feature: the presence of saddle points at the origin and at other states, such as at the resource carrying capacity in the absence of consumers (Murdoch et al. 2003). The term 'saddle point' in this context refers to a population equilibrium that exists but is unstable in a specific way: from certain starting values,

population levels will tend to approach the saddle point, but eventually they will deviate away from these equilibrium values. If the consumer-resource cycle brings the population trajectory close enough to these saddle points, then a long but ultimately transient epoch within each cycle during which the populations remain near the saddle will result (for an introduction to this role of saddle points, Abbott 2020); this effect has been termed a saddle crawlby (Hastings et al. 2018). The duration of the crawlby scales with the inverse of the distance from the saddle (Cushing et al. 1998, Jäger et al. 2008, Morozov et al. 2020). Therefore, we conjecture that stochasticity that results in more variability from cycle to cycle in how close the populations get to these saddles will produce higher variability in cycle timing. To explore this possibility, we will partition variance in the total duration of a cycle into variances in the durations of its constituent parts, with particular interest in the times spent passing by and moving between the saddles.

Our contributions come from acknowledging two important features of ecological dynamics. The first is the potential role of saddle crawlbys in cycle period variation. Second, a frustrating problem in population ecology is that there are many different aspects of an underlying ecological process that may be influenced by stochasticity, but different stochastic models can lead to very different dynamics (Allen and Allen 2003, Nisbet and Gurney 2003, Coulson et al. 2004, Vadillo 2019). In this work, we use a comparison of two parameters to highlight the fact that environmental stochasticity in different demographic parameters causes different levels of cycle period variability, and we discuss how these parameters' effects on certain features of the cycling dynamics leads to their impacts on variability.

This paper presents our analysis of the roles of saddle crawlbys and different sources of environmental noise in generating variably timed population cycles, based on a series of complementary approaches. We begin with a brief characterization of cycle period variation in a large collection of ecological time series and assess the plausibility of saddle crawlbys as a source of variability in some of these instances. We then present the general theory that links saddle crawlbys to variation in cycle duration in the presence of environmental stochasticity. Specifically, the mechanism involved is that stochastic parameter variation causes variation across successive population cycles in the population densities that arise at stages where saddle crawlbys occur, which in turn induces variation in crawlby durations and overall cycle lengths. The bulk of our paper is devoted to a detailed numerical study of the stochastic Rosenzweig-MacArthur predator-prey model (Rosenzweig and MacArthur 1963), through which we can understand how these theoretical effects manifest in an ecological setting. We close by zooming back out to consider what our numerical study reveals more generally. Our study of the effect of saddle crawlbys on cycle period leads both to a new proposed explanation for some of the irregularly timed cycles observed in nature and to new insights about why past work has seen little effect of environmental stochasticity on cycle period variability.
# Data and analysis

## Data show evidence of cycle period variability and putative saddle crawlbys

While a complete meta-analysis of cycle period variation is well beyond the scope of this study, we do need a reasonable expectation for what levels of variation are realistic. Following Dwyer et al. (2004), we measure cycle period variability by its coefficient of variation (CV), or standard deviation divided by mean. Dwyer et al. (2004) report cycle period CVs ranging from 0.16 to 0.67 for various outbreaking forest insects. We supplement these data with our own analysis of ecological time series from the Global Population Dynamics Database (GPDD; Prendergast et al. 2010). We defined the duration of a cycle as the time from one local maximum of the population time series to the next, and we used wavelet analysis to assess the periodicity of these cycles. Wavelet analysis decomposes a time series into a set of periodic oscillations of different frequencies in order to, among other things, identify the dominant frequencies in the data (see Cazelles et al. 2008, for further information). Beginning with the 664 GPDD time series analyzed by Sibly et al. (2007), we found that 459 had at least 2 cycles and showed some evidence of periodicity (the peak in the true wavelet power spectrum – indicating the strength of the dominant frequency – exceeded the average observed in 100 bootstrap trials). Of these 459, 47 (10.2%) showed perfectly regular cycle periods (CV = 0). The remaining cycle period CVs ranged from 0.10 to 0.98. Overall, the cyclic GPDD data had a mean cycle period CV of 0.30. This analysis confirms that the cycle period CVs reported for forest insects in Dwyer et al. (2004) are representative of other taxa as well.

Identifying a saddle crawlby in data with reasonable certainty requires a carefully validated dynamical model (Hastings et al. 2018) and while this is infeasible in bulk, we can survey the GPDD time series for dynamics that are consistent with our hypothesized crawlby effect. Because the saddle points in common ecological models often occur where at least one species is extinct, we focus on dynamics at low abundances. Populations that become small and stay small for an extended period of time may be crawling by a saddle that lies at zero abundance. Defining ‘low abundance’ as ≤ 5% of the maximum observed abundance (though note that the patterns we show are not sensitive to the exact percentage used; see the Supporting information), we recorded the longest continuous stretch of time that each GPDD time series dropped to low abundance and found this to be positively related to cycle period CV (Fig. 1A, B). Time series with only short (or no) time at low density showed examples of both high and low CVs (blue histogram in Fig. 1B; examples in C and D), suggesting that cycle period variability in these cases is due to a range of mechanisms not necessarily related to saddle crawlbys. (The fact that the relationship in Fig. 1A, though significant, explains only a small amount of the variance in cycle period CV likewise suggests that other mechanisms are represented in these data.) However,

for cases in which durations at low abundance exceeded a threshold (three time series time steps), CVs were relatively large (Fig. 1B), consistent with what we would predict if these populations were influenced by saddle crawlbys at their troughs. Some of these putative crawlbys are almost certainly long transients caused by other phenomena (e.g. human impacts on the fisher (Aubry and Lewis 2003), Fig. 1E), while others are conceivably due to the type of tightly-coupled consumer-resource interactions we study here (e.g. larch budmoth cycles (Turchin et al. 2003), Fig. 1F). In sum, this coarse analysis of 459 ecological time series supports the plausibility of the crawlby mechanism as one among several reasons that real cycle periods vary, and provides a benchmark against which to compare our theoretical results.

## Saddle crawlbys convert variable population levels at specific cycle phases into variable cycle periods

Mathematically, a saddle point is an unstable equilibrium state (e.g. combination of population densities) that, though ultimately repelling, is approached from some states. Ecologically, the influence of a saddle can be seen in the approach toward joint extinction before eventual resource recovery during consumer-resource cycles. Such behavior also appears in the approach from jointly low densities toward the carrying capacity of a resource or fast-growing competitor species, before the successful recovery or reinvasion of a predator or second competitor. A saddle crawlby is dominated by a phase of ‘lingering’ near the saddle point during relatively slow passage toward and then away from it (Hastings et al. 2018), within each population cycle.

In a dynamical system without added noise, the amount of time the system spends near a saddle – that is, the duration of the saddle crawlby – is well understood (Morozov et al. 2020). Using a predator-prey system as an example, we can characterize the crawlby associated with the joint extinction state by considering how the populations enter and exit the proximity of the extinction saddle point at (0, 0). Using $x$ to denote the prey population density and $y$ the predator population density, we can consider the crawlby to start once predator density drops below some threshold, $y_{0,in}$, as illustrated in Fig. 2A. Use $x_{0,in}$ to denote the prey density when the $y=y_{0,in}$ threshold is crossed. In Fig. 2A, we show the crawlby trajectories for two different choices of $x_{0,in}$, one in blue and one in red. We will say that the (0, 0) crawlby is complete once prey density grows from $x_{0,in}$ to a prey threshold value, $x_{0,out}$, also shown in Fig. 2A.

During the crawlby, the populations are sufficiently close to (0, 0) that typical consumer-resource dynamics will be well approximated by a linear system of equations, such as

$$ \frac{dx}{dt} = rx \qquad (1) $$

$$ \frac{dy}{dt} = -my $$
<figure><img src="image_4.png" /><figcaption>(C) Atlantic herring (*Clupea harengus*)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(E) fisher (*Martes pennanti*)</figcaption></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>(D) wood mouse (*Apodemus sylvaticus*)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(F) larch budmoth (*Zeiraphera diniana*)</figcaption></figure>

<figure><img src="image_9.png" /></figure>

Figure 1. Summary of Global Population Dynamics Database (GPDD) time series. (A) Relationship between the longest consecutive timespan spent at low (≤ 5% maximum) abundance and the coefficient of variation (CV) in cycle period, with regression line. Each dot is one of the 459 examples analyzed; colored and labeled dots denote examples shown in panels (C)–(F). (B) Histograms of cycle period CV for time series with short (0–3 time step) visits to low abundance (blue bars, summarized by blue boxplot above) versus those with longer intervals at low abundance (red). (C–F) Examples of time series from the GPDD spanning the range shown in (A). The extended time at low abundance in (E)–(F) is consistent with a crawlby of a saddle point at zero abundance.

with *r*, *m* > 0. System (1) has solutions $x(t) = x(0)e^{rt}$, $y(t) = y(0)e^{-mt}$. Taking the initial conditions as the entry point to the crawlby, $(x(0), y(0)) = (x_{0,in}, y_{0,in})$, and defining a stopping time $t_{pass}$ by the condition $(x(t_{pass}), y(t_{pass})) = (x_{0,out}, y_{0,out})$, we obtain $x_{0,out} = x_{0,in}e^{rt_{pass}}$. This expression yields $t_{pass} = (1/r) \ln(x_{0,out}/x_{0,in}) = (1/r)[\ln(x_{0,out}) - \ln(x_{0,in})]$, which is positive because $x_{0,out} > x_{0,in}$. That is, the crawlby duration $t_{pass}$ depends linearly on $\ln(x_{0,in})$. Similarly, if we substitute this result into

the equation $y_{0,out} = y_{0,in}e^{-mt_{pass}}$, then we find that the logarithm of the predator density at the end of the crawlby also depends linearly on $\ln(x_{0,in})$.

Although we use the (0, 0) crawlby as an illustrative example in making these arguments, the same mathematical features arise for a crawlby past any other saddle point with exponential growth and decay rates in any two-dimensional system. That is, in any such two-dimensional crawlby, there
<figure><img src="image_10.png" /><figcaption>Figure 2. (A) Illustration of a crawlby past a saddle at (0, 0). We define the vicinity of the saddle (gray region) by an entry condition, $y=y_{0,in}$, and an exit condition, $x=x_{0,out}$. The crawlby begins when predator density $y$ crosses the $y_{0,in}$ threshold. Prey density at this point ($x_{0,in}$) is expected to be variable, according to some distribution (black overlain curve). When prey density upon entry is lower (red trajectory, compared to teal), the saddle crawlby will be slower, with a longer duration, because the trajectory passes closer to the saddle. (B) Components of a typical predator–prey cycle, A: passage by the origin when both species are rare, B: recovery of the prey population, C: passage by the prey’s carrying capacity, K, D: predator recovery and prey decline, E: initiation of predator decline and finally F: sharp predator decline back toward the origin. Segments A and C involve saddle crawlbys. Panel (A) is general, but can be thought of as a zoom-in of region A (where the crawlby near the origin occurs) in panel (B).</figcaption></figure>

is contraction in the deviation of one variable, say $y$, from its value at the saddle point while there is growth in the deviation of the other variable, say $x$; moreover, both the crawlby duration and the deviation of the logarithm of $y$ at the end of the crawlby depend logarithmically on the deviation of $x$. A similar scaling of passage time and population sizes results from any crawlby near a saddle point in a system with more than two variables as well.

All of this is true regardless of whether the saddle crawlby is part of a larger, repeating cycle, but if the populations are cycling (Fig. 2B), then we can carry these arguments further. The linear dependence of the crawlby duration $t_{pass}$ on $\ln(x_{0,in})$ implies that if there is variability in the $x_{0,in}$ values from cycle to cycle (as illustrated in the black $x_{0,in}$ distribution curve in Fig. 2A, for example), then there will also be variability in crawlby durations. Intuitively, the most intensive variability in durations will result from an $x_{0,in}$ distribution that features close approaches to the saddle (low $x_{0,in}$) and correspondingly long passage times, along with enough variance to yield some more distant approaches and short passage times. As an example, if the values of $x_{0,in}$ are lognormally distributed with parameters $\mu$ and $\sigma^2$, then the crawlby times from $y=y_{0,in}$ to $x=x_{0,out}$ are normally distributed. Explicitly, if we choose the same fixed value, call it $c$, as both $y_{0,in}$ and $x_{0,out}$ for convenience, so that the grey shaded crawlby region in Fig. 2A is a square, then by applying the formula for $t_{pass}$ to the whole distribution of $x_{0,in}$ values, we find that $t_{pass}$ is distributed normally; specifically, $t_{pass} \sim N(\ln(c)/r - \mu/r, (\sigma/r)^2)$ (we use $N(\mu, \sigma^2)$ to denote the normal distribution with mean $\mu$ and variance $\sigma^2$), with larger values in the $x_{0,in}$ distribution leading on average to smaller passage times $t_{pass}$. Similarly, the distribution of predator densities at the end of the crawlby is log-normal, $\ln(y_{0,out}) \sim N((1 - m/r)\ln(c) + m\mu/r, (m\sigma/r)^2)$, with smaller

$x_{0,in}$ leading on average to smaller $y_{0,out}$. (Note that because $c$ is constant, the way that it shows up in these distributions is not important for understanding how population size variance at the start of the crawlby generates variance in crawlby duration.) After completing a saddle crawlby (i.e. crossing the $x=x_{0,out}$ boundary), the populations continue around the cycle, during which various features, including other crawlbys (Fig. 2B, region C), may occur en route to the next entry through the $y=y_{0,in}$ boundary. Importantly, parameters that do not explicitly affect the saddle crawlby itself can still impact the variability in passage times by impacting the trajectory around the cycle and hence the distribution of entry positions (here, $x_{0,in}$ values) into the saddle crawlby.

In summary, with log-normally distributed prey densities upon entry, the passage times past the saddle at the origin will be normally distributed, the exit predator density will be log-normally distributed, and the variance of both of these distributions will grow with that of the prey distribution upon entry to the neighborhood of the origin. Other types of prey density entry distributions will also translate into predator density exit distributions of the same type. Understanding the stochastic processes that shape the distribution of prey densities at the start of each crawlby is therefore an important goal of our detailed numerical study below.

## Numerical case study: saddle crawlbys generate variable cycle periods in a consumer-resource model

### Formulation and behavior of the underlying deterministic model

To explore the implications of this theory in an ecological context, we analyze the cycles in a stochastic version of the
Rosenzweig–MacArthur predator–prey model. The underlying deterministic model is

$$ \frac{dx}{dt} = rx \left( 1 - \frac{x}{K} \right) - \frac{axy}{x+h} \qquad (2) $$

$$ \frac{dy}{dt} = \frac{abxy}{x+h} - my $$

where *x* is prey population density and *y* is predator population density. The prey population grows at maximum rate *r* and would equilibrate to carrying capacity *K* in the absence of predation. Predation occurs according to a type 2 functional response (Holling 1965) with maximum attack rate *a* and half-saturation constant *h*. Predator population growth occurs at a rate proportional (described by *b*) to its prey consumption rate, and predator death is density independent with rate *m*.

Our baseline parameter values are listed in Table 1. These values were selected to provide certain dynamical features and do not represent a specific biological scenario. Similar results hold for a wide range of other parameter choices that give rise to similar dynamical features, as discussed below. The implied time unit is most easily interpreted via the mean predator lifetime, which is $1/m=1.67$ time units. These values result in a limit cycle featuring close, but not too close, passage to the $x=0$ and $y=0$ axes and thus to the saddle points of system (2) at $(x, y) = (0, 0)$ and $(K, 0)$. Close passage sets the stage for saddle crawlbys, but passage too close to the axes results in unrealistic dynamics, since models such as (2) allow recovery from infinitesimally small population densities at which real populations would almost certainly go extinct. We used $10^{-5}$ times the prey carrying capacity as the extinction threshold, and chose our baseline parameter set so that this threshold was not crossed by the deterministic model with reasonable initial conditions, and crossed only rarely in stochastic simulations unless the noise amplitude was quite high. We omit from our analyses any simulated cycles in which prey and/or predator dropped below the extinction threshold. Imposing this threshold ensures that crawlby effects that we observe in the model are not merely artifacts of letting the system drop to unrealistically low population sizes.

For stochastic *r* or *h*, the values listed refer to $\bar{r}$ or $\bar{h}$, respectively in the Cox–Ingersoll–Ross (CIR) process (3).

To study saddle crawlbys, we should first characterize when system (2) generates oscillations and how parameters impact the paths of oscillations, when present, past saddle points. Depending on the parameter values, the deterministic Rosenzweig–MacArthur model can exhibit a stable equilibrium corresponding to co-existence of the two species or a stable limit cycle. The system undergoes a supercritical Andronov–Hopf (AH) bifurcation at $h=K(ab-m)/(ab+m)$, and cycles arise for *h* below the bifurcation value; note that we only consider parameter values such that $ab > m$. As long as we restrict to parameter values with *h* below the bifurcation value and $ab > m$, whether we choose our baseline parameter set or some other values, system (2) will exhibit periodic orbits, and varying any parameter value individually will affect the amplitudes of these orbits and their proximity to the x- and y-axes, as shown for the parameter *h* in Fig. 3A and B.

The interaction of oscillations with saddle points can be studied by focusing on the parameters *r* and *h*. Because the bifurcation condition does not involve the parameter *r*, changing *r* cannot induce or eliminate cycles in system (2). Nevertheless, when cycles exist, lowering *r* lengthens the oscillation period by slowing the recovery of the prey population; the periods arising for a range of *r* values are displayed as the solid black curve in Fig. 3C. This effect prolongs the crawlbys past both the extinction point (0, 0) and the carrying capacity (*K*, 0) (Supporting information provides the boundary conditions used to define the regions, while the actual crawlby durations appear as the red dashed and blue dotted lines, respectively, in Fig. 3C). However, as *r* decreases, the slowing of the crawlbys is matched by comparable slowing of other parts of the cycle, so the fraction of time the system spends in crawlbys does not increase significantly (Fig. 3E, solid black line).

As a contrasting example, we see that reducing *h* can induce an AH bifurcation and subsequently drive the system deeper into the cyclic regime. Decreasing the value of *h* pushes trajectories closer to both axes (Fig. 3A, B). This prolongs the two saddle crawlbys (Fig. 3D, red and blue lines) by placing the populations closer to the saddles. The saddle

<table><caption>Table 1. Baseline parameters of the model defined by Eq. 2 and 3. For stochastic *r* or *h*, the values listed for *r* and *h* are used as the values for $\bar{r}$ and $\bar{h}$, respectively, in the Cox–Ingersoll–Ross (CIR) process (3).</caption><thead><tr><th>Parameter</th><th>Meaning</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>r</td><td>Prey intrinsic growth rate</td><td>1</td><td>time<sup>-1</sup> per capita</td></tr><tr><td>K</td><td>Prey carrying capacity</td><td>1</td><td>Prey biomass</td></tr><tr><td>a</td><td>Maximum prey consumption rate</td><td>2</td><td>time<sup>-1</sup> per predator</td></tr><tr><td>h</td><td>Prey density at half maximum consumption rate</td><td>0.15</td><td>Prey</td></tr><tr><td>b</td><td>Prey to predator conversion efficiency</td><td>0.5</td><td>—</td></tr><tr><td>m</td><td>Predator death rate</td><td>0.6</td><td>time<sup>-1</sup> per capita</td></tr><tr><td>p</td><td>Parameter subject to CIR process (r or h)</td><td></td><td></td></tr><tr><td>γ</td><td>Speed of adjustment to mean p̄</td><td>1</td><td>time<sup>-1</sup></td></tr><tr><td>σ</td><td>σ√p̄/2γ is standard deviation of the CIR process</td><td>0–0.3</td><td>√units of p</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>Figure 3. Oscillations in system (2) without noise. (A) The system exhibits stable oscillations for $h \lesssim 0.25$ and baseline values of the other parameters. For smaller values of $h$, periodic orbits pass closer to the coordinate axes. Shapes denote points 1 time unit apart on each orbit and thus show the rates of passage through different parts of the cycle. (B) As $h$ is reduced from 0.4 toward its baseline value (0.15), system (2) undergoes an Andronov-Hopf (AH, black circle) bifurcation at which the coexistence equilibrium switches from stable (solid black) to unstable (dashed black) and a family of stable periodic orbits emerges. The black curve with the blue side branches shows the $x$ value at the equilibrium point for each $h$, while the blue branches show the maximal and minimal values of $x$ along the periodic orbits. The other black curve and its red branches are analogous for $y$. (C) As $r$ (baseline value 1.0) decreases, the oscillation period grows (solid black line), partially due to prolonged crawlbys past the saddles at (0, 0) (blue dotted line) and $(K, 0)$ (red dashed line). The purple dashed line shows the total amount of time during each cycle that is spent in a saddle crawlby (red plus blue). (D) As $h$ decreases, the oscillation period grows, but less extremely than in (C). The longer cycle duration is even more strongly influenced by saddle crawlbys; for low $h$, most of the time needed to complete a cycle is spent crawling by the saddle points. The bifurcation diagram and periods displayed in this figure were computed using XPPAUT (Ermentrout 2002). (E-F) The dotted (upper) and solid (lower) black lines show the time spent passing through parts of the cycle other than crawlbys and the proportion of the cycle spent in crawlbys, respectively, as $r$ (E) and $h$ (F) are varied.</figcaption></figure>

crawlbys are not only slower with lower $h$, but they also take up proportionately more of the time needed to complete a cycle (compare the solid black curve showing overall cycle periods to the dashed purple curves in Fig. 3D).

In sum, slowing the overall dynamics (by decreasing the prey's population growth rate, $r$) and pushing trajectories closer to the two saddle points (by decreasing the predator's half-saturation constant, $h$) will each extend the cycle period, but they do so through different mechanisms. Decreasing the prey's growth rate extends the duration of saddle crawlbys, but only by slowing the dynamics as a whole and prolonging

the entire cycle period. Decreasing the half-saturation constant enhances the degree to which predators can overexploit their prey, leading to more pronounced population crashes. These crashes exaggerate the crawlby effect by placing the populations closer to the two saddle points at joint extinction and at the prey carrying capacity. The result is longer cycles that are made up of relatively more time spent for prey recovery (passage near (0, 0)) and then predator recovery (passage near $(K, 0)$). These deterministic results hint at how environmental stochasticity that causes variability in specific model parameters could modulate the influence of saddle crawlbys
and hence, as discussed in the ‘Saddle crawllbs convert variable population levels at specific cycle phases into variable cycle periods’ section, might impact cycle durations. With this insight attained, we are now ready to explore model behavior when stochastic effects are incorporated.

## Including environmental stochasticity

To model different possible ways for environmental stochasticity to affect the predator–prey interaction, we replace individual parameters with a Cox–Ingersoll–Ross (CIR) stochastic process (Cox et al. 1985). That is, we can make any parameter *p* (where *p* is a stand-in for *r*, *h*, ...) stochastic by modeling it as

$$dp = \gamma(\bar{p} - p)dt + \sigma\sqrt{p}dW \quad (3)$$

where *W(t)* is a Wiener process (i.e. a continuous stochastic process for which successive time increments, *W(t+Δt) – W(t)*, are independently normally distributed, with mean 0 and variance Δt) and $\bar{p}$ is the mean (i.e. baseline) value for the parameter. This formulation prevents parameters from becoming negative and has been shown in past work to have various mathematically desirable and biologically reasonable properties: it is mean-reverting and has a positive temporal autocorrelation that decays exponentially (Allen 2016).

All of our numerical results on the stochastic dynamics of our model represent averages derived from simulations run for 20 000 time units. For each run, we discard an initial period extending from time 0 up until asymptotic cycling behavior has been reached (defined as the time of the first crossing of one of the boundaries given in the Supporting information). Each sample on which averaging is performed includes at least 400 cycles during which population sizes stay above the extinction threshold. We have also used small-noise approximations to derive analytical relationships between certain model features and passage times along the axes and near the saddle points. Although the details of these calculations are outside the scope of this article, we provide a summary of relevant results in the Supporting information.

We focus here on contrasting the effects of a noisy prey growth rate, *r*, versus a noisy predator half-saturation constant, *h*. This pair of parameters spans many interesting contrasts (i.e. appearance in a linear versus only in a nonlinear term, a multiplier versus a denominator, a rate in one species' equation versus both), represents distinct modes of influencing cycle period in the deterministic model ('Formulation and behavior of the underlying deterministic model' section), and thus serves well to illustrate the impacts of different sources of noise. We recognize that natural species interactions are not characterized by a single stochastic parameter, but modeling one parameter at a time as stochastic allows us to isolate different putative effects of noise. Modeling noise only in *r* approximates the situation where the prey's intrinsic growth rate is more sensitive to environmental perturbations (e.g. fluctuations in temperature) than are other aspects of the system. Noise in *h* represents a situation where predator

behavior or physiology is the system feature most sensitive to perturbations. Simulations with noise in other parameters that, like *r*, appear linearly in Eq. 2, such as *a*, *b*, *m*, yield similar results to those with noise in *r*, consistent with a mathematical averaging argument (Supporting information).

## Stochasticity can cause variance in cycle period, but not all stochasticity has the same impact

Adding Cox–Ingersoll–Ross parameter noise (Eq. 3) to the Rosenzweig–MacArthur model results in noisy cycles such as the ones displayed in Fig. 4. Figure 4A and B shows an example in which noise in *r* causes cycles to vary in amplitude but maintain a nearly constant period. Figure 4C and D shows that for the same noise amplitude (σ value in Eq. 3), the cycle period is much more variable with noise in *h*.

Figure 4 shows just two examples, but they are representative of the systematic differences we found when we explored different intensities of environmental noise acting on different model parameters. For the same absolute noise magnitudes (σ values), the cycle period becomes much more variable when noise affects the half-saturation constant of the functional response, *h*, than when it affects prey growth rate, *r* (Fig. 5A, B). Specifically, for noise in *r*, the mean cycle duration changed little with increasing σ from the deterministic value of 24.83 time units, and the standard deviation of these durations depended only weakly on σ (Fig. 5A). Conversely, mean cycle durations and variability increased substantially with σ when noise entered through *h* (Fig. 5B), with the distribution developing a tail of long cycle durations. For comparison with the biological data, note that the results plotted here yield CVs from 0.04 to 0.12 with noise in *r* and from 0.11 to 0.33 with noise in *h*. Results for σ=0.2 (red bars in Fig. 5) are also summarized in Table 2 (cases 2–3). Finally, if noise enters into both *h* and *r* simultaneously (Table 2 case 4, implemented via two copies of Eq. 3 that either share (correlated) or have independent (uncorrelated) instantiations of the Wiener process *W*), the results are similar to the case with noise in *h* alone (case 3).

There are several possible explanations for the stronger effect of noise in *h* than noise in *r* in generating variable cycle periods. This is the pattern we expect if cycle period variability is driven by variation in the time it takes to complete saddle crawllbs, because changes in *h* exaggerate the influence of crawllbs while changes in *r* do not (Fig. 3C, D). Alternatively, or perhaps in addition, noise of a particular magnitude may have a larger impact when it enters through *h* simply because our baseline value of *h* is lower. In other words, a noise amplitude of, say, σ=0.2 acting on an $\bar{r}$ of 1.0 might reasonably have less effect than the same σ=0.2 acting on an $\bar{h}$ of 0.15. We tested this idea in two ways. First, we reduced $\bar{r}$ to 0.15 to eliminate the difference in magnitudes between *r* and *h*. This results in unrealistically long cycles, but allows a straightforward comparison of standard deviation and CV. Again, the CV in cycle period is much higher when noise enters through *h* (Table 2, cases 6–7). Second, we returned to the baseline parameter set and increased the noise amplitude
<figure><img src="image_12.png" /><figcaption>Figure 4. Typical behavior of system (2)–(3) with noise amplitude $\sigma=0.3$ and noise in $r$ (A–B) or $h$ (C–D). (A, C) show the dynamics over 2000 time units in the $(x, y)$ phase plane, during which several population cycles, each progressing counterclockwise, occur. (B, D) show time series of $x$ (black) and $y$ (red) for the first 500 of the 2000 time units shown in (A, C). Because noise is present, each cycle is defined from the time of crossing through the entry boundary into region D (Fig. 2B) up to the next such crossing time, under the requirement that all other regions are traversed between these crossings. This definition produces results that largely agree with the peak-to-peak (or trough-to-trough) times but avoid the issue of stochastic fluctuations that produce multiple local peaks and troughs.</figcaption></figure>

on $r$ (by rescaling by $\sqrt{h}$) to allow the prey population to become quite small at times, approaching the x-axis and now allowing the exaggerated crawlby effect to be seen via noise in $r$. The result is a cycle period CV that is now quite comparable to the baseline case with noise in $h$ (Table 2, case 5 compared to case 3).

We thus found evidence of true differences (Table 2, cases 2–3 and 6–7) in how noise in the prey growth rate and noise in the predator's half-saturation constant affect variability in the cycle period. We attribute these differences to the way changes in $h$ affect the proximity of cycles to saddle points and hence the relative speed of saddle crawlbys, and we investigate this hypothesis in more detail below ('Cycle variability is driven by variable passage times near extinction and prey carrying capacity' section). We also found evidence that at some sufficiently high level, noise in either parameter can lead to exaggerated crawlbys that cause cycle periods to become more variable (Table 2, cases 3 and 5).

Finally, we note that the half-saturation constant, $h$, is not unique in its ability to drive the system through the AH bifurcation and produce cycles of increasingly large magnitude, as shown in Fig. 3B; in fact, all parameters except $r$ have that property in this model. Therefore, we consider one

more scenario before moving on: we increase $m$ from its baseline value of 0.6 to 0.65. In the absence of stochasticity, this change moves the periodic orbit of the deterministic system (2) from the outer cycle in Fig. 3A to a path with a less proximal passage past the axes. That is, while the system is still cycling with $m=0.65$, the expected (mean) trajectory does not pass as close to the saddle points. Unexpectedly, changing $m$ in this way increases the CV of cycle periods (Supporting information; Table 2, cases 8–9 compared to 2–3), and the trend continues for $m=0.7$. Once again, this effect is particularly strong for noise in the half-saturation constant (Table 2). Although we might have expected more cycle period variability deeper into the crawlby region (smaller $m$), it appears that being closer to the bifurcation increases the flexibility in the paths that orbits can take. The long durations are still associated with close crawlbys, but crawlbys now occur within a more variable mix of long and short cycles. Recall that CV is given by standard deviation divided by mean. As $m$ increases, if we were to lose the long duration cycles, then the mean and standard deviation could both decrease in such a way that CV could drop. In reality, what we find is that increasing $m$ does introduce some short cycles but also maintains enough long crawlby cycles that mean does not change much while standard deviation grows, and hence CV does increase. These
<figure><img src="image_13.png" /><figcaption>Figure 5. Distributions of relative cycle durations for various noise amplitudes σ. (A–B) Results from simulations with model (2) and (3). (A) With noise in *r*, the mean and median duration remain nearly constant, while the distribution widens and remains symmetric, as σ increases. (B) With noise in *h*, the mean and median duration increase and the distribution widens and becomes skewed, with a tail of long duration cycles, as σ increases. Note that the two panels show different ranges of cycle durations. (C–D) Results from Gillespie simulations of model (2) with carrying capacity *K*=8000 and environmental noise in *r* (C) or in *h* (D) are similar to those from the simulations with only environmental noise shown in (A, B). These plots were derived from equal time simulations across noise conditions, which led to the following numbers of cycles. For the Gillespie case: σ<sub>r</sub>=0.1, 978 cycles; σ<sub>h</sub>=0.2, 915 cycles; σ<sub>v</sub>=0.3, 839 cycles; σ<sub>p</sub>=0.1, 1018 cycles; σ<sub>v</sub>=0.2, 1018 cycles; σ<sub>p</sub>=0.3, 1003 cycles. For simulations of (2), (3): σ<sub>h</sub>=0.1, 775 cycles; σ<sub>v</sub>=0.2, 668 cycles; σ<sub>p</sub>=0.3, 410 cycles; σ<sub>v</sub>=0.1, 809 cycles; σ<sub>p</sub>=0.2, 808 cycles; σ<sub>v</sub>=0.3, 798 cycles.</figcaption></figure>

trends persist if *m* is increased further, until we eventually run
into difficulty defining distinct cycles. The trends are reversed
if *m* is decreased, although as *m* becomes smaller, we start to
observe more cycles on which *x* and *y* drop below our mini-
mum size cutoff of 10<sup>-5</sup>.

***Demographic noise does not overwhelm environmental noise effects***

We have included noise in model (2) as Cox–Ingersoll–Ross
parameter noise (Eq. 3), representing environmental vari-
ability. In reality, population dynamics necessarily include
demographic noise as well. In theory, this variability across
individuals could induce dynamic effects that dominate the
overall cycle paths, reducing the importance of the environ-
mental fluctuations that we have described. We tested this
possibility by performing simulations in which all aspects
of the population dynamics are treated as events with pre-
scribed expected rates, and events are implemented sequen-
tially at random, with event likelihoods and timings based
on these expected rates. The details of these simulations,

which are based on the Gillespie algorithm, are provided in
the Supporting information. We found that these simula-
tions yielded very similar distributions of cycle durations to
our direct simulations of system (2)–(3), matching the differ-
ences between noise in *r* and noise in *h* as well as those across
noise amplitudes *σ* (Fig. 5; Supporting information). This is
important for confirming that none of our main findings are
expected to disappear in the presence of demographic noise,
and that they are not due simply to specific modeling or meth-
odological choices.

***Cycle variability is driven by variable passage times near extinction and prey carrying capacity***

We have shown how parameter changes can affect the prox-
imity of cycles to the saddle points (Fig. 3) and how sto-
chastic parameter variation generates variable cycle periods
(Table 2, Fig. 5). One step that remains is to connect these
results and determine whether cycle period variation can
indeed be traced to variable passage by the saddle points. To
do this, we begin by partitioning the cycle into a sequence of
<table><caption>Table 2. Mean, standard deviation and coefficient of variation of cycle periods under different scenarios.</caption><thead><tr><th>Case</th><th>Noisy parameter(s)</th><th>Parameter set</th><th>Noise amplitude, σ</th><th>Mean period</th><th>Standard deviation</th><th>CV</th><th>Explanation</th></tr></thead><tbody><tr><td>1</td><td>None</td><td>Baseline</td><td>0</td><td>24.83</td><td>0</td><td>0</td><td>Deterministic</td></tr><tr><td>2</td><td>r</td><td>Baseline</td><td>0.2</td><td>24.74</td><td>2.07</td><td>0.084</td><td>Base stochastic case to compare</td></tr><tr><td>3</td><td>h</td><td>Baseline</td><td>0.2</td><td>28.88</td><td>5.59</td><td>0.19</td><td>Base stochastic case to compare</td></tr><tr><td>4</td><td>r and h</td><td>Baseline</td><td>σ<sub>r</sub>=σ<sub>h</sub>=0.2</td><td>28.70</td><td>5.47</td><td>0.19</td><td>Shared Wiener process</td></tr><tr><td>4'</td><td>r and h</td><td>Baseline</td><td>σ<sub>r</sub>=σ<sub>h</sub>=0.2</td><td>29.76</td><td>6.01</td><td>0.20</td><td>Independent Wiener processes</td></tr><tr><td>5</td><td>r</td><td>Baseline</td><td>0.2/√h = 0.516</td><td>25.26</td><td>5.06</td><td>0.20</td><td>Scaled up amplitude</td></tr><tr><td>6</td><td>r</td><td>r&#773; = h = 0.15</td><td>0.2</td><td>68.04</td><td>4.30</td><td>0.063</td><td>Set r&#773; equal to baseline h</td></tr><tr><td>7</td><td>h</td><td>r&#773; = h&#773; = 0.15</td><td>0.2</td><td>61.81</td><td>7.58</td><td>0.123</td><td>Set r equal to baseline h</td></tr><tr><td>8</td><td>r</td><td>m=0.65</td><td>0.2</td><td>23.79</td><td>2.30</td><td>0.097</td><td>Increase m, which decreases baseline orbit amplitude</td></tr><tr><td>9</td><td>h</td><td>m=0.65</td><td>0.2</td><td>28.95</td><td>7.80</td><td>0.27</td><td>Increase m, which decreases baseline orbit amplitude</td></tr></tbody></table>

Results are based on simulations of (2) with our baseline parameter set given in Table 1, except where noted. In all cases except for the deterministic case (case 1), either *r* or *h* (or both in cases 4 and 4') is made stochastic by applying Eq. 3; the Wiener processes in cases 4 and 4' appear in (3).

six segments (Fig. 2B), **A**: passage by (0, 0); **B**: prey growth from near 0 to near *K*; **C**: passage by (*K*, 0); **D**: predator growth and prey decline (that is, movement from near (*K*, 0) toward the *y*-axis); **E**: the transition from predator growth to predator decline (i.e. passage across the predator's nullcline); and finally, **F**: sharp predator decline to again approach (0, 0). Our hypothesis that saddle crawlbys drive cycle period variability leads us to predict that significant variance should be observed in the passage times near (0, 0) and (*K*, 0). We calculate times spent transitioning through each of the six segments as the populations undergo a cycle, based on times when trajectories cross partition entrance and exit boundaries at locations that we specify in the (*x*, *y*) plane. The locations used for these boundaries for most of our simulations are listed in the Supporting information; for a few special cases with altered cycle paths (cases 6–9 in Table 2), we used other boundaries. The boundary locations ensure that fluctuations rarely cause spurious crossings, while even those orbits with relatively large fluctuations do not miss any boundaries. Passage times through partition segments are computed sequentially, with the exit boundary from one region also serving as the entrance boundary for the next region. For example, once a trajectory exits a region through its exit boundary, re-entry through the exit boundary due to stochasticity does not count as another passage through the original region; rather, it is ignored, and the calculation of the time of passage through the next region continues until that region's exit boundary is crossed. This subdivision allows us to decompose the variability in the full cycle period based on the variability in the time needed to complete each step.

The distributions of times spent in these regions differ significantly in their variability. By a clear margin, the most variability arises in the passage times through regions **A** and

**C**, corresponding to crawlbys near the saddles (Supporting information). Moreover, substantially more variability occurred in these regions with noise in *h* than with noise in *r*; note the differences in *x*-axis details between different histograms in Supporting information as well as between the histograms in Fig. 5A and those in Fig. 5B. The Gillespie simulations shown in Fig. 5C and D yielded qualitatively, and in many cases quantitatively, similar results. Supporting information includes a thorough comparison of the durations and standard deviations of passage times through the individual regions, for various noise amplitudes, across noise in *r* and *h* and across the two simulation experiments.

While substantial variability in the durations of the crawlbys (that is, in passage times near (0, 0) and (*K*, 0)) occurs, the times spent within crawlbys remain long enough to dominate the overall cycle period. For example, the positive correlation between period and these passage times, but not between period and passage time along the *x*-axis, is readily apparent in Fig. 6. These relationships, especially for the passage time near (0, 0), are even tighter with noise in *r* than with noise in *h*, although the passage times themselves are less widely distributed when the noise is in *r*.

### Loss of prey strongly influences each cycle period, with little history dependence between cycles

We have shown strong evidence that environmental noise has different impacts on the distribution of cycle periods, depending on the origin of the noise, and that variability in cycle periods is dominated by variability in the duration of the two saddle crawlbys near (0, 0) and (*K*, 0). To conclude our characterization of these effects, we ask what factors best predict these crawlby durations. Longer crawlbys near origin or boundary saddle points elevate the risk of stochastic extinction
<figure><img src="image_14.png" /><figcaption>Figure 6. Correlations between the different durations for (A–B) noise in the prey growth rate (r) or (C–D) noise in the predator's half-saturation constant (h) for baseline parameter values and noise amplitude σ=0.2. (A, C) The full cycle period (per) is tightly correlated with passage times near the (0, 0) and (K, 0) saddle points (blue triangles and red circles, respectively). In contrast, it is effectively independent of the time it takes to transition between the saddles ('x-axis passage', cyan markers). (B, D) Cycle period is highly correlated with the total time needed to complete both saddle crawlbys ('sum', black dots). It is uncorrelated with the time needed to complete all other parts of the cycle ('difference', i.e. the total cycle period minus the sum of the two saddle crawlby durations; green dots).</figcaption></figure>

by prolonging the time during which one or both popula-
tions are small. Longer crawlbys also delay the time until the
next population peak, with important consequences for some
cyclic species like outbreaking forest insects. We therefore ask,
what should we look for in a population's history to help us
anticipate the durations of upcoming saddle crawlbys?

For deterministic dynamics, we saw that the time to com-
plete a crawlby of the extinction state is linearly related to
ln(x₀ₗ), the logarithm of the prey density when the popula-
tions first enter the vicinity of the extinction saddle ('Saddle
crawlbys convert variable population levels at specific cycle
phases into variable cycle periods' section). In the presence
of environmental stochasticity, we find that a similar relation
still holds, with a strong correlation between ln(x₀ₗ) and the
passage time by the (0, 0) saddle (Fig. 7A, B). There is also a
strong correlation between log prey density, ln(x₀ₗ), and log
predator density upon entry into the neighborhood of the (K,
0) saddle (Fig. 7C, D), and between the time of passage near
(K, 0) and this log predator density (Fig. 7E, F). Thus, even in
the presence of noise, the prey's population density near the
end of the predator-decline portion of the cycle, which can be
considered the start of a saddle crawlby, strongly influences

the durations of both saddle crawlbys, which in turn strongly
influence the time it will take to complete the current cycle.

The effects of noise during the passage around the rest
of the cycle eliminates the effects of these entry positions,
so that the positions of entry into the neighborhoods of the
saddles in the next cycle are effectively independent of past
history (Supporting information). This independence means
that we expect no meaningful autocorrelation in the periods
of consecutive cycles. It also validates our decision (Methods)
to omit individual cycles with unrealistically low population
sizes from our analyses. The consistency of Gillespie simula-
tion results with those from simulating Eq. 2 and 3 provides
additional evidence that successive cycles are effectively inde-
pendent. The timescale of the CIR stochastic process, γ, likely
plays a role in this effect and the others that we describe.
Specifically, if γ were too large, then stochastic parameters
would be pinned close to their mean values, reducing vari-
ability, while if γ were too small, then significant history-
dependence could emerge; for a wide range of γ between
these extremes, the results that we report are preserved.

Finally, referring back to Fig. 7 and comparing the two
columns, we notice that the predator and prey densities upon
<figure><img src="image_15.png" /><figcaption>Figure 7. Strong dependence of passage times through saddle neighborhoods on entry positions into those neighborhoods. In (A, C, E) the prey growth rate, *r*, is the stochastic parameter and in (B, D, F) the predator's half-saturation constant, *h*, is stochastic. In all cases, the quantity on the horizontal axis is shown on a logarithmic scale; a logarithmic scale is also used on the vertical axis in (C, D). (A–B) show the linear dependence of the (0, 0) saddle crawlby duration on the log of the prey's density upon entry into that neighborhood, (ln(*x*<sub>0,in</sub>)). (C–D) show linear dependence of log predator density upon entry into the neighborhood of (*K*, 0), ln(*y*<sub>*K*,in</sub>), on ln(*x*<sub>0,in</sub>). Finally, (E–F) show the linear dependence of the duration of the (*K*, 0) saddle crawlby on ln(*y*<sub>*K*,in</sub>). Results for noise amplitude σ=0.1, 0.2 and 0.3 are shown in (A, B, E, F), while (C, D) show σ=0.2.</figcaption></figure>

entry into the two saddle neighborhoods differ, depend-
ing on the source of noise. When noise enters through the
prey's population growth rate, *r*, the entry values seem to
be symmetrically distributed on the logarithmic scale, while
when noise enters through the half-saturation constant, they
become skewed toward smaller values even on the logarith-
mic scale, corresponding to the skew toward longer passage
times (Supporting information). A mathematical argument
for this disparity is provided in the Supporting information.

Generalization

We have demonstrated a range of phenomena related to
cycle duration variability associated with saddle crawlbys in
the Rosenzweig–MacArthur model (2) with environmental
stochasticity of the form (3). There are two natural ways to
think about generalization from this specific case study. First,
from a mathematical perspective, we have already identified
key elements that can give rise to cycle period variability: 1)
the presence of a saddle point; 2) an underlying, deterministic
limit cycle that allows the system to repeatedly pass near the saddle, creating the opportunity for crawlbys; and 3) a source of environmental noise that combines with the deterministic cycle to induce a distribution of entry positions into the neighborhood of the saddle, to generate both close and more distant passages. In any two-dimensional system, near any saddle equilibrium point featuring exponential approaches of population levels toward the equilibrium and growth of population levels away from the equilibrium, the dynamics is qualitatively similar to system (1), with natural generalizations to higher dimensions. As an additional test of this mechanism, we took the set of $x_{0,in}$ values observed in our previous simulations of the Rosenzweig–MacArthur model with $\sigma=0.2$ and used these in the initial conditions of the more general, deterministic system (1). We simulated (1) and recorded the time needed to reach $x=0.1$ from points $(x, y)=(x_{0,in}, 0.1)$, then calculated the CV in these passage times across the different $x_{0,in}$ values. When we obtained our $x_{0,in}$ values from Rosenzweig–MacArthur simulations with noise in $h$, the resulting CV in passage time was 0.62. In contrast, when our $x_{0,in}$ values came from simulations with noise in $r$, the CV was only 0.19. This provides additional confirmation of the central role of the population distribution upon approach to a saddle crawlby in producing variability in cycle durations. An interesting challenge for future work is to determine how different noise processes (i.e. other than Eq. 3) affect the shape of this pre-crawlby population distribution.

A second perspective on generalization from system (2) is that this model represents certain biological features such as logistic resource population growth in the absence of consumption that are present in a wide range of consumer-resource interactions, but also neglects some other features. Thus, we can ask what happens if we incorporate alternative or additional biological effects. We consider this question in the Supporting information. In brief, we find that our results on the link between saddle crawlbys and cycle duration variability do generalize, with other dynamical features making secondary contributions to this variability.

## Discussion

Although real cyclic populations can display considerable variability in cycle period, previous authors have observed that stochastic models of cycling populations typically predict much more regular periods (Dwyer et al. 2004). Inspired by recent work emphasizing the contribution of saddle crawlbys to non-equilibrium ecological dynamics (Hastings et al. 2018), we hypothesized that stochasticity that causes variation in how close populations come to saddle points during their cycles might promote variability in cycle periods. Our numerical study of stochastic Rosenzweig–MacArthur predator–prey cycles confirms this idea. For the cycles that we observe in model simulations, variability in period is dominated by deviation in the passage times past the saddle points (Supporting information); further evidence for the central impact of saddle crawlbys arises in the relation of these passage

times to the proximity of trajectories to the stable manifolds of the two saddles (i.e. the x- and y-axes; Fig. 7). Our simulation results reflect averages over many more cycles than will be present in empirical time series. However, the independence of crawlby durations from one cycle to the next suggests that although there will be a drop in precision with smaller samples, there should not be any systematic bias that would prevent the effects we illustrate here from appearing in empirical data. Our results thus update our understanding of the conditions under which environmental stochasticity can contribute meaningfully to observed cycle period variation in nature.

Our work was motivated by the paradox posed by Dwyer et al.'s (2004) observation that forest insect cycle durations are much more variable than predicted by standard models, and our specific focus on consumer-resource interactions reflects the long history of understanding forest insect dynamics through consumer-resource models (Anderson and May 1980). However, the theoretical basis for our result is general: processes that generate variability in the conditions present when a saddle crawlby begins will generate variability in the duration of whatever dynamic pattern that crawlby is a part of, whether cyclic or otherwise. For example, variable times to invasion by a novel competitor could be due to variability in the initial state of the community or in its state following management activities, if the native community's equilibrium is a saddle (as in standard Lotka–Volterra competition; see Francis et al. (2021)). In short, the situation we chose to study here is just one example of how saddle crawlbys might meaningfully impact the ecological patterns we observe.

The crawlby-induced mechanism that we present has important features in common with some previously-proposed explanations for variable cycle periods. In the gypsy moth outbreaks studied by Dwyer et al. (2004) and the midge outbreaks studied by Ives et al. (2008), for example, variable timing of fluctuations arises from the effect of stochasticity on multiple coexisting attractors, separated by saddle dynamics, in the underlying deterministic 'skeleton', or model of density-dependent interactions on top of which stochasticity introduces randomness. In both these models and ours, the saddles are a feature of the deterministic skeleton, and stochasticity allows for movement past the saddles to require variable passage times. However, the important distinction is that in the attractor-switching mechanism, the transit past the saddle itself is relatively fast (Boettiger and Hastings 2013) and not the source of variability, which instead arises from variable residence times at one attractor or the other. Nevertheless, these examples and ours jointly demonstrate that a broader consideration of the role of saddles and other non-equilibrium phenomena, and of how they interact with stochasticity, is likely to continue deepening our ecological understanding (Hastings et al. 2018, 2021, Abbott 2020). For any particular system with variably timed cycles, fitting models to data and using dynamical systems techniques to study the resulting models for evidence of saddle crawlbys, multiple stable attractors, resonance or other dynamic features would provide the clearest view of which putative mechanism may be acting to induce variability; for example,
models derived from data may have saddles at states where the system tends to linger, suggesting a crawlby mechanism, whereas those from other data may instead have attractors at such states, supporting the attractor-switching mechanism.

In brief, the first ingredient in the crawlby-induced variability mechanism that we describe is that consumer-resource dynamics sets up a log-normal distribution of population sizes at return times to the same neighborhood of a saddle across successive cycles (Supporting information). This distribution of population sizes results in a broad distribution of passage times through the saddle crawlby, and hence of overall cycle durations. Within this qualitatively general phenomenon, quantitative differences arise depending on the source and amplitude of the stochasticity in the dynamics of population interactions, which determines the shape of this log-normal distribution, as well as the extent to which fluctuations are impactful during the saddle crawlby itself (note the differences in passage time variability for fixed entry values in Fig. 7A, B, E, F). We suspect that this general framework is important in a variety of contexts and warrants further study. In our work, we find that CIR noise (Eq. 3) in the half-saturation constant, which causes variation in the efficiency with which predators exploit their prey and hence in the extent of the population crashes experienced on each passage around the cycle, induces more variability in cycle period than does noise in the prey's population growth rate (Fig. 4, 5, Supporting information; Table 2, cases 2-3). Rescaling our noise sources to equalize the CV of the fluctuations from these sources, rather than their amplitude, seems to mitigate this difference (Table 2, case 3 versus case 5). Yet, when we vary parameters to induce closer passage near the axes and saddles, corresponding to more extreme population losses during cycle troughs, and also to equalize the noise source CVs, the difference between these noise sources re-emerges (Table 2, case 6 versus case 7). A heuristic idea of the cause of this distinction comes from the observation that for fixed predator level $y$, noise in $r$ enters the dynamic equation for $x$ in Eq. 2 in a linear manner, with small deviations in $r$ multiplying the small term $x(1 - x/K)$ that behaves like $x$ itself near $x=0$. On the other hand, noise in $h$ enters the equations like $1/h$, such that even small deviations in $h$ can have an amplifying effect (see the first two sections of the Supporting information for more details).

While we propose that saddle-crawlbys may occur in a range of contexts, it is important to note that the simple presence of a saddle point is not sufficient for our proposed mechanism to play out. A population, cyclic or otherwise, whose trajectory does not pass near the saddle will not experience the dynamical slowing that creates a crawlby. That said, an interesting, and perhaps initially surprising, quantitative observation emerging from our study is that crawlby-induced variability is strongest when the mean cycle trajectory is not too close to the saddle points (Fig. 3A, Supporting information), such that individual cycles represent a mix of crawlbys and non-crawlbys. Indeed, even though crawlbys are responsible for the overall cycle duration variability, the most heterogeneous mixture of long and short passage times past saddles results when crawlbys are only influential on a subset of cycles.

This result also highlights the point that even though several consecutive cycles may occur in which cyclic populations do not come close to extinction, their population dynamics may nonetheless hold the ingredients to produce significant crashes under environmental fluctuations that are within the range of expected occurrence. An interesting direction for future mathematical analysis is to explain the observation from our simulations that as a parameter moves closer to an AH bifurcation point and the deterministic limit cycle correspondingly moves away from the saddle crawlby, the distribution of cycle durations appears to expand in both directions to include both shorter cycles and, counterintuitively, longer cycles (e.g. Table 2, case 2 versus case 8 and case 3 versus case 9).

As with any modeling study, the effects we consider in this work do not provide a complete representation of reality. Real predators may disperse, change behaviors or evolve in ways that change how closely their prey, and thus they, approach extinction (van Baalen and Sabelis 1995). Alternatively, in biological systems, local populations with high-amplitude cycles may routinely go extinct and then undergo recolonization. Our results were conditioned on the populations not going extinct, and we do not model dispersal in or out of the local community nor do we model complex behaviors, alternative prey and other such factors that are likely to influence the dynamics of real populations. Nevertheless, the elements of our model, like population cycles and saddle crawlbys, and their connections to real dynamics are well studied (Murdoch et al. 2003). Our contribution, therefore, is not to replicate or explain the dynamics of any specific single population, but to extend our understanding of cycles and crawlbys as mechanisms through which environmental stochasticity can influence predator-prey dynamics. In doing so, we reconcile a long-recognized mismatch between models and data in their representation of how variable the periods of cyclic populations may be.

In cases in which populations become small during the trough of each cycle, one might reasonably expect demographic stochasticity to be an important factor in cycle duration variability. However, our Gillespie simulations reveal that the effects of environmental stochasticity on cycle period variation are robust to the inclusion of demographic noise. It is nevertheless possible that other dynamic mechanisms not considered in this study may also contribute to temporal variability in cycle durations, and of course even more extreme variability can arise from larger perturbations such as environmental collapse or human interventions. Our intent is not to argue that saddle crawlbys are the only or even the predominant source of cycle period variability, but rather to highlight their putative role.

The observations in this paper represent interesting topics for future mathematical analysis, particularly for models such as Eq. 2 and 3 to which the stochastic averaging methods of Skorokhod et al. (2002) cannot be applied. Previous analytical and numerical work, some done in other biological contexts, has discussed some of the differences in dynamics that arise from the details of how stochasticity appears in a dynamical system (Lande 1993, Goldwyn and Shea-Brown
2011, Allen 2016, Vadillo 2019). In contrast to the individ-
ual-level impact of demographic noise, environmental noise
is shared across entire populations, which is conducive to
larger variability and can increase the likelihood of extreme
events such as extinction (Lande 1993). In theory, we could
have implemented another stochastic process, such as mean-
reverting Ornstein-Uhlenbeck (OU) dynamics, as an alter-
native to Eq. 3. The specific representation of stochasticity
that we have implemented in Eq. 3 features biologically desir-
able properties: it causes the stochastic parameter to remain
positive, to be mean-reverting, and to have some history-
dependence, with a temporal autocorrelation that is nonzero
and decays exponentially (cf. Allen 2016); in our simulations,
it also yields more close passes near the saddle than does
mean-reverting OU dynamics. Our findings reveal that cycle
duration variability can strongly relate to what distribution of
prey densities arises across repeated approaches into a saddle
crawlby, and similar distributions will yield similar variability,
irrespective of which stochastic processes generate them.

To our knowledge, analytical work to date has not addressed how such stochasticity impacts cycle durations and other properties in systems with saddle crawlbys. A related direction that has been pursued analytically (Ashwin and Postlethwaite 2016) is how stochasticity affects the passage of trajectories along a heteroclinic cycle (i.e. a union of orbits that each form a connection from one saddle point to another, which together form a closed loop). These systems may require noise to oscillate, unlike model (2), but nonetheless the resulting oscillations do feature passage near saddle points. Ashwin and Postlethwaite (2016) analyzed residence times near equilibria and switches between neighborhoods of attractors (including heteroclinic cycles), while Giner-Baldó et al. (2017) derived analytical approaches to compute the power spectra of observables in a planar system with a heteroclinic cycle. Both of these analyses were done on systems with stochasticity in the form of additive Gaussian white noise, however, whereas the CIR noise that we consider poses a new challenge for future work.

Our results update our understanding of how environmen-
tal stochasticity can contribute meaningfully to observed cycle
period variation in nature. While our results for stochastic
prey growth rate, r, largely confirm past results that stochastic-
ity has little effect on cycle period, we have shown that some
noise sources – specifically those that induce significant vari-
ability in the distance from saddle points as part of the normal
cycle – can lead to substantial cycle period variation. We find
that this effect is strongest when the mean cycle trajectory is
not too close to the saddle points, such that individual cycles
represent a mix of crawlbys and non-crawlbys. In these cases,
the coefficient of variation in cycle period approaches 0.3,
the mean value estimated from time series data. We therefore
conclude that variance in the duration of saddle crawlbys pro-
vides a viable explanation for some, though not all, observed
instances of variability in cycle periods.

*Funding* – This work was initiated at the ‘Transients in Biological Systems’ Investigative Workshop at the National Institute for Mathematical and Biological Synthesis (NIMBioS), supported

by the National Science Foundation through NSF Award DBI-
1300426, with additional support from The University of
Tennessee, Knoxville. Additional progress was made via a meeting
at the University of Pittsburgh supported by the Dietrich School of
Arts and Sciences through the Mathematics Research Center. JER
acknowledges support from NSF Award DMS-1951095. DJDE
was supported by the Natural Sciences and Engineering Research
Council of Canada (NSERC).

Author contributions

**Jonathan Rubin:** Conceptualization (equal); Formal analysis (supporting); Investigation (lead); Methodology (lead); Project administration (lead); Software (lead); Visualization (lead); Writing – original draft (equal); Writing – review and editing (lead). **David Earn:** Conceptualization (equal); Investigation (supporting); Methodology (supporting); Software (supporting); Writing – review and editing (supporting). **Priscilla Greenwood:** Conceptualization (equal); Formal analysis (supporting); Investigation (supporting); Writing – review and editing (supporting). **Todd Parsons:** Formal analysis (lead); Investigation (supporting); Methodology (supporting); Writing – review and editing (supporting). **Karen Abbott:** Conceptualization (equal); Data curation (lead); Investigation (supporting); Methodology (supporting); Visualization (supporting); Writing – original draft (equal); Writing – review and editing (supporting).

Data availability statement

Codes are available at the GitHub repository jonathanerubin/irregularPopulationCycles (https://github.com/jonathanerubin/irregularPopulationCycles).

**Supporting information**

The Supporting information associated with this article is
available with the online version.

References

Abbott, K. C. 2020. Non-equilibrium dynamics and stochastic pro-
cesses. – In: McCann, K. S. and Gellner, G. (eds), Theoretical
ecology: concepts and applications. Oxford Univ. Press, pp. 40–52.

Allen, E. 2016. Environmental variability and mean-reverting pro-
cesses. – Discrete Contin. Dyn. Syst. B 21: 2073–2089.

Allen, L. J. S. and Allen, E. J. 2003. A comparison of three differ-
ent stochastic population models with regard to persistence
time. – Theor. Popul. Biol. 64: 439–449.

Anderson, R. M. and May, R. M. 1980. Infectious diseases and
population cycles of forest insects. – Science 210: 658–661.

Ashwin, P. and Postlethwaite, C. 2016. Quantifying noisy attrac-
tors: from heteroclinic to excitable networks. – SIAM J. Appl.
Dyn. Syst. 15: 1989–2016.

Aubry, K. B. and Lewis, J. C. 2003. Extirpation and reintroduction
of fishers *Martes pennanti* in Oregon: implications for their
conservation in the Pacific states. – Biol. Conserv. 114: 79–90.

Barraquand, F., Louca, S., Abbott, K. C., Cobbold, C. A., Cor-
doleani, F., DeAngelis, D. L., Elderd, B. D., Fox, J. W., Green-
wood, P., Hilker, F. M., Murray, D. L., Stieha, C. R., Taylor, R.
A., Vitense, K., Wolkowicz, G. S. K. and Tyson, R. C. 2017. Moving forward in circles: challenges and opportunities in modelling population cycles. – *Ecol. Lett.* 20: 1074–1092.

Benedetti-Cecchi, L. 2003. The importance of the variance around the mean effect size of ecological processes. – *Ecology* 84: 2335–2346.

Bengfort, M., van Velzen, E. and Gaedke, U. 2017. Slight phenotypic variation in predators and prey causes complex predator–prey oscillations. – *Ecol. Complex.* 31: 115–124.

Benincà, E., Ballantine, B., Ellner, S. P. and Huisman, J. 2015. Species fluctuations sustained by a cyclic succession at the edge of chaos. – *Proc. Natl Acad. Sci. USA* 112: 6389–6394.

Boettiger, C. and Hastings, A. 2013. No early warning signals for stochastic transitions: insights from large deviation theory. – *Proc. R. Soc. B* 280: 20131372.

Cazelles, B., Chavez, M., Berteaus, D., Ménard, F., Vik, J. O., Jenouvrier, S. and Stenseth, N. C. 2008. Wavelet analysis of ecological time series. – *Oecologia* 156: 257–304.

Coulson, T., Rohani, P. and Pascual, M. 2004. Skeletons, noise and population growth: the end of an old debate? – *Trends Ecol. Evol.* 19: 359–364.

Cox, J. C., Ingersoll, J. E. and Ross, S. A. 1985. A theory of the term structure of interest rates. – *Econometrica* 53: 385–407.

Cushing, J. M., Dennis, B., Desharnais, R. A. and Costantino, R. F. 1998. Moving toward an unstable equilibrium: saddle nodes in population systems. – *J. Anim. Ecol.* 67: 298–306.

Dwyer, G., Dushoff, J. and Yee, S. H. 2004. The combined effects of pathogens and predators on insect outbreaks. – *Nature* 430: 341–345.

Ermentrout, B. 2002. Simulating, analyzing and animating dynamical systems: a guide to XPPAUT for researchers and students. – *SIAM*.

Francis, T. B., Abbott, K. C., Cuddington, K., Gellner, G., Hastings, A. M., Lai, Y.-C., Morozov, A., Petrovskii, S. and Zeeman, M. L. 2021. Management implications of long transients in ecological systems. – *Nat. Ecol. Evol.* 5: 285–294.

Giner-Baldó, J., Thomas, P. J. and Lindner, B. 2017. Power spectrum of a noisy system close to a heteroclinic orbit. – *J. Stat. Phys.* 168: 447–469.

Goldwyn, J. H. and Shea-Brown, E. 2011. The what and where of adding channel noise to the Hodgkin–Huxley equations. – *PLoS Comput. Biol.* 7: e1002247.

Greenman, J. and Pasour, V. 2011. Phase control of resonant systems: interference, chaos and high periodicity. – *J. Theor. Biol.* 278: 74–86.

Hastings, A., Abbott, K. C., Cuddington, K. M., Francis, T., Gellner, G., Lai, Y.-C., Morozov, A., Petrovskii, S., Scranton, K. and Zeeman, M. L. 2018. Transient phenomena in ecology. – *Science* 361: eaat6412.

Hastings, A., Abbott, K. C., Cuddington, K. M., Francis, T., Lai, Y.-C., Morozov, A., Petrovskii, S. and Zeeman, M. L. 2021. Effects of stochasticity on the length and behavior of ecological transients. – *J. R. Soc. Interface* 18: 20210257.

Holling, C. S. 1965. The functional response of predators to prey density and its role in mimicry and population regulation. – *Memoirs Entomol. Soc. Can.* 97: 5–60.

Holyoak, M. and Wetzel, W. C. 2020. Variance-explicit ecology: a call for holistic study of the consequences of variability at multiple scales. – In: Dobson, A., Tilman, D. and Holt, R. D. (eds), *Unsolved problems in ecology*. Princeton Univ. Press, pp. 25–42.

Ives, A. R., Einarsson, Á., Jansen, V. A. A. and Gardarsson, A. 2008. High-amplitude fluctuations and alternative dynamical states of midges in Lake Myvatn. – *Nature* 452: 84–87.

Jäger, C. G., Diehl, S., Matauschek, C., Klausmeier, C. A. and Stibor, H. 2008. Transient dynamics of pelagic producer–grazer systems in a gradient of nutrients and mixing depths. – *Ecology* 89: 1272–1286.

Kendall, B. E., Briggs, C. J., Murdoch, W. W., Turchin, P. V., Ellner, S. P., McCauley, E., Nisbet, R. M. and Wood, S. N. 1999. Why do populations cycle? A synthesis of statistical and mechanistic modeling approaches. – *Ecology* 80: 1789–1805.

Lande, R. 1993. Risks of population extinction from demographic and environmental stochasticity and random catastrophes. – *Am. Nat.* 142: 911–927.

Morozov, A., Abbott, K. C., Cuddington, K. M., Francis, T., Gellner, G., Hastings, A., Lai, Y.-C., Petrovskii, S., Scranton, K. and Lou Zeeman, M. 2020. Long transients in ecology: theory and applications. – *Phys. Life Rev.* 32: 1–40.

Murdoch, W. W., Briggs, C. J. and Nisbet, R. M. 2003. Consumer-resource dynamics. – *Princeton Univ. Press*.

Myers, J. H. 2018. Population cycles: generalities, exceptions and remaining mysteries. – *Proc. R. Soc. B* 285: 20172841.

Nguyen, H. T. H. and Rohani, P. 2007. Noise, nonlinearity and seasonality: the epidemics of whooping cough revisited. – *J. R. Soc. Interface* 5: 403–413.

Nisbet, R. M. and Gurney, W. C. S. 2003. Modelling fluctuating populations. – *Blackburn Press*.

Prendergast, J., Bazeley-White, E., Smith, O., Lawton, J., Inchausti, P., Kidd, D. and Knight, S. 2010. The global population dynamics database. – *urn:node:KNB.doi:10.5063/F1BZ63Z8*.

Rosenzweig, M. L. and MacArthur, R. H. 1963. Graphical representation and stability conditions of predator–prey interactions. – *Am. Nat.* 97: 209–223.

Shoemaker, L. G., Sullivan, L. L., Donohue, I., Cabral, J. S., Williams, R. J., Mayfield, M. M., Chase, J. M., Chu, C., Harpole, W. S., Huth, A., HilleRisLambers, J., James, A. R. M., Kraft, N. J. B., May, F., Muthukrishnan, R., Satterlee, S., Taubert, F., Wang, X., Wiegand, T., Yang, Q. and Abbott, K. C. 2020. Integrating the underlying structure of stochasticity into community ecology. – *Ecology* 101: e02922.

Sibly, R. M., Barker, D., Hone, J. and Pagel, M. 2007. On the stability of populations of mammals, birds, fish and insects. – *Ecol. Lett.* 10: 970–976.

Skorokhod, A. V., Hoppensteadt, F. C. and Salehi, H. 2002. Random perturbation methods with applications in science and engineering. – Springer.

Turchin, P. V. and Hanski, I. 2001. Contrasting alternative hypotheses about rodent cycles by translating them into parameterized models. – *Ecol. Lett.* 4: 267–276.

Turchin, P., Wood, S., Ellner, S. P., Kendall, B. E., Murdoch, W. W., Fischlin, A., Casas, J., McCauley, E. and Briggs, C. J. 2003. Dynamical effects of plant quality and parasitism on population cycles of larch budmoth. – *Ecology* 84: 1207–1214.

Vadillo, F. 2019. Applied mathematics and computation. – *Appl. Math. Comput.* 360: 181–189.

van Baalen, M. and Sabelis, M. W. 1995. The milker-killer dilemma in spatially structured predator–prey interactions. – *Oikos* 74: 391–400.

Violle, C., Enquist, B. J., McGill, B. J., Jiang, L., Albert, C. H., Hulshof, C., Jung, V. and Messier, J. 2012. The return of the variance: intraspecific variability in community ecology. – *Trends Ecol. Evol.* 27: 244–252.