Irregular population cycles driven by environmental stochasticity and saddle crawlbys

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.


Introduction
Recurrent population cycles -eruptive increases and crashes in population sizeare widespread (Kendall et al. 1999, Barraquand et al. 2017. Consumer-resource and host-enemy interactions are famously known to promote cycles 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 Irregular population cycles driven by environmental stochasticity and saddle crawlbys 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 susceptibleinfected 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 andPasour 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 . 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 a 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 tightlycoupled 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 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 x e rt 0, 0, = out i n 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  , 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 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 µ and σ 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~ N(ln(c)/r − µ/r, (σ/r) 2 ) (we use N(µ, σ 2 ) to denote the normal distribution with mean µ and variance σ 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 )~ N((1 − m/r)ln(c) + mµ/r, (mσ/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 lognormally 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.

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 . 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). Rosenzweig-MacArthur predator-prey model 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 r or 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 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 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. and hence, as discussed in the 'Saddle crawlbys 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 p p dt p dW 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 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 crawlbys, because changes in h exaggerate the influence of crawlbys 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 r of 1.0 might reasonably have less effect than the same σ = 0.2 acting on an h of 0.15. We tested this idea in two ways. First, we reduced 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 on r (by rescaling by 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 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. 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 minimum size cutoff of 10 −5 .

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 variability. 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 environmental 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 prescribed expected rates, and events are implemented sequentially 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 simulations yielded very similar distributions of cycle durations to our direct simulations of system (2)-(3), matching the differences 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 methodological choices.

Cycle variability is driven by variable passage times near extinction and prey carrying capacity
We have shown how parameter changes can affect the proximity of cycles to the saddle points (Fig. 3) and how stochastic 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 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 by prolonging the time during which one or both populations 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 complete a crawlby of the extinction state is linearly related to ln(x 0,in ), the logarithm of the prey density when the populations 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 0,in ) and the passage time by the (0, 0) saddle (Fig. 7A, B). There is also a strong correlation between log prey density, ln(x 0,in ), 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 simulation results with those from simulating Eq. 2 and 3 provides additional evidence that successive cycles are effectively independent. 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 variability, while if γ were too small, then significant historydependence 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 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). entry into the two saddle neighborhoods differ, depending 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 logarithmic 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 σ = 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 consumerresource 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, 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 . 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 individual-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 meanreverting Ornstein-Uhlenbeck (OU) dynamics, as an alternative to Eq. 3. The specific representation of stochasticity that we have implemented in Eq. 3 features biologically desirable properties: it causes the stochastic parameter to remain positive, to be mean-reverting, and to have some historydependence, 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 environmental 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 stochasticity has little effect on cycle period, we have shown that some noise sources -specifically those that induce significant variability 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 provides 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).

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

S1 | EFFECTS OF STOCHASTICITY ON THE APPROACH TO THE SADDLE CRAWLBY
The Rosenzweig-MacArthur equations, which we rewrite here: are unfortunately analytically intractable, but the stable limit cycle has segments that are approximately parallel to the x and y axes, which approach closer and closer to the axes as r or h decreases. This allows us to approximate the trajectory via the method of boundary layers (O'Malley, 1991;Kevorkian and Cole, 1996), linearizing along the axes.
To do so, we introduce a small, dimensionless parameter ϵ, and substitute x (t ) = ϵξ (t ) or y (t ) = ϵψ (t ), as appropriate, into system (S1). We then take the limit as ϵ → 0. This is equivalent to focusing on the part of the orbit where the transformed variable(s) are small, and allows us to obtain approximate solutions along the x -axis and y -axis and at the origin. Similarly, letting x (t ) = K + ϵξ (t ) and y (t ) = ϵψ (t ), we can approximate the trajectory near the We illustrate the method by computing the solution along the y -axis: substituting x (t ) = ϵξ (t ) into (S1) gives In the limit as ϵ → 0, these equations simplify to which may be solved to yield y (t ) = y (0)e −mt and ξ (t ) = ξ (0)e r t − a y (0) mh (1−e −mt ) , where y (0) and x (0) = ϵξ (0) are the (unspecified) point of entry into a neighborhood of the y -axis. We may also return to the original variables, taking When h is allowed to vary stochastically according to the Cox-Ingersoll-Ross (CIR) process, which we write here for an arbitrary parameter p and Wiener process W (t ) as we need to modify our solution to (S3) to account for the fact that h is a (random) function of time, giving us Proceeding similarly, we can obtain approximate local solutions for stochastically varying r or h, which we summarize in Table S2.
Comparing the solutions along the y -axis for varying r vs. varying h gives some intuition for the greater impact of noise in the latter: in this region, the noise enters as 1/h, whereas h itself fluctuates about a mean ofh = 0.15, and thus will regularly approach values very close to 0. If this happens near the y -axis, x (t ), which is proportional to can become quite small. This can cause the perturbed trajectory to make a very close approach to the unstable node at the origin and the consequent slow-down of the passage. By contrast, in the same region, noise in r enters linearly, and only a rare, large fluctuation in r will cause the trajectory to have a similar excursion towards the crawlby region near the origin.

S2 | AVERAGING ARGUMENT FOR WHY NOISE THAT ENTERS SYSTEM (S1) NONLINEARLY YIELDS DIFFERENT EFFECTS THAN NOISE THAT ENTERS LINEARLY
Here, we make some informal asymptotic analyses providing further evidence for qualitative differences in the consequences of noise occurring linearly vs. nonlinearly in the equations governing the dynamics. Let P (t ) be any parameter (which could be one or more of a, b, r , m or h) that we treat as stochastic, with time evolution governed by the CIR equation (S4), which we restate here: Let ϵ be a small dimensionless parameter and set P ϵ (t ) = P t ϵ . P ϵ (t ) satisfies a rescaled stochastic differential equation where, as previously, W (t ) is a Wiener process. P (t ) and P ϵ (t ) have the same stationary distribution, which is gamma distributed with parameters α = 2γp σ 2 and β = 2γ σ 2 and probability density function The rescaling accelerates the convergence to the stationary distribution.
Informally, in the limit as ϵ → 0, the process P ϵ (t ) becomes a family (indexed by t ∈ ) of independent draws from the stationary distribution for P (t ); in particular, given an arbitrary differentiable function f , for sufficiently small ∆t , appealing to Birkhoff's ergodic theorem (Birkhoff, 1931), we have as ϵ → 0, where, as previously, π (p ) is the probability density function of the stationary distribution for P (t ).
Intuitively, based on (S7), we might expect that the solution to an arbitrary system of ordinary differential equa-3 tions driven by P ϵ (t ), say should be well approximated in the limit as ϵ → 0 by the averaged system Indeed, this line of reasoning can be made rigorous provided that the equation (S8) has a unique solution (Kurtz, 1992).
If we apply this approximation to our predator-prey system (S1), we see that after replacing any of a, b, r , m by a rapidly mixing CIR noise, say A ϵ , B ϵ , etc. with meansā,b, etc., respectively, in the limit as ϵ → 0 the noise term is simply replaced by its mean value (so a is replaced byā and so forth) and we recover the original system of ordinary differential equations.
By contrast, if we replace h by a rapidly mixing CIR noise H ϵ (t ) with meanh, then we can use (S6) to write the average of the terms involving h as where α = 2γh σ 2 and β = 2γ σ 2 . The latter integral equals (see e.g., NIS, 2020, equation 8.6.4) which may be more compactly expressed using confluent hypergeometric equations (see e.g., Abramowitz and Stegun, 1965) as βU (1, 2 − α, β x ). While the latter approaches an asymptote of 1/(x +h ) for large values of x , the two do not coincide; for example, Thus in the specific case of noise in h, averaging applied to (S1) yields the novel system of averaged equations which provides another argument for why noise in a term appearing nonlinearly in (S1), such as h or K , will impact the dynamics differently than noise that enters linearly.

S3 | SIMULATIONS WITH DEMOGRAPHIC NOISE
For our simulations with demographic noise, we used the standard Gillespie algorithm (Gillespie, 1976(Gillespie, , 1977 to recast the system (S1), (S4) as a continuous-time Markov jump process, whereby exponential holding times alternate with Markov events (each term in the differential equations (S1) is treated as an event rate for changes in the numbers of predators and prey). The degree of demographic stochasticity in the Gillespie simulations is determined by the 4 prey carrying capacity K , which we set to 8000 in all our simulations. Note that in this framework, K is a population scale. Thus, no other parameters need to be adjusted when we tune K , but the initial conditions for the simulations do need to be scaled with K . The prey and predator densities (i.e., the numbers re-scaled by the carrying capacity K ) obtained from the Gillespie algorithm are very well approximated by the equations (S1) when K is large (Kurtz, 1971), and preliminary simulations (with initial conditions represented as proportions of K ) showed K = 8000 to be sufficiently large without being too computationally expensive. To avoid extinctions, we included an additional, small, constant immigration rate given by the product of parameters a and b; at each immigration event, one predator and ⌈1/b ⌉ prey were added. The time to the next event is drawn from an exponential distribution with mean equal to the inverse of the total rate at which events of any type occur. The probability that the event is of a given type is the ratio of its event rate to the total event rate. These simulations emulate demographic noise effects. To increase simulation efficiency, we used the adaptive tau-leaping approximation to the exact Gillespie algorithm (Gillespie, 2007).
The standard Gillespie algorithm implicitly assumes that the event rates are constant, whereas we needed to be able to modulate either the prey growth rate (r ) or the predator's half-saturation constant (h) stochastically. Since we have found that autocorrelated noise in the input (r or h) affects the variance in the output (cycle length), we used an Euler-Maruyama numerical scheme to evolve the focal parameter between the jump times in the Markov chain.
More precisely, if t n is the time of the nth event (or nth τ-leap), ∆t n = t n − t n −1 , and N n is a draw from a standard Normal(0,1) distribution, then the value of focal parameter p at time t n (p n = p (t n )) is obtained via the following algorithm (cf. Table 1 in the main text for parameter meanings): n ← 0 t n ← 0 • Draw initial value p 0 from the stationary (CIR) distribution while (t n < t max ) do n ← n + 1 • Generate time ∆t n to next jump using rate p n −1 in the τ-leap algorithm • Draw N n from a Normal(0, 1) distribution p n ← p n −1 + γ (p − p n −1 )∆t n + σ p n ∆t n N n t n ← t n −1 + ∆t n end while Note that √ ∆t n N n is equivalent to a draw from the Wiener process W n ∼ Normal(0, ∆t n ). The time steps ∆t n in the Euler-Maruyama algorithm are assumed to be small, which we have verified in our simulations (e.g., 99% of time steps ∆t n are smaller than 0.01 and the median ∆t n is less than 0.001).

S4 | ENTRY DENSITY DISTRIBUTIONS
The densities of predators and prey at entry to the neighborhoods of saddle points, as shown in Figure 8 in the main text, can be empirically fit by a lognormal distribution. This result may seem counter to reasonable naïve expectations, which might suggest that these densities should reflect the stationary distribution for the CIR process, which is a gamma distribution. This difference may result because population dynamics act as a complex filter, transforming the input noise in unexpected ways. We conjecture that the lognormal fit arises as a consequence of the central limit theorem, which states roughly that "the sum of many independent random variables will be approximately normally 5 distributed if each summand has high probability of being small" (Billingsley, 1995). To briefly explain this logic, we write system (S1) abstractly as (S10) Using 0, T 1 , T 2 , . . . to denote successive crossing times for a specific entry boundary, we can formally express the x -component of the solution of system (S10) on the time interval [0, T k ] for any natural number k as Taking a logarithm of both sides yields If the noise is not too strong, then the time elapsed between subsequent crossings, T k − T k −1 , is approximately equal to a period and each integral f (x (s ), y (s ) ) ds is small; indeed, with no noise, the integral over the deterministic limit cycle through one period vanishes. Thanks to these properties, along with the exponentially fast decay in noise correlation and the exponential stability of the limit cycle, which give rise to the rapid decay of correlation between subsequent boundary crossings, the sum (S11) approximately satisfies the conditions for a central limit theorem that would imply that {ln(x (T k ) ) } are normally distributed. We leave a rigorous treatment for future work.

S5 | GENERALIZATION TO OTHER MODELS
Although the hypothetical possibilities for alternative models to system (S1) are vast, past authors (Turchin and Batzli, 2001;Krebs, 2013) have identified 6-8 major classes of models that encompass the predominant types of vegetationherbivore population interactions, four of which also apply to predator-prey dynamics. The first model in this classification is exactly system (S1), and another includes time-dependent seasonality, which is outside the scope of this work. We therefore focus on the remaining two classes of predator-prey model.
The Bazykin model features predator (or herbivore) self-limitation. The model equations take the form with d > 0. In fact, algebraic calculations show that under the same assumption, ab > m, that we have already imposed, this system has a unique coexistence equilibrium in the positive (x, y ) quadrant, which can undergo a supercritical AH bifurcation just as occurs in model (S1). The predator self-limitation term −d y 2 does not affect the existence of the (0, 0) and (K , 0) predator-free equilbria and, since it is nonlinear, also does not impact calculations based on linearization about these equilibria, which remain saddles. In short, the dynamics of system (S12) agree qualitatively with those of (S1) in the parameter regime of interest, and we observe similar effects when we impose environmental stochasticity of the form (S4) in individual parameters of model (S12) as well. 6 The remaining model class yields results that are consistent with these and also provides some new observations. This variable territory model also includes predator (or herbivore) self-limitation but the intensity of this effect depends on prey (or vegetation) density, fading away when food supplies are abundant. The model equations are again with d > 0. For this model, (K , 0) is again a saddle point. Dynamics with x = 0 is not well-defined, but this can be addressed by replacing d y 2 /x with d y 2 /(x + ε ) for a small value ε > 0, in which case the origin also becomes a saddle and the dynamics is otherwise similar. Henceforth we will ignore dynamics with x ≤ 0, since the positive quadrant is the biologically relevant part of the phase space and is invariant.
Model (S13) has stable oscillations in some parameter regimes, including our baseline parameter set, with the one change of setting m = 0.5, with c = 0.1. We introduced noise separately in r and h in model (S13) using equation (S4) with σ = 0.1 and 0.2 and performed numerical experiments analogous to those done for model (S1).
We found that, as with model (S1), the mean cycle duration was similar for all cases but the corresponding CV was much greater with noise in h and increased with σ (noise in r , σ = 0.1: mean 16.09, CV 0.032; noise in h, σ = 0.1: mean 16.24, CV 0.066; noise in r , σ = 0.2: mean 16.05, CV 0.048; noise in h, σ = 0.2: mean 17.23, CV 0.13). In addition, some novel insights emerged from our numerical simulations with system (S13). First, this model showed higher variability of passage times through region F in Figure 2B in the main text than seen with model (S1). This effect appears to emerge from the sensitivity of the term d y 2 /x (or d y 2 /(x + ε )), and hence of the rate of decay of y in region F, to small changes in x . Thus, while saddle crawlbys remain as important sources of variability here, this system illustrates how other effects can also contribute to cycle duration variability. Second, the variability of passage times through region A depends significantly on the position of the entry boundary, {x = x 0,in }, to the region. With noise in h, the CV of the region A passage time distributions from our simulations with σ = 0.1 and σ = 0.2 approximately doubled if we changed x 0,in from 0.1 to 0.08. This observation highlights the point that saddle crawlbys are local, and the neighborhood of the saddle on which they exert a dominant influence on system dynamics may have different sizes in different systems. Because the term d y 2 /x in the variable territory model (S13) with predator self-limitation dependent on prey density exerts its own strong influence when x is small, the impact of (0, 0) becomes more limited, while the saddle point at (K , 0), away from extinction, still strongly impacts overall cycle duration variability. 8 TA B L E S 1 Conditions used to define the phase space regions shown in Figure 2B in the main text. These values were used for all results except cases 6-9 in Table 2 in the main text, for which sections were slightly adjusted to reflect altered orbit paths. 9 TA B L E S 2 Solutions obtained by linearizing the model defined by equations (S1), (S4) and in different regions of phase space.  Figure 1A-B, we arbitrarily considered a population to be at "low" abundance when it fell below 5% of its maximum observed abundance. Here, we show that the patterns are insensitive to the exact cutoff used (A: 1%, B: 5% as in the main text, and C: 10%). As in Figure 1A-B, the blue histogram shows time series with only short (0-3 time steps) forays to low abundance, and red shows histograms with longer visits to low abundance; vertical blue and red lines mark the corresponding means.  Figure 2: A: near the origin, B: along the x -axis, C: near the (K , 0) saddle point, D: along the unstable manifold of (K , 0), E: traversing the y -nullcline, F: along the y -axis. For noise in h, note that the distributions in regions B and F have sharp minima corresponding to the passage times generated by dx /d t = r x (1 − x /K ) and d y /d t = −m y , respectively, which are approached on cycles for which the stochastic ax y /(x + h ) term from system (S1) (the functional response) becomes so small that it has very limited influence over the passage times along the axes. The numbers of cycles incorporated in these histograms are as in Figure 5.
F I G U R E S 4 Mean passage times through regions A-F (see Figure 2) for noise in r (top) and noise in h (bottom). In each plot, blue bars show data obtained from simulation of the stochastic differential equation system (S1)-(S4) while red bars are from Gillespie simulations. Error bars represent standard deviations. From left to right in each group, noise amplitude σ increased from 0.1 to 0.2 to 0.3. Note the different y-axis scales in the different plots. F I G U R E S 5 Independence from one cycle to the next of the population densities upon entry into the neighborhoods of the saddles. The stochastic parameter is prey growth rate, r , in A,C and predator's half-saturation constant, h, in B,D. A-B show the lack of a correlation between prey density upon entry into the (0, 0) neighborhood, ln(x 0,in ), from the n th to the n + 1 st cycle. C-D show the same for predator density upon entry into the (K , 0) neighborhood, ln(y K ,in ). In all cases, σ = 0.2 and axes are on logarithmic scales. Note that results in B include some values of x n+1 below the extinction cutoff of 10 −5 ; these were not used in any of the duration calculations discussed in the main text.
14 A B C D E F F I G U R E S 6 Distributions of prey densities, ln(x 0,in ), at crossing of a Poincaré section of constant y = y 0,in , which marks entry into the neighborhood of (0, 0); see Figure 2A in the main text. Results are shown for trajectories obtained by simulation of (S1), (S4); the noisy parameter is r in the left panels (with y 0,in = 0.07) and h in the right panels (with y 0,in = 0.12). A,B: σ = 0.1, C,D: σ = 0.2, E,F: σ = 0.3. In each panel, the red curve shows the best lognormal fit to each distribution.

A B
C D E F F I G U R E S 7 Distributions of predator densities, y K ,in , at crossing of a Poincaré section of constant x = x K ,in , marking entry into the neighborhood of (K , 0). Results are shown for trajectories obtained by simulation of (S1), (S4); the noisy parameter is r in the left panels (with x K ,in = 0.6) and h in the right panels (with x K ,in = 0.5). A,B: σ = 0.1, C,D: σ = 0.2, E,F: σ = 0.3. In each panel, the red curve shows the best lognormal fit to each distribution.