# No sensitivity to functional forms in the Rosenzweig-MacArthur model with
strong environmental stochasticity

Frédéric Barraquand  Institute of Mathematics of Bordeaux, CNRS & University
of Bordeaux, Talence, France

###### Abstract

The classic Rosenzweig-MacArthur predator-prey model has been shown to
exhibit, like other coupled nonlinear ordinary differential equations (ODEs)
from ecology, worrying sensitivity to model structure. This sensitivity
manifests as markedly different community dynamics arising from saturating
functional responses with nearly identical shapes but different mathematical
expressions. Using a stochastic differential equation (SDE) version of the
Rosenzweig-MacArthur model with the three functional responses considered by
Fussmann & Blasius (2005), I show that such sensitivity seems to be solely a
property of ODEs or stochastic systems with weak noise. SDEs with strong
environmental noise have by contrast very similar fluctuations patterns,
irrespective of the mathematical formula used. Although eigenvalues of
linearised predator-prey models have been used as an argument for structural
sensitivity, they can also be an argument against structural sensitivity.
While the sign of the eigenvalues’ real part is sensitive to model structure,
its magnitude and the presence of imaginary parts are not, which suggests
noise-driven oscillations for a broad range of carrying capacities. I then
discuss multiple ways to evaluate structural sensitivity in a stochastic
setting, for predator-prey or other ecological systems.

Keywords: predator-prey; functional response; structural sensitivity;
stochastic differential equation; multiplicative noise

Correspondence to<EMAIL_ADDRESS>
## Introduction

Structural sensitivity, as defined by Wood & Thomas (1999) and Fussmann &
Blasius (2005), occurs when small changes in model structure can induce large
changes in model behaviour. As Fussmann & Blasius (2005) demonstrated for the
continuous-time predator-prey model of Rosenzweig & MacArthur (1963), picking
a particular equation for the functional response (i.e., the kill rate of
individual predators as a function of prey density), out of a set of functions
that are barely different to the eye, can result in the long run in either a
limit cycle or a stable fixed point. Surprisingly, the phenomenon can occur
for sizeable sections of parameter space, which makes us reflect upon the
arbitrary choices made by modellers: the commonly used Holling (also called
Michaelis-Menten) form, for instance, always tends to generate cycles more
easily. Structural sensitivity has been judged to be one of the main limits to
the predictive power of dynamic models (Fussmann & Blasius, 2005; Adamson &
Morozov, 2013), and it has even been suggested that such sensitivity may
reflect the inherent complexity of ecological processes (Adamson & Morozov,
2014a; Adamson _et al._ , 2016). Several mathematical approaches have emerged
to both detect such structural sensitivity of nonlinear ordinary differential
equations (Adamson & Morozov, 2013, 2014a, 2014b), and mitigate its
consequences, by averaging the predictions of multiple plausible models
(Adamson _et al._ , 2016; Aldebert _et al._ , 2018).

The perspective adopted here consists instead of viewing structural
sensitivity from a stochastic standpoint. Nisbet & Gurney (1976, 1982)
stressed early on that an underdamped oscillator sufficiently perturbed by
noise could exhibit quasi-cycles, i.e., noise-sustained oscillations with a
dominant periodicity. There is therefore reason to believe that an oscillator
that presents a stable focus with a pair of complex eigenvalues (instead of a
stable node, which only presents real eigenvalues) will often show cycles in a
stochastic environment. Much of the research on quasi-cycles has focused on
the influence of demographic stochasticity (McKane & Newman, 2004, 2005;
Pineda-Krch _et al._ , 2007), which means relatively small stochastic
perturbations in large populations. With demographic perturbations, effects of
noise will certainly be visible, but they may be somewhat limited, and
consequently not strong enough to counteract the sensitivity of deterministic
models to functional forms—unless populations are very small. This creates an
incentive to investigate what happens in the presence of large environmental
perturbations, that are also believed to dominate the stochasticity
experienced by most large natural populations (Lande _et al._ , 2003;
Mutshinda _et al._ , 2009). An overlooked phenomena in random dynamical
systems is that of “noisy precursors” (Wiesenfeld, 1985): close to and below a
supercritical Hopf bifurcation in the deterministic skeleton of the model, the
behaviour of its stochastic counterpart will make it look like the bifurcation
has already been crossed (Wiesenfeld, 1985; Neiman _et al._ , 1997; Wyse _et
al._ , 2022). Or, in other words, environmental noise smoothes over the
transition between fixed point stability and cyclicity. This also suggests
that strong stochasticity may remove some of the near-pathological behaviour
of deterministic models.

In the following, I introduce strong environmental stochasticity in the
population growth rates of prey and predator in a stochastic differential
equation (SDE) version of the Rosenzweig-MacArthur model (Rosenzweig &
MacArthur, 1963; Rosenzweig, 1971). I show numerically that the amplitude of
cycles reacts to the bifurcation parameter $K$ in smooth and similar manner
for the three functional responses considered by Fussmann & Blasius (2005), so
that the sensitivity to structure is essentially removed from the (strongly)
stochastic models. A reinterpretation of how real and imaginary eigenvalues of
the deterministic skeleton of the model change with $K$ further illustrates
why such stochasticity smoothes dynamical transitions and largely forbids
structural sensitivity.

## Model and methods

The stochastic Rosenzweig-MacArthur (RMA) model is formulated as a system of
coupled stochastic differential equations (SDEs) in Itô form, defined by eqs.
1–2 for the prey ($X$) and the predator ($Y$) densities, that are random
variables evolving over time:

$\displaystyle dX={}$
$\displaystyle\left(g(X)-f(X)Y\right)dt+\sigma_{E}XdW_{x}$ (1) $\displaystyle
dY={}$ $\displaystyle\left(\epsilon f(X)Y-mY\right)dt+\sigma_{E}YdW_{y}.$ (2)

We use in this article the exact same deterministic skeleton as Fussmann &
Blasius (2005), with three functional forms for the functional response,
$f_{H}(x)=\frac{a_{H}x}{1+b_{H}x}$,$f_{I}(x)=a_{I}\left(1-\exp(-b_{I}x)\right)$,
$f_{T}(x)=a_{T}\tanh(b_{T}x)$ and a logistic prey growth rate function
$g(x)=rx\left(1-\frac{x}{K}\right)$. The functional response parameters are
identical to Fussmann & Blasius (2005), $a_{H}=3.05,b_{H}=2.68$;
$a_{I}=1,b_{I}=2$; $a_{T}=0.99,b_{T}=1.48$. We also use the same other
parameters values $r=1,\epsilon=1,m=0.1$ (the conversion efficiency was added
to be able to also consider more realistic $\epsilon<1$ values, though we
stick to $\epsilon=1$ for main text simulations). $W_{x}(t)$ and $W_{y}(t)$
are uncorrelated Wiener processes. The environmental noise standard deviation
is set to $\sigma_{E}=0.25$ ($\sigma_{E}=0.1$ is considered in Appendix A).

The numerical integration is performed using the simplest possible Euler-
Maruyama scheme, with steps of size $\delta t=0.05$ for $K\in[0,1.5]$, with
integration over $2000$ timesteps, and $\delta t=0.005$ for $K\in[0,15]$
(Appendix A). The ODE integration uses an Euler scheme for the same step size;
Matlab’s `ode45` provides a good match, suggesting that it is a sufficiently
small. The eigenvalues of the nontrivial fixed point of the deterministic
model have been computed through an analytical expression of the Jacobian,
computed numerically at the equilibrium for the chosen parameter values. The
computer code is available at
https://github.com/fbarraquand/sensitivity_stochRMA. The original code, used
for the figures presented here, has been written and executed in Matlab
(R2017a). As always with numerical results, there is a risk of a different
kind of sensitivity to the computational implementation of the models.
Therefore, the main result of Fig 2 has been replicated in C++ for some
parameter values (see code repository).

## Results

Comparisons of the deterministic and stochastic trajectories (example for
$K=0.5$ in Fig. 1) reveal that while the dynamical behaviour of the ODEs for
Holling, Ivlev and tanh functional responses can be markedly different, the
stochastic trajectories are not too dissimilar. For this value of $K$, only
the Holling functional response model would have crossed the Hopf, with
relatively mild-amplitude oscillations (Fig. 1a,b). All three models show
wide-amplitude oscillations though. This can be related to oscillatory
transients in the ODEs, that are excited by noise. In other words, oscillatory
trajectories that would converge to a stable focus in a deterministic setup
(dashed lines in Fig. 1c,d,e,f) are recurrently perturbed under environmental
stochasticity, which favours sustained, broad-amplitude oscillations.

Figure 1: Trajectories of the ODEs (dashed black lines) and SDEs (plain blue
lines) for the 3 functional responses. Top (a-b), Holling; middle (c-d),
Ivlev; bottom (e-f), tanh. Functional response parameters are given in the
text. Other parameters $r=1,\;K=0.5,\;\epsilon=1,\;m=0.1\;,\sigma_{E}=0.25$.

The deterministic and stochastic versions of the bifurcation diagram with
respect to the carrying capacity $K$ reveal that while structural sensitivity
is very present in the deterministic model, it disappears almost entirely in
the stochastic versions (Fig. 2). SDEs with a little less environmental noise
still show some differences in the stochastic bifurcation diagrams (Appendix
A), but these differences are much fainter than in the deterministic model.
These results are also robust to lower values of $\epsilon$ (see code
repository), and the power spectra are similar for the three functional
responses under strong noise (Appendix B).

Figure 2: Bifurcation diagrams for the ODEs (red dots) and SDEs (black dots).
The grey line is the 97.5th percentile. Top (a-b), Holling; middle (c-d),
Ivlev; bottom (e-f), tanh functional responses. Left column, logarithmic
scale; right column, untransformed scale. Functional response parameters are
given in the text. Other parameters
$r=1,\;\epsilon=1,\;m=0.1\;,\sigma_{E}=0.25$.

Despite the apparent increased complexity of stochastic models, these
differences between deterministic and stochastic models can be readily
explained by the eigenvalues of the linearized model at the positive fixed
point. In the case of a focus (oscillatory convergence or divergence to the
fixed point), the solutions of the linearized ODEs and SDEs are proportional
to $\exp(\mu t\pm i\omega t)$ where $\mu$ is the real part and $\omega$ is the
imaginary part of the eigenvalues. When presenting how deterministic
structural sensitivity can arise, Fussmann & Blasius (2005), in their Fig. 2,
zoomed in on the section of $\text{Re}(\lambda)$ values ($\mu$ here, $\tau$ in
Fussmann & Blasius 2005) that crosses zero. In a deterministic and asymptotic
perspective this makes sense: one wants to pinpoint exactly at which $K$ value
the fixed point becomes a limit cycle and vice versa. They noticed that the
zero crossing for $\text{Re}(\lambda)$ was obtained at quite different $K$
values for the various functional forms (reproduced in Fig 3a here, though
Fussmann & Blasius (2005) consider also larger $K$ values). However, in a
stochastic setting we need to think not only of the sign of
$\text{Re}(\lambda)$ but also its magnitude. If we change the y-axis scale of
$\text{Re}(\lambda)$ and, unlike Fussmann & Blasius (2005), decide to consider
all observed $\text{Re}(\lambda)$ values from $K=0$ to $1.5$ (Fig 3b), we see
that (i) $\text{Re}(\lambda)$, even if negative, stops to be a main pulling
force to the equilibrium for relatively low $K$ ($\approx 0.2$, Fig 3b), long
before the bifurcation points. Indeed, the return time after perturbation
being of the order $\frac{1}{\text{Re}(\lambda)}$, we can deduce much longer
return times for $K>0.2$. Considering as well the imaginary part (Fig 3c),
$\text{Im}(\lambda)$ emerges as nonzero for very low $K$ as well ($K<0.2$)—and
nonzero $\text{Im}(\lambda)$ means oscillations. Combining information on both
$\text{Re}(\lambda)$ and $\text{Im}(\lambda)$, we conclude that there is a
broad range of $K$ values below the bifurcation point, for all three
functional responses, where large noise-sustained oscillations are possible
because $\text{Im}(\lambda)$ is nonzero and $\text{Re}(\lambda)$ is very low.

To sum up, the asymptotic limit cycle behaviour, controlled by
$\text{sign}(\text{Re}(\lambda))$, is very sensitive to the functional form
(Fig 3a). However, the magnitude of the real part (Fig 3b), that governs the
response time to perturbations, and the imaginary part (Fig 3c), that governs
the oscillatory tendency, are not that sensitive to $K$. This is what explains
why all three stochastic models start to oscillate in a cyclic manner for very
low $K$ values, with a smooth increase in amplitude as $K$ increases.

Figure 3: Eigenvalues of the nontrivial fixed point for the three functional
responses, as a function of the carrying capacity $K$. Panel (a) shows for
which functional responses the real part of the dominant eigenvalue is
positive, as presented by Fussmann and Blasius (2005); panel (b) zooms out to
show the lowest eigenvalues for very low $K$, the vertical dashed line
separating fast from slow return times to equilibrium; (c) presents the
imaginary part of the eigenvalue. In (c), the vertical dashed line
materializes the $K$ value for which eigenvalues associated to all functional
responses are complex.

Let us note one small caveat: the $\tanh$ functional response predator-prey
model is a bit more difficult to interpret through eigenvalues of the fixed
point, because this assumes a single fixed point and there is a actually
multistability (Fussmann & Blasius, 2005; Seo & Wolkowicz, 2018) with a
coexisting limit cycle for some $K$ values.

In passing, from Fig. 3 it is easy to deduce the average period of the quasi-
cycles, $\text{frequency}=\frac{\omega}{2\pi}\approx\frac{0.3}{2\pi}$, so
$\text{period}=1/\text{frequency}=20$ time units for most $K$ values here.
Overall, the imaginary part of the eigenvalue at the fixed point is as rich in
information as the real part. The actual dominant frequencies in the signal
are slightly below $1/20$ with strong noise (Appendix B).

## Discussion

Although the asymptotic behaviour of the deterministic Rosenzweig-MacArthur
predator-prey model is very sensitive to the type of functional forms used for
the functional response, we have seen that such predictions are radically
altered when strong environmental stochasticity is introduced in all three
versions of this model. The stochastic bifurcations diagrams for similar
Holling, Ivlev, and tanh functional responses were indeed difficult to
differentiate, in stark contrast to the deterministic ones.

While the results may appear at first sight surprising because of the large
difference between stochastic and deterministic model behaviour, such strong
interaction between transients and stochasticity is by no means unexpected. In
fact, the idea has been around for decades (see Pielou 1977, p. 109), together
with the mathematical machinery to analyse it (see Nisbet & Gurney 1982, p.
98). Yet the workhorse of theoretical ecology seems to remain nonlinear
ODEs—given this, it is important to get an interpretation of deterministic
stability analyses that is compatible with a stochastic worldview. Two
interpretation pitfalls when studying with ODEs (and their linearizations)
ecological systems that are actually stochastic are (i) focusing solely on the
sign of the real part of the dominant eigenvalue, rather than its actual value
(e.g., is $\text{Re}(\lambda_{\text{dom}})$ _very_ negative or just below
zero?) and (ii) not reporting the imaginary part of the dominant eigenvalues.
Indeed, as soon as the imaginary part is nonzero the solution oscillates, and
if $\text{Re}(\lambda_{\text{dom}})$ is small the return times will be long,
so that the system may oscillate for a very long time and never reach in
practice a stable equilibrium. A similar conclusion that a too narrow focus on
asymptotic stability of deterministic systems can be misleading can be found
on the abundant literature on ecological transients (see e.g. Hastings 2004;
Hastings _et al._ 2018).

Several authors have noted that stochastic nonlinear dynamical systems tended
to “cross the Hopf bifurcation earlier” (Wiesenfeld, 1985; Neiman _et al._ ,
1997; Gellner _et al._ , 2016; Barraquand _et al._ , 2017; Wyse _et al._ ,
2022), in the sense that the behaviour right after the bifurcation (limit
cycles) is emulated by the forced system right before the bifurcation. This
was the case here as well. These “noisy precursors” are not exactly the same
as the quasi-cycles that would occur in a damped Lotka-Volterra system (where
no limit cycle can be reached), but instead constitute an early transition
towards noisy limit cycle behaviour, induced by the fluctuations in the
bifurcation parameter. In the case of a (supercritical) Hopf, this phenomenon
is not clear-cut because the transition between the excited oscillatory regime
and forced limit cycle regime is relatively smooth (Fig. 2, see also Dutta _et
al._ 2018). In fact, in models with a supercritical Hopf, the transients
become more and more pronounced in the deterministic model as we get near the
bifurcation, so that in a stochastic environment, it becomes difficult to
differentiate between noise-sustained oscillations from limit cycles without a
detailed examination of the power spectrum (Louca & Doebeli, 2014). Other
kinds of bifurcations might make noise-induced transitions more evident, such
as period-doubling (unless noise destroys the periodicity). Overall, it useful
to remember that noise can have effects both relatively far away from
bifurcations, by exciting the transients to a deterministic fixed point, and
effects close to a bifurcation point, by making the system cross the Hopf
bifurcation for different parameter values than in the deterministic setting.

These results, together with the existing body of stochastic population and
community dynamics literature (Nisbet & Gurney, 1982; Nolting & Abbott, 2016;
Arnoldi _et al._ , 2019), some of which pertains to finding signals of regime
transitions (Kéfi _et al._ , 2013; Dutta _et al._ , 2018), can help us to make
suggestions for further work on the structural sensitivity of ecological
models. First, environmental noise effects on population and community
dynamics are usually larger than anticipated, especially in the case of
cycling populations, as recognized early on by Nisbet & Gurney (1976). In case
of sensitivity to structure detected in a deterministic model, a healthy check
is therefore to introduce (multiplicative) noise in various amounts in the
model to verify whether such structural sensitivity is solely a property of
the deterministic model or whether it applies to the stochastic model as well
(Wood & Thomas, 1999). In the simple example with one predator and one prey
that we have considered, structural sensitivity was shown to disappear under
strong environmental noise on the population growth rates. This remains partly
true when introducing noise elsewhere in the model (see Appendix C), although
noise on the functional response’s half-saturation constant provides slightly
less similar bifurcation diagrams than multiplicative noise on the growth
rates. This may explain the difference between this study’s main results and
Wyse _et al._ (2022), who independently found reduced—but not
absent—structural sensitivity in predator-prey models where only the
functional response was stochastic. There are, however, a broad range of
deterministic models that display asymptotic structural sensitivity (Adamson &
Morozov, 2014b). It might well be that some still exhibit structural
sensitivity under strong environmental noise.

Second, although qualitative conclusions with ODE models are completely fine
(e.g., when increasing $K$ or the attack rate, a Hopf bifurcation eventually
occurs), a too heavy focus on the bifurcation structure associated with the
_asymptotic deterministic behaviour of the model_ might distract modellers
from investigating transient behaviours (Hastings, 2004; Hastings _et al._ ,
2018). Studying transients that follow a given perturbation, even in a purely
deterministic model, is a first step towards taking into account the
environmental variation that characterizes real ecosystems (Arnoldi _et al._ ,
2018). Investigating more ecologically meaningful structural sensitivity in a
ODE framework would probably imply developing (or adapting) metrics to compare
transients. Most structural sensitivity approaches focus instead on asymptotic
behaviour, though the possibility of studying transients has been mentioned
(Adamson & Morozov, 2014b). However, the interactions between transients and
noise suggest that SDEs or other random dynamical systems will be paramount to
compare the dynamics generated by alternative model specifications in a
stochastic world.

Properly quantifying the difference in model behaviours generated by SDEs with
different functional forms is challenging though. An idea to quantify more
formally the distance between stochastic attractors could be to use quasi-
potential techniques, introduced in ecology by Nolting & Abbott (2016), and
derived from Freidlin-Wentzell theory. Stochastic differential equations with
the environmental (and thus multiplicative) noise considered above may need to
be converted first to an additive scale through Itô’s lemma (Moore _et al._ ,
2015), but this poses no problem. Nolting & Abbott (2016) note, however, that
the approximation used to produce the quasi-potential works well for small
noise, which differs from the choices made above. A simulation approach might
therefore often be part of the answer to study systems affected by strong
noise. In the (stochastic) bifurcation diagrams showed in this paper, I have
chosen to plot all values of the simulated stochastic process, so that Fig. 2
shows the upper and lower bounds for each $K$ value over a long period of
time, all the values in-between, and a 97.5th percentile (in grey). The
advantage of defining a stochastic bifurcation diagram as such, outside of
sheer simplicity, is that the definition would apply equally well to a chaotic
system with a strange attractor (such as the logistic or Ricker map in
discrete time). Yet, a difference between the deterministic and the stochastic
bifurcation diagram is that as time goes to infinity, the minima and maxima of
the stochastic process still grow increasingly apart—by contrast deterministic
attractors have sharper boundaries. In other words, the tails of the steady-
state distribution of the stochastic process modelled by coupled SDEs are
increasingly revealed as one considers longer time frames. For this reason, it
may be that evaluating in detail the structural sensitivity of coupled
nonlinear SDEs requires to compare the full steady-state distributions under
different model specifications or the moments of such distributions.

### Acknowledgements

I thank the anonymous referee of a previous publication who kindly pointed
towards the literature on noisy precursors of instabilities, and Matthew
Adamson for the discussion in 2018 that motivated me to pursue this work, as
well as for multiple suggestions to improve it. I also thank Hannes Erbis for
an independent replication of some results in 2020, and Frank Hilker and two
referees for comments on the manuscript.

## References

  * Adamson & Morozov (2013) Adamson, M. & Morozov, A.Y. (2013).  When can we trust our model predictions? unearthing structural sensitivity in biological systems.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 469, 20120500.
  * Adamson & Morozov (2014a) Adamson, M. & Morozov, A.Y. (2014a).  Bifurcation analysis of models with uncertain function specification: how should we proceed?  _Bulletin of Mathematical Biology_ , 76, 1218–1240.
  * Adamson & Morozov (2014b) Adamson, M. & Morozov, A.Y. (2014b).  Defining and detecting structural sensitivity in biological models: developing a new framework.  _Journal of Mathematical Biology_ , 69, 1815–1848.
  * Adamson _et al._ (2016) Adamson, M., Morozov, A.Y. & Kuzenkov, O. (2016).  Quantifying uncertainty in partially specified biological models: how can optimal control theory help us?  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 472, 20150627.
  * Aldebert _et al._ (2018) Aldebert, C., Kooi, B.W., Nerini, D. & Poggiale, J.C. (2018).  Is structural sensitivity a problem of oversimplified biological models? Insights from nested Dynamic Energy Budget models.  _Journal of Theoretical Biology_ , 448, 1–8.
  * Arnoldi _et al._ (2018) Arnoldi, J.F., Bideault, A., Loreau, M. & Haegeman, B. (2018).  How ecosystems recover from pulse perturbations: A theory of short-to long-term responses.  _Journal of Theoretical Biology_ , 436, 79–92.
  * Arnoldi _et al._ (2019) Arnoldi, J.F., Loreau, M. & Haegeman, B. (2019).  The inherent multidimensionality of temporal variability: how common and rare species shape stability patterns.  _Ecology Letters_ , 22, 1557–1567.
  * Barraquand & Gimenez (2021) Barraquand, F. & Gimenez, O. (2021).  Fitting stochastic predator–prey models using both population density and kill rate data.  _Theoretical Population Biology_ , 138, 1–27.
  * Barraquand _et al._ (2017) Barraquand, F., Louca, S., Abbott, K.C., Cobbold, C.A., Cordoleani, F., DeAngelis, D.L., Elderd, B.D., Fox, J.W., Greenwood, P., Hilker, F.M. _et al._ (2017).  Moving forward in circles: challenges and opportunities in modelling population cycles.  _Ecology Letters_ , 20, 1074–1092.
  * Dutta _et al._ (2018) Dutta, P.S., Sharma, Y. & Abbott, K.C. (2018).  Robustness of early warning signals for catastrophic and non-catastrophic transitions.  _Oikos_ , 127, 1251–1263.
  * Fussmann & Blasius (2005) Fussmann, G.F. & Blasius, B. (2005).  Community response to enrichment is highly sensitive to model structure.  _Biology Letters_ , 1, 9–12.
  * Gellner _et al._ (2016) Gellner, G., McCann, K.S. & Hastings, A. (2016).  The duality of stability: towards a stochastic theory of species interactions.  _Theoretical Ecology_ , 9, 477–485.
  * Hastings (2004) Hastings, A. (2004).  Transients: the key to long-term ecological understanding?  _Trends in Ecology & Evolution_, 19, 39–45.
  * Hastings _et al._ (2018) Hastings, A., Abbott, K.C., Cuddington, K., Francis, T., Gellner, G., Lai, Y.C., Morozov, A., Petrovskii, S., Scranton, K. & Zeeman, M.L. (2018).  Transient phenomena in ecology.  _Science_ , 361, eaat6412.
  * Kéfi _et al._ (2013) Kéfi, S., Dakos, V., Scheffer, M., Van Nes, E.H. & Rietkerk, M. (2013).  Early warning signals also precede non-catastrophic transitions.  _Oikos_ , 122, 641–648.
  * Lande _et al._ (2003) Lande, R., Engen, S. & Saether, B.E. (2003).  _Stochastic population dynamics in ecology and conservation_.  Oxford University Press.
  * Louca & Doebeli (2014) Louca, S. & Doebeli, M. (2014).  Distinguishing intrinsic limit cycles from forced oscillations in ecological time series.  _Theoretical Ecology_ , 7, 381–390.
  * McKane & Newman (2004) McKane, A. & Newman, T. (2004).  Stochastic models in population biology and their deterministic analogs.  _Physical Review E_ , 70, 41902.
  * McKane & Newman (2005) McKane, A. & Newman, T. (2005).  Predator-prey cycles from resonant amplification of demographic stochasticity.  _Physical Review Letters_ , 94, 218102.
  * Moore _et al._ (2015) Moore, C.M., Stieha, C.R., Nolting, B.C., Cameron, M.K. & Abbott, K.C. (2015).  Qpot: an r package for stochastic differential equation quasi-potential analysis.  _arXiv preprint arXiv:1510.07992_.
  * Mutshinda _et al._ (2009) Mutshinda, C.M., O’Hara, R.B. & Woiwod, I.P. (2009).  What drives community dynamics?  _Proceedings of the Royal Society B: Biological Sciences_ , 276, 2923–2929.
  * Neiman _et al._ (1997) Neiman, A., Saparin, P.I. & Stone, L. (1997).  Coherence resonance at noisy precursors of bifurcations in nonlinear dynamical systems.  _Physical Review E_ , 56, 270.
  * Nisbet & Gurney (1976) Nisbet, R. & Gurney, W. (1976).  A simple mechanism for population cycles.  _Nature_ , 263, 319–320.
  * Nisbet & Gurney (1982) Nisbet, R. & Gurney, W. (1982).  _Modelling fluctuating populations_.  John Wiley & Sons.
  * Nolting & Abbott (2016) Nolting, B.C. & Abbott, K.C. (2016).  Balls, cups, and quasi-potentials: quantifying stability in stochastic systems.  _Ecology_ , 97, 850–864.
  * Pielou (1977) Pielou, E. (1977).  _Mathematical ecology.[2d ed.]._ 2nd edn.  John Wiley and Sons, New York, USA.
  * Pineda-Krch _et al._ (2007) Pineda-Krch, M., Blok, J., Dieckmann, U. & Doebeli, M. (2007).  A tale of two cycles-distinguishing quasi-cycles and limit cycles in finite predator-prey populations.  _Oikos_ , 116, 53.
  * Rosenzweig (1971) Rosenzweig, M. (1971).  Paradox of enrichment: destabilization of exploitation ecosystems in ecological time.  _Science_ , 171, 385–387.
  * Rosenzweig & MacArthur (1963) Rosenzweig, M. & MacArthur, R. (1963).  Graphical representation and stability conditions of predator-prey interactions.  _American Naturalist_ , 97, 209–223.
  * Seo & Wolkowicz (2018) Seo, G. & Wolkowicz, G.S. (2018).  Sensitivity of the dynamics of the general rosenzweig–macarthur model to the mathematical form of the functional response: a bifurcation theory approach.  _Journal of Mathematical Biology_ , 76, 1873–1906.
  * Wiesenfeld (1985) Wiesenfeld, K. (1985).  Noisy precursors of nonlinear instabilities.  _Journal of Statistical Physics_ , 38, 1071–1097.
  * Wood & Thomas (1999) Wood, S.N. & Thomas, M.B. (1999).  Super–sensitivity to structure in biological models.  _Proceedings of the Royal Society of London. Series B: Biological Sciences_ , 266, 565–570.
  * Wyse _et al._ (2022) Wyse, S.K., Martignoni, M.M., Mata, M.A., Foxall, E. & Tyson, R.C. (2022).  Structural sensitivity in the functional responses of predator–prey models.  _Ecological Complexity_ , 51, 101014.

##  Appendix A Additional bifurcation diagrams

We first zoom out to a large parameter space where $K$ is allowed to go up to
15 (Fig. A1) as in Fussmann & Blasius (2005), with the same variance of the
noise as in the main text. The results are similar to those of the main text,
that is, no sensitivity to model structure in the stochastic model.

Figure A1: Bifurcation diagrams for the ODEs (red dots) and SDEs (black dots),
up to large $K$, $\sigma_{E}=0.25$. Top (a-b), Holling; middle (b-c), Ivlev;
bottom (e-f), tanh functional responses. Some extinction occurred in the
Holling simulations for large K, hence the small gap; this might be remedied
with finer-scale integration.

Now we go back to small values of $K$ as in Fig. 2 of the main text, but we
reduce the amount of stochasticity from $\sigma_{E}=0.25$ to $\sigma_{E}=0.1$
(Fig. A2). We see that the main result does hold—there is a smooth transition
to large amplitude cycles for all three functional responses—although cycle
amplitudes now differ a little between the three functional responses,
especially with tanh. Still structural sensitivity has massively decreased
with the introduction of stochasticity.

Figure A2: Bifurcation diagrams for the ODEs (red dots) and SDEs (black dots),
$\sigma_{E}=0.1$. Top (a-b), Holling; middle (b-c), Ivlev; bottom (e-f), tanh
functional responses. Other parameters otherwise identical to those of the
main text.

##  Appendix B Power spectrum

We present here the power spectrum for the time series of Fig. 1, in a raw
periodogram (Fig. B1) and smoothed (Fig. B2) format. We see that the spectrum
does not appear to be sensitive to details of the representation chosen.
Although some differences appear in the spectrum for different values of $K$,
no large differences between the three functional responses are noticeable
(small differences at low frequencies should not be interpreted due to
difficulties in smoothing the spectrum close to zero). Other techniques to
reconstruct the Fourier spectrum (e.g. Welch, AR($p$)) tend to produce redder
spectra and in some cases the dominant frequencies disappear, but (1) visual
inspection of spacings between local maxima of the time series match the
dominant frequencies presented here, which are below 1/20, the frequency of
quasi-cycles suggested by the linearized model, and (2) computation of the
autocorrelation (logically) shows the same frequencies. A lower noise variance
gets the dominant frequency closer to 1/20.

Figure B1: Power spectrum of $N$. Raw periodogram computed using the FFT with
the function spectrum in R. Figure B2: Power spectrum of $N$. Smoothed using a
modified Daniell(11,11,11) kernel.

##  Appendix C Large noise on the functional response

Here we make a stochastic functional response by introducing log-normal,
autocorrelated noise on the half-saturation constant, a trick previously used
by Barraquand & Gimenez (2021) in a discrete-time context. The model now
writes

$\displaystyle dX=$
$\displaystyle\left(g(X)-f(X,Z)Y\right)dt+\sigma_{E}XdW_{x}$ (C1)
$\displaystyle dY=$ $\displaystyle\left(\epsilon
f(X,Z)Y-mY\right)dt+\sigma_{E}YdW_{y}$ (C2) $\displaystyle dZ=$
$\displaystyle-\gamma Zdt+\sigma_{\text{FR}}dW_{z}$ (C3)

where the Holling type functional response is now

$f_{H}(x,z)=\frac{Cx}{De^{z}+x}$ (C4)

and the Ivlev-type

$f_{I}(x,z)=C\left(1-\exp\left(-\frac{\ln(2)}{De^{z}}x\right)\right).$ (C5)

The noisy functional responses generate time series shown in Fig. C1. We used
$\sigma_{E}=0.1$ and $\sigma_{\text{FR}}=0.25$ to model moderate noise on the
growth rates but large multiplicative noise on the functional response’s half-
saturation constant. We used $C=1$ and $D=\frac{0.5}{\ln(2)}$ for both
functional responses.

Figure C1: Trajectories of the ODEs (dashed black lines) and SDEs (plain blue
lines) for 2 stochastic versions of the functional responses. Top (a-b),
Holling; bottom (b-d), Ivlev. $K=0.5$.

We only compare Holling and Ivlev functional responses for which the half-
saturation is clearly defined as a parameter. Because it can be logical to
assume that the functional response’s half-saturation constant will not be
fully white noise (it will vary with autocorrelated factors like other
species’ dynamics), we used an Ornstein-Uhlenbeck process, with parameter
$\gamma=0.1$. Noise on the functional response seem to make the dynamical
system shift between alternate dynamical states, with small vs large
oscillations, which might warrant further exploration. Although there is a
little more amplitude in the Holling-type model, structural sensitivity has
again massively decreased.

Figure C2: Bifurcation diagrams for the ODEs (red dots) and SDEs (black dots),
$\sigma_{E}=0.1$ and $\sigma_{\text{FR}}=0.25$. Top (a-b), Holling; bottom
(b-d), Ivlev functional response.