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

Frédéric Barraquand

► To cite this version:

Frédéric Barraquand. No sensitivity to functional forms in the Rosenzweig–MacArthur model with strong environmental stochasticity. Journal of Theoretical Biology, 2023, 572, pp.111566. 10.1016/j.jtbi.2023.111566. hal-04293523

HAL Id: hal-04293523

https://hal.science/hal-04293523v1

Submitted on 2 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# 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 fluctuation 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 frederic.barraquand@u-bordeaux.fr

Published in *Journal of Theoretical Biology* with doi:10.1016/j.jtbi.2023.111566
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, 2014b; Adamson et al.,
2016). Several mathematical approaches have emerged to both detect such structural sensitivity of nonlinear
ordinary differential equations (Adamson & Morozov, 2013, 2014a,b), 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 stand-
point. 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 envi-
ronment. 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 pertur-
bations 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 determin-
istic 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 phenomenon 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:

$$dX = (g(X) - f(X)Y) dt + \sigma_E X dW_x \quad (1)$$

$$dY = (\epsilon f(X)Y - mY) dt + \sigma_E Y dW_y. \quad (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 (1 - \exp(-b_I x))$, $f_T(x) = a_T \tanh(b_T x)$ and a logistic prey growth rate function $g(x) = rx(1 - \frac{x}{K})$. 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><img src="image_2.png" /><figcaption>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$.</figcaption></figure>

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><img src="image_3.png" /><figcaption>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$.</figcaption></figure>

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 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 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 Re(λ) but also its magnitude. If we change the y-axis scale of Re(λ) and, unlike Fussmann & Blasius (2005), decide to consider all observed Re(λ) values from *K* = 0 to 1.5 (Fig 3b), we see that (i) Re(λ), even if negative, stops to be a main pulling force to the equilibrium for relatively low *K* (≈ 0.2, Fig 3b), long before the bifurcation points. Indeed, the return time after perturbation being of the order $1/\text{Re}(\lambda)$, we can deduce much longer return times for *K* > 0.2. Considering as well the imaginary part (Fig 3c), Im(λ) emerges as nonzero for very low *K* as well (*K* < 0.2)—and nonzero Im(λ) means oscillations. Combining information on both Re(λ) and Im(λ), 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 Im(λ) is nonzero and Re(λ) is very low.

To sum up, the asymptotic limit cycle behaviour, controlled by sign(Re(λ)), 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.

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, frequency $= \frac{\omega}{2\pi} \approx \frac{0.3}{2\pi}$, so period $= 1/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
<figure><img src="image_4.png" /><figcaption>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.</figcaption></figure>

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 Re($\lambda_{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 Re($\lambda_{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 fluctuations in a 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 asso-
ciated 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 (see e.g. 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, 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, 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, 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, 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, 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, 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, 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, F. & Gimenez, O. (2021). Fitting stochastic predator–prey models using both population density and kill rate data. *Theoretical Population Biology*, 138, 1–27.

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, P.S., Sharma, Y. & Abbott, K.C. (2018). Robustness of early warning signals for catastrophic and non-catastrophic transitions. *Oikos*, 127, 1251–1263.

Fussmann, G.F. & Blasius, B. (2005). Community response to enrichment is highly sensitive to model structure. *Biology Letters*, 1, 9–12.

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, A. (2004). Transients: the key to long-term ecological understanding? *Trends in Ecology & Evolution*, 19, 39–45.

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, 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, R., Engen, S. & Saether, B.E. (2003). *Stochastic population dynamics in ecology and conservation*. Oxford University Press.

Louca, S. & Doebeli, M. (2014). Distinguishing intrinsic limit cycles from forced oscillations in ecological time series. *Theoretical Ecology*, 7, 381–390.

McKane, A. & Newman, T. (2004). Stochastic models in population biology and their deterministic analogs. *Physical Review E*, 70, 41902.

McKane, A. & Newman, T. (2005). Predator-prey cycles from resonant amplification of demographic stochasticity. *Physical Review Letters*, 94, 218102.

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, 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, A., Saparin, P.I. & Stone, L. (1997). Coherence resonance at noisy precursors of bifurcations in nonlinear dynamical systems. *Physical Review E*, 56, 270.

Nisbet, R. & Gurney, W. (1976). A simple mechanism for population cycles. *Nature*, 263, 319–320.

Nisbet, R. & Gurney, W. (1982). *Modelling fluctuating populations*. John Wiley & Sons.

Nolting, B.C. & Abbott, K.C. (2016). Balls, cups, and quasi-potentials: quantifying stability in stochastic systems. *Ecology*, 97, 850–864.

Pielou, E. (1977). *Mathematical ecology.* [2d ed.]. 2nd edn. John Wiley and Sons, New York, USA.

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, M. (1971). Paradox of enrichment: destabilization of exploitation ecosystems in ecological time. *Science*, 171, 385–387.
Rosenzweig, M. & MacArthur, R. (1963). Graphical representation and stability conditions of predator-prey interactions. *American Naturalist*, 97, 209-223.

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, K. (1985). Noisy precursors of nonlinear instabilities. *Journal of Statistical Physics*, 38, 1071–1097.

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, 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><img src="image_5.png" /><figcaption>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.</figcaption></figure>

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><img src="image_6.png" /><figcaption>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.</figcaption></figure>
# 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><img src="image_7.png" /><figcaption>Figure B1: Power spectrum of *N*. Raw periodogram computed using the FFT with the function *spectrum* in R.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure B2: Power spectrum of *N*. Smoothed using a modified Daniell(11,11,11) kernel.</figcaption></figure>
# 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

$$dX = (g(X) - f(X, Z)Y) dt + \sigma_E X dW_x \quad (C1)$$

$$dY = (\epsilon f(X, Z)Y - mY) dt + \sigma_E Y dW_y \quad (C2)$$

$$dZ = -\gamma Z dt + \sigma_{FR} dW_z \quad (C3)$$

where the Holling type functional response is now

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

and the Ivlev-type

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

The noisy functional responses generate time series shown in Fig. C1. We used $\sigma_E = 0.1$ and $\sigma_{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><img src="image_9.png" /><figcaption>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.</figcaption></figure>

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 (see also Adamson & Morozov (2014b) on this point). Although there is a little
more amplitude in the Holling-type model, structural sensitivity has again massively decreased (see Fig. C2).
<figure><img src="image_10.png" /><figcaption>Figure C2: Bifurcation diagrams for the ODEs (red dots) and SDEs (black dots), $\sigma_E = 0.1$ and $\sigma_{FR} = 0.25$. Top (a-b), Holling; bottom (b-d), Ivlev functional response.</figcaption></figure>