# General statistical scaling laws for stability in ecological systems

Adam Thomas Clark, Jean-francois Arnoldi, Yuval R Zelnik, György Barabas, Dorothee Hodapp, Canan Karakoç, Sara König, Viktoriia Radchuk, Ian Donohue, Andreas Huth, et al.

► To cite this version:

Adam Thomas Clark, Jean-francois Arnoldi, Yuval R Zelnik, György Barabas, Dorothee Hodapp, et al.. General statistical scaling laws for stability in ecological systems. Ecology Letters, 2021, 24 (7), pp.1474-1486. 10.1111/ele.13760. hal-03574740

HAL Id: hal-03574740

https://ut3-toulouseinp.hal.science/hal-03574740v1

Submitted on 17 Oct 2024

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

# General statistical scaling laws for stability in ecological systems

Adam Thomas Clark<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> ID

Jean-François Arnoldi<fnref n="4" /> ID

Yuval R. Zelnik<fnref n="5" />,<fnref n="6" />

György Barabas<fnref n="7" />,<fnref n="8" />

Dorothee Hodapp<fnref n="9" />,<fnref n="10" />

Canan Karakoç<fnref n="3" />,<fnref n="11" />

Sara König<fnref n="12" />

Viktoriia Radchuk<fnref n="13" /> ID

Ian Donohue<fnref n="4" /> ID

Andreas Huth<fnref n="14" />

Claire Jacquet<fnref n="15" />,<fnref n="16" /> ID

Claire de Mazancourt<fnref n="6" />

Andrea Mentges<fnref n="3" />,<fnref n="17" />

Dorian Nothaaß<fnref n="1" />,<fnref n="14" />

Lauren G. Shoemaker<fnref n="18" /> ID

Franziska Taubert<fnref n="14" />

Thorsten Wiegand<fnref n="3" />,<fnref n="14" />

Shaopeng Wang<fnref n="19" /> ID

Jonathan M. Chase<fnref n="3" />,<fnref n="17" /> ID

Michel Loreau<fnref n="6" />

Stanley Harpole<fnref n="1" />,<fnref n="3" />,<fnref n="20" />

**Correspondence**

Adam Thomas Clark, Department of Physiological Diversity, Helmholtz Centre for Environmental Research (UFZ), Leipzig, Germany.

Email: adam.clark@uni-graz.at

**Funding information**

TULIP French Laboratory of Excellence, Grant/Award Number: ANR-10-LABX-41 and ANR-11-IDEX-0002-02; Swedish Research Council, Grant/Award Number: VR 2017-05245; Deutsche Forschungsgemeinschaft, Grant/Award Number: FZT 118; NSF, Grant/Award Number: EPS-1655726; Horizon 2020 Framework Programme, Grant/Award Number: 666971

**Abstract**

Ecological stability refers to a family of concepts used to describe how systems of interacting species vary through time and respond to disturbances. Because observed ecological stability depends on sampling scales and environmental context, it is notoriously difficult to compare measurements across sites and systems. Here, we apply stochastic dynamical systems theory to derive general statistical scaling relationships across time, space, and ecological level of organisation for three fundamental stability aspects: resilience, resistance, and invariance. These relationships can be calibrated using random or representative samples measured at individual scales, and projected to predict average stability at other scales across

<footnote n="1">¹Department of Physiological Diversity, Helmholtz Centre for Environmental Research (UFZ), Leipzig, Germany</footnote>

<footnote n="2">²Institute of Biology, University of Graz, Graz, Austria</footnote>

<footnote n="3">³German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig, Leipzig, Germany</footnote>

<footnote n="4">⁴Zoology Department, Trinity College Dublin, Dublin, Ireland</footnote>

<footnote n="5">⁵Department of Ecology, Swedish University of Agricultural Sciences, Uppsala, Sweden</footnote>

<footnote n="6">⁶Centre for Biodiversity Theory and Modelling, Theoretical and Experimental Ecology Station, CNRS, Moulis, France</footnote>

<footnote n="7">⁷Division of Theoretical Biology, Department of Physics, Chemistry, and Biology, Linköping University, Linköping, Sweden</footnote>

<footnote n="8">⁸MTA-ELTE Theoretical Biology and Evolutionary Ecology Research Group, Budapest, Hungary</footnote>

<footnote n="9">⁹Helmholtz Institute for Functional Marine Biodiversity (HIFMB), Oldenburg, Germany</footnote>

<footnote n="10">¹⁰Alfred-Wegener-Institute Helmholtz-Centre for Polar and Marine Research (AWI), Bremerhaven, Germany</footnote>

<footnote n="11">¹¹Department of Environmental Microbiology, Helmholtz Centre for Environmental Research (UFZ), Leipzig, Germany</footnote>

<footnote n="12">¹²Department of Soil System Science, Helmholtz Centre for Environmental Research (UFZ), Halle (Saale), Germany</footnote>

<footnote n="13">¹³Department of Ecological Dynamics, Leibniz Institute for Zoo and Wildlife Research (IZW), Berlin, Germany</footnote>

<footnote n="14">¹⁴Department of Ecological Modelling, Helmholtz Centre for Environmental Research (UFZ), Leipzig, Germany</footnote>

<footnote n="15">¹⁵Department of Evolutionary Biology and Environmental Studies, University of Zurich, Zürich, Switzerland</footnote>

<footnote n="16">¹⁶Department of Aquatic Ecology, Swiss Federal Institute of Aquatic Science and Technology, Eawag, Dübendorf, Switzerland</footnote>

<footnote n="17">¹⁷Department of Computer Sciences, Martin Luther University, Halle, Germany</footnote>

<footnote n="18">¹⁸Botany Department, University of Wyoming, Laramie, WY, USA</footnote>

<footnote n="19">¹⁹Institute of Ecology, College of Urban and Environmental Science, and Key Laboratory for Earth Surface Processes of the Ministry of Education, Peking University, Beijing, China</footnote>

<footnote n="20">²⁰Institute of Biology, Martin Luther University, Halle, Germany</footnote>

<footnote>This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.</footnote>

<footnote>© 2021 The Authors. *Ecology Letters* published by John Wiley & Sons Ltd.</footnote>
Editor: Tim Coulson

a wide range of contexts. Moreover deviations between observed vs. extrapolated scaling relationships can reveal information about unobserved heterogeneity across time, space, or species. We anticipate that these methods will be useful for cross-study synthesis of stability data, extrapolating measurements to unobserved scales, and identifying underlying causes and consequences of heterogeneity.

### KEYWORDS

community, disturbance, diversity, invariability, invariance, population, resilience, resistance, spatial, temporal

## INTRODUCTION

Ecological stability refers to a range of concepts that describe how interacting systems of species and their environments vary over time (Donohue et al., 2013; Grimm & Wissel, 1997; Holling, 1973; Ives & Carpenter, 2007; Lewontin, 1969; May, 1973; Pimm, 1984; Pimm et al., 2019). Stability metrics are primary tools used to study how systems withstand and recover from disturbances, and are therefore vital for predicting how anthropogenic pressures, such as land use change, global warming, and species loss are likely to influence natural systems (Carpenter et al., 2001; Donohue et al., 2016; Scheffer et al., 2009; Zelnik et al., 2018). However, measures of stability are highly context dependent, and vary with the spatiotemporal scale and ecological level of organisation at which measurements are made (Clark et al., 2019; Domínguez-García et al., 2019; Kéfi et al., 2019; Levin, 1992; Stommel, 1963). Measurements taken across different scales are therefore not directly comparable, which limits opportunities for cross-system comparison and synthesis (Clark et al., 2019; Csillag et al., 2000; Levin, 1992; Wang et al., 2019). This scale dependence is also a major challenge for conservation and management, because the scales that are most relevant for decision making often differ from those at which ecological systems are measured (Carpenter et al., 2001; Chesson, 2000; Isbell et al., 2018; Leibold & Chase, 2018; Levin, 1992).

To address this challenge, our goal here is to identify general statistical scaling laws for stability in ecological systems. Specifically, we seek to identify relationships that can: (1) be fitted using empirical data observed at one set of scales; (2) be extrapolated to accurately describe conditions across other scales; and (3) that are valid regardless of the underlying processes governing dynamics.

Throughout, we focus on three common stability metrics—resilience, resistance, and invariance—and explore how these change across levels of temporal, spatial, and ecological organisation (see *Stability and scale definitions* below) (Donohue et al., 2016; Grimm & Wissel, 1997). These metrics are especially relevant to ecological studies because they describe how systems respond to regularly occurring, externally imposed pulse perturbations, as might be expected from, for example,

extreme weather events or anthropogenic disturbances such as logging, mowing, or fishing (Holling, 1973; Isbell et al., 2018; Ives & Carpenter, 2007; Zelnik et al., 2018). Moreover recent advances in stochastic dynamical systems theory have identified links between these stability metrics and basic statistical properties that can be estimated from empirical data (Arnoldi et al., 2018; Haegeman et al., 2016; Lee et al., 2020). Because these statistical properties arise directly from the mathematical definitions of variance and covariance, the resulting scaling relationships are exceedingly general (Arnoldi et al., 2019; Wang et al., 2017; Zelnik et al., 2018).

We proceed in four parts. First, we present our working definitions for resilience, resistance, and invariance, and for temporal, spatial, and ecological scale. Second, we derive the scaling relationships, and explain how they arise from the basic mathematical definitions of variance and covariance. Third, we introduce three models that simulate species abundance dynamics across a range of scales, and use these to demonstrate the scaling functions. Finally, we introduce methods for fitting the scaling relationships to observational data. These methods can be implemented via our accompanying *ecostatscale* R package, available via CRAN and archived on GitHub (github.com/adamtclark/ecostatscale) and Zenodo (https://doi.org/10.5281/zenodo.4626672).

## METHODS

The statistical scaling functions developed here are rooted in stochastic dynamical systems theory, which extends the general concept of dynamical systems modelling (e.g. ordinary differential equations, such as in logistic growth or Lotka–Volterra interactions) to include effects of stochastic perturbations on dynamic variables (Black & McKane, 2012; Gillespie, 1992; Wilkinson, 2011, 2018). In particular, we leverage the fact that when variability is grouped across different subsets of data, it adheres to rules that follow directly from the mathematical definitions of variance and covariance. These rules hold for any set of variables regardless of their distributions, so long as their mean, variance, and covariance exist and can be measured.
As an example, imagine a set of n timeseries {$y_1, y_2, \dots, y_n$}, each representing a variable with its own mean and variance. For our purposes, each $y_i$ might represent abundance dynamics for different species. By the definition of covariance, the variance of the summed abundance across all species is

$$ \begin{aligned} \mathrm{var}(\Sigma_i y_i) &= \Sigma_i (\mathrm{var}(y_i) + \Sigma_{j \neq i} \mathrm{cov}(y_i, y_j)) \\ &= n \langle \mathrm{var}(y_i) \rangle + (n^2 - n) \langle \mathrm{cov}(y_i, y_j) \rangle \end{aligned} \quad (1a) $$

where $\langle \mathrm{var}(y_i) \rangle$ and $\langle \mathrm{cov}(y_i, y_j) \rangle$ indicate, respectively, average species-level variance and between-species covariance. Following Wang et al. (2017), this relationship can be written as

$$ \mathrm{var}(\Sigma_i y_i) = n \langle \mathrm{var}(y_i) \rangle (1 + (n-1) \langle \rho(y_i, y_j) \rangle) \quad (1b) $$

where $\langle \rho(y_i, y_j) \rangle = \langle \mathrm{cov}(y_i, y_j) \rangle / \sqrt{\langle \mathrm{var}(y_i) \rangle \langle \mathrm{var}(y_j) \rangle}$, which is roughly analogous to the Pearson correlation coefficient (see Appendix A for more details). Critically, the relationships in Equations 1a–1b hold regardless of the underlying structure of the groups, or the processes by which data were generated.

## Stability and scale definitions

Our next step is to define stability and scale in terms that can be related to the statistical rules in Equations 1a–1b. Specifically, we assume that species abundances, and the effects of disturbances on those abundances, can both be defined in terms of statistical distributions with a measurable mean and variance. In practice, this assumption implies that: (1) species abundances fluctuate around a fixed average value; (2) disturbances can be abstracted into discrete pulses that lead to near-instantaneous

changes in species abundances; and (3) disturbances do not alter the processes underlying abundance dynamics (e.g. reproductive rates). We will eventually relax some of these assumptions, though even in their strict form, they apply broadly across abundance dynamics that lack strong monotonic trends, and to disturbances that act quickly relative to species growth rates.

We choose to focus on resilience, resistance, and invariance both because they are widely used in ecology (Donohue et al., 2013; Grimm & Wissel, 1997; Ives & Carpenter, 2007; Pimm, 1984; Pimm et al., 2019), and because they are particularly well suited to statistical abstraction (see examples in Figure 1; Arnoldi et al., 2018; Lande 1993***; Turelli 1986***). Resilience (r) describes the per-capita rate at which abundances recover following a disturbance – high resilience thus implies fast recovery. Given N(t) as abundance at time t, we can therefore calculate r = log(N(t + τ)/N(t))/τ for a sufficiently small unit of time τ. Resistance (σ⁻¹) is inversely related to the effect that disturbances have on abundances—thus, high resistance implies either that disturbances are weak, or that they weakly impact species. To quantify disturbance impacts, we define σ = sd(N(t_i) - N(t_i - τ)), where t_i is the moment of time at which a disturbance occurred, t_i - τ is the time immediately before that disturbance, and sd is the standard deviation. Invariance (var(x)⁻¹) is inversely related to how widely abundances fluctuate over time—high invariance thus implies that abundances remain relatively constant, and var(x) is simply temporal variance calculated across timeseries x. To avoid confusion, we use the term “invariance” rather than the more common “invariability,” as the latter often signifies a specific metric standardised by mean abundance.

We focus on three types of scaling, chosen both because they are widely used in empirical studies, and because they

<figure><img src="image_3.png" /><figcaption>**FIGURE 1** Working definitions of different aspects of stability for this paper. Black line shows dynamics for a basic linear process with equilibrium at x = 0 and subjected to repeated disturbance events, as introduced in the single-patch, single-species model in Equation 3. Disturbance impacts on x are drawn from a normal distribution with mean μ = 0 and standard deviation σ, and disturbances occur with average frequency λ. In this framework, resilience (r) describes rate of return towards equilibrium following a disturbance; resistance (σ⁻¹) describes ability to avoid displacement by disturbances; and invariance (var(x)⁻¹) describes the tendency to stay near the equilibrium value (blue dashed lines show standard deviation of the abundance fluctuations). Here, and in subsequent figures, we use resistance⁻¹ and invariance⁻¹ to signify these inverse relationships</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIGURE 2 Working definitions of temporal, spatial, and ecological scales for this paper. Black shaded regions and lines highlight individual measurements at each scale. (a–c) Temporal scale describes the sampling interval (e.g. daily vs. annual sampling). (d–f) Spatial scale describes the total area across which available information is aggregated (e.g. plot size). (g–i) Ecological scale describes the number of species across which information is combined (e.g. species vs. functional groups). Note that although the scale of measurements varies among panels, the total simulated extent is constant within each type of scaling scenario (i.e. total time span, spatial area, or community size). For easier visualisation, results in panels (d–i) are standardised by the area or the number of species sampled (a and m, respectively)</figcaption></figure>

relate to grouping data into different subsets, which lends itself to the definitions in Equations 1a–1b (see examples in Figure 2). **Temporal scale** describes the rate at which observations are made over time (e.g. measurements per year). Thus, increasing temporal scale leads to larger time gaps between observations. **Spatial scale** describes the grain at which spatial information is aggregated (e.g. plot size, sample volume). Thus, increasing spatial scale involves grouping information across a larger number of replicates. **Ecological scale** describes the number of species across which measurements are aggregated (e.g. species vs. functional groups). Thus, increasing ecological scale involves grouping information across a larger number of species.

## Statistical scaling relationships

Because they are primarily related to how data are grouped across observations, spatial and ecological scales have similar relationships with stability, which

follow directly from Equations 1a–1b, and are described below. Conversely, because our definitions of stability are themselves functions of time, temporal scale has somewhat different effects, primarily related to measurement biases. These are described in the *Measuring stability* section, below.

The simplest scaling relationship to arise from Equations 1a–1b is that for invariance. Imagine that we have a series of observations of invariance measured at scale *b*, and wish to extrapolate these to some larger aggregate scale *B*. For example, *B* might represent a landscape (e.g. Figure 2f), and *b* could represent individual plots within that landscape (e.g. Figure 2d). Substituting var(*y*) with plot-level variance in Equation 1b yields

$$ \text{var}(x_B) = \langle \text{var}(x_b) \rangle (B/b) (1 + \langle \rho(x_b) \rangle (B/b - 1)) \quad (2a) $$

where $\langle \text{var}(x_b) \rangle$ and $\langle \rho(x_b) \rangle$ represent, respectively, the average variance and average pairwise correlation measured at scale *b*.
A similar relationship arises for resistance. If we replace var(y) in Equation 1b with $\sigma^2$, that is, the observed variance in disturbance impacts on abundances, we find

$$\sigma_B = \langle \sigma_b \rangle = \sqrt{(B/b) (1 + \langle \rho(\xi_b) \rangle (B/b - 1))} \quad (2b)$$

where $\langle \sigma_b \rangle$ is the mean standard deviation of disturbance impacts observed at scale b, and $\langle \rho(\xi_b) \rangle$ is the average pairwise correlation between observed disturbance impacts. For example, in a plant community, $\sigma_b$ might represent the average distance that extreme weather events push species away from their average abundances, whereas $\rho(\xi_b)$ would indicate whether species respond similarly ($\rho(\xi_b) > 0$) or differently ($\rho(\xi_b) < 0$) to those events.

Finally, resilience scaling relationships are somewhat more complicated. Given strong effects of dispersal or species interactions, r can vary greatly across space and time, even within a single measurement scale (Arnoldi et al., 2018). While the precise relationship varies with context (see Appendix B), the median expectation for an aggregate spatial or ecological scale B (e.g. total landscape or total community abundance) can be approximated as

$$\begin{aligned} \text{median} (r_B) & \approx (\lambda/2) (\langle \sigma_b^2 \rangle (1 + \langle \rho(\xi_b) \rangle (B/b - 1))) / (\langle \text{var}(x_b) \rangle (1 + \langle \rho(x_b) \rangle (B/b - 1))) \\ & \approx (\lambda/2) (\sigma_B^2) / \text{var}(x_B) \end{aligned} \quad (2c)$$

where $\lambda$ is average disturbance frequency. Note that Equation 2c is effectively the ratio of the variance vs. resistance scaling relationships in Equations 2a–2b. Unlike Equations 2a–2b, the relationship in Equation 2c formally applies only to systems with linear dynamics (see Appendix F). However, in practice, the relationship tends to hold across more complex systems (Arnoldi et al., 2019; Zelnik et al., 2018), as we will demonstrate in the modelling tests below.

## Modelling tests

To illustrate the statistical scaling relationships presented above, we consider three model structures, describing: (1) a single species occupying a single patch; (2) a single species occupying multiple patches; and (3) multiple species occupying a single patch. These roughly correspond to the simplest structures needed to explore temporal, spatial, and ecological scaling, respectively. Although less general than the statistical scaling functions themselves, we choose these basic models because they allow us to analytically relate expected scaling relationships to simple functions of model parameters, against which the predictions of the statistical scaling functions can be clearly compared. Detailed model derivations, along with comparisons to the classic logistic

and Lotka-Volterra models, are available in Appendices B and C.

Our first model of a single species and patch describes a linear stochastic dynamical system near equilibrium (Arnoldi et al., 2018, 2019; Lande et al., 1999; Lee et al., 2020). We consider dynamics of standardised abundance $x(t) = N(t) - K$, where $N(t)$ is species abundance at time t, and K is carrying capacity. In the absence of disturbances, the system is drawn towards the equilibrium $x(t) = 0$ (i.e. $N(t) = K$). Thus, $x(t)$ describes the difference between current and equilibrium abundance (Figure 1). In this model, dynamics follow

$$\frac{dx}{dt} = -rx + \xi(t) \quad (3)$$

where r is the intrinsic growth rate, and $\xi(t)$ is a stochastic function representing the effects of disturbance events. Disturbances occur at discrete moments in time, $t_i$, with time between disturbances drawn from an exponential distribution with frequency $\lambda$ (Gillespie, 1976, 1992). The effects of disturbances on species abundances are drawn from a normal distribution with mean $\mu = 0$, and standard deviation $\sigma$.

Our second model extends Equation 3 to consider dynamics of a single species in multiple patches. Here, dynamics follow

$$\frac{dx_k}{dt} = -r_k x_k + \xi_k(t) - D_k N_k + (a/A) D_l \Sigma_l N_l \quad (4)$$

where $x_k = N_k - K_k$ describes standardised abundance in patch k, $D_k$ is dispersal rate out of patch k, a is patch area, A is the area of the total focal region, and $r_k$ and $K_k$ are patch-level growth rate and carrying capacity. Disturbances occur simultaneously across patches, but disturbance impacts on each patch, $\xi_k(t)$, are drawn from a multivariate normal distribution with inter-patch covariance $cov(\xi_k, \xi_l)$. Thus, correlation in stochastic forcing describes, for example, multi-patch responses to a shared environmental disturbance, such as a flood or drought.

Finally, our third model considers multiple interacting species within a single patch. Here, for a community of M total species, dynamics for species i follow

$$\frac{dx_i}{dt} = -r_i x_i + \xi_i(t) + \sum_{j \neq i} \alpha_{ij} x_j \quad (5)$$

where $\alpha_{ij}$ describes the effect of species j on species i, and $\xi_i(t)$ describes disturbance impacts on species i. In this model, we linearise dynamics around the multi-species equilibrium achieved by each of the M coexisting species, such that $x_i = N_i - N_i^*$ (i.e. distance from the multi-species equilibrium abundance). Under these circumstances, parameters r, K, and $\alpha$ can be interpreted in the same way as in the classic Lotka-Volterra equations. Effects of disturbances in this model also occur simultaneously across species, but their strengths differ, and are drawn from a multivariate normal distribution with covariances $cov(\xi_i, \xi_j)$.
Following the definitions above, stability for the model with one species and patch in Equation 3 relates directly to model parameters, with *r* describing resilience, resistance inversely related to *σ*, and invariance inversely related to var(*x*). As we will discuss, stability in the multi-patch and multi-species models in Equations 4 and 5 varies across spatial and ecological scales, following Equations 2a–2c. Additionally, as discussed in the *Measuring stability* section below, measurement biases arise in all three models as a function of temporal scale, especially for slow sampling rates.

## Simulating dynamics

For each model, we simulated 1000 independent iterations of 100 time-units each, and, after allowing a 20 time-unit “burn-in” period, measured resilience, resistance, and invariance following the empirical methods described below. Simulations were conducted in R version 3.6.1 (R Development Core Team, 2019). Source code is archived on Zenodo (https://doi.org/10.5281/zenodo.4626668).

For the single species, single patch model, we varied temporal scale by “sampling” the simulated timeseries at frequencies ranging from every 0.01 to every 10 time-units, with growth rate *r* = 1, standard deviation of disturbance effects *σ* = √0.1, and average frequency of disturbances *λ* = 1. Note that because both this model and the multi-species model in Equation 5 can be expressed fully in terms of standardised abundances *x* or *x*<sub>*i*</sub>, dynamics can be simulated without directly including effects of carrying capacity *K*.

For spatial scale, we used Equation 4 to simulate a single species in 30 separate patches, and “sampled” regions that included anywhere from 1 to all 30 patches, with *r*<sub>*k*</sub> and *K*<sub>*k*</sub> drawn from normal distributions with mean 1 and standard deviation 0.1, and inter-patch disturbance covariance cov(ξ<sub>*k*</sub>, ξ<sub>*l*</sub>) = σ<sup>2</sup>/2 = 0.05 for all patches (i.e. positive correlations among disturbances, with &lt;ρ(ξ<sub>*k*</sub>)&gt; = 1/2). For the results presented in the main text, we set dispersal rate *D* = *D*<sub>*k*</sub> = 1 for all patches (i.e. global dispersal). In the supplement, we also present results for global dispersal rates ranging from *D* = 0 to *D* = 2, directional dispersal where *D*<sub>*k*</sub> differs among patches, and dispersal in an open system, where we included an additional loss rate, −*LN*<sub>*k*</sub>, on the right-hand side of Equation 4 to represent net loss of propagules dispersing out of the focal region.

For ecological scale, we used Equation 5 to simulate a single patch with a community of 30 species, and “sampled” groupings containing anywhere from a single species to all 30, with between-species disturbance covariance cov(ξ<sub>*i*</sub>, ξ<sub>*j*</sub>) = 0 for species pairs (i.e. no correlation in effects of disturbances). Growth rates *r*<sub>*i*</sub> were again drawn from a normal distribution with mean zero and standard deviation 0.1, and interspecific interaction strengths *α*<sub>*ij*</sub> were drawn from a censored normal distribution with mean −1/2, standard deviation 0.1, and upper limit zero (n.b. we

avoid positive interactions in our model because these tend to cause runaway growth). For both the spatial and ecological scale simulations, we set disturbance frequency *λ* = 1 and used a “fast” sampling interval of one sample per 0.1 time-units in order to avoid measurement biases.

## Measuring stability

If sampling rates are slow relative to system dynamics, empirical estimates of resilience and resistance can be biased. This bias occurs because repeated disturbances between observations drive the system away from equilibrium, therefore leading to underestimates of resilience (i.e. *r*), whereas lags between disturbance events and observations lead to underestimation of disturbance effects (i.e. *σ*). In contrast, because sample-size corrected variance is an unbiased estimator, different sampling frequencies do not bias measurements of invariance (i.e. var(*x*)<sup>−1</sup>).

Importantly, we can leverage information from unbiased estimates of variance to correct biased estimates of resilience and resistance. Following the properties of point processes (Arnoldi et al 2019; Zelnik et al., 2018), average temporal variance of *x* in the single species, single patch model in Equation 3 is

$$ \mathrm{var}(x) = \sigma^2 \lambda / (2r) \qquad (6) $$

Thus, variance decreases, and invariance increases, with higher resistance (lower *σ*), higher resilience (larger *r*), and slower disturbance frequency (lower *λ*). Building on this formula, we can derive an estimate for how the distance of *x* from equilibrium changes over time, yielding

$$ E[x(t + \tau)^2] = x(t)^2 \exp[-2r\tau] + \sigma^2 \sum_{i=1}^{p(t+\tau)} (\exp[-2ri/\Lambda]) \quad (7) $$

where *t* represents the time of an observation, *t* + *τ* is the time of the observation made immediately after *t*, *p*(*t* + *τ*) is the number of disturbances that occurred between *t* and *t* + *τ*, and 1/Λ is the average waiting time between disturbances. For our simulations, we applied this function to observations of individual species and patches to obtain unbiased estimates of *r* and *σ*. Because this expression is non-separable, we estimated parameters using a nonlinear optimiser (either nls in the stats package, or gnls in the nlme package in R; Pinheiro & Bates, 2000). More details, and automated functions for applying these corrections, are available in Appendix E and in the ecostatscale R package.

## RESULTS

Matching theoretical expectations for our single-species single-patch model, simulations showed no
temporal scaling relationships for invariance, nor for bias-corrected estimates of resilience and resistance based on Equation 7 (Figure 3). Raw estimates of both were biased (Figure 3a,b), and for temporal scales above ~3, raw resistance could not be computed, as no disturbance-free reference time steps were available. For all metrics, uncertainty increased at slower sampling rates.

For the spatial model with global dispersal $D_k = 1$ for all patches, both $\sigma_B$ and var($x_B$) increased (and thus, resistance and invariance decreased) with spatial scale,

<figure><img src="image_5.png" /><figcaption>**FIGURE 3** Effects of temporal scale on observed stability, following the working definitions in Figures 1 and 2. Recall that resilience is directly proportional to r, whereas resistance and invariance are inversely related to σ and var(x), respectively (signified with "-1" superscript). Shaded regions in (a-b) show ±one standard deviation from the mean for corrected parameter estimates, following Equation 7. Hatched regions show raw estimates of r and σ, calculated by comparing abundances in sequential time windows with differing numbers of disturbances—note resulting bias (see main text for details). Shaded region in (c) shows temporal variance calculated directly from the raw timeseries—note that this estimate is unbiased. Intervals are calculated based on 1000 model iterations simulated over 100 time-units, with sampling interval specified on the horizontal axis, and parameter values r = 1, D = 0, μ = 0, σ² = 0.1, and λ = 1. Dashed lines show true values for r and σ, and analytical expectation for var(x) based on Equation 6</figcaption></figure>

matching expected relationships from Equations 2a–2b (Figure 4b,c). Following the expected scaling relationship in Equation 2c, resilience declined slightly with spatial scale, eventually settling on the average inter-patch growth rate $r_k = 1$ (Figure 4a). There was no effect of changing global dispersal rates, nor of allowing variable dispersal rates among patches, largely because of compensatory effects at the patch vs. landscape levels (see Figure S1a,b, and more detailed discussion in Zelnik et al., 2018). In contrast, for the model with an additional

<figure><img src="image_6.png" /><figcaption>**FIGURE 4** Effects of spatial scale on observed stability, following the working definitions in Figures 1 and 2. All results are for a system with 30 patches, but measured at different spatial scales. (a) Resilience declines somewhat as a function of spatial scale, due to stronger buffering effects among plots at smaller scales. (b-c) Resistance and invariance differ across scales due to changes in total abundance and between-patch covariance. In all panels, dashed lines show analytical expectation for scaling relationships, following Equations 2a–2c. Sampling interval is 0.1 time-units for all simulations, with global dispersal ($D_k = D = 1$ for all patches) and disturbances with positive between-patch covariance (cov($\xi$) = $\sigma^2/2$). Patch-level growth rate $r_k$ and carrying capacity $K_k$ are drawn from normal distributions with mean 1 and standard deviation 0.1. See Figures S1a-c in the supplement for examples of other types of dispersal. Otherwise, intervals and parameters are as described in the legend for Figure 3</figcaption></figure>
loss rate from dispersal outside of the focal area, mean abundance declined as a function of loss rate, leading to large changes in stability estimates (Figure S1c).

<figure><img src="image_7.png" /><figcaption>**FIGURE 5** Effects of ecological scale on observed stability, following the working definitions in Figures 1 and 2. All results are for a system with 30 interacting species, but measured at different ecological scales. (a) Resilience varies with ecological scale in this model because of compensatory dynamics among species. (b) Resistance depends only on instantaneous responses to the disturbance regime, and therefore varies with ecological scale, but not with species interactions (n.b. dashed and dotted lines overlap perfectly). (c) Invariance varies quadratically with ecological scale as a function of average species-level variance in abundance, and between-species correlation—negative correlation implies a concave-down relationship (as shown here), whereas positive correlation implies a concave-up relationship (see Figure S2 in the appendix for other examples from empirical datasets). Note that the vertical axis is log-transformed for easier visualisation. For all panels, dashed lines show analytical expectation for scaling relationships, following Equations 2a–2c. For reference, dotted lines show hypothetical scaling relationships for a system without interactions. Sampling interval is 0.1 time-units for all simulations, with species interactions ($\alpha_{ij}$) drawn from a truncated normal distribution (mean = -1/2; standard deviation = 0.1; upper limit = 0), and species growth rates $r_i$ drawn from a normal distribution with mean 1 and standard deviation 0.1. Otherwise, intervals and parameters are as described in the legend for Figure 3</figcaption></figure>

Finally, for the multi-species model (Figure 5), both $r$ and $\sigma_B$ increased (and thus, resilience increased, and resistance decreased) with ecological scale following the expected scaling relationships in Equations 2b–2c. Invariance also matched the expected analytical relationship Equation 2a, resulting in a concave-down hump-shaped relationship for $var(x_B)$.

## DISCUSSION

Our results demonstrate that resilience, resistance, and invariance can be related to simple statistical properties of dynamical systems, which allow general and robust scaling of these stability metrics across time, space, and level of ecological organisation. Recall that following our definitions, resistance describes the immediate effect of disturbances on abundance, resilience describes the rate at which abundances recover from disturbances, and invariance describes the joint effects of these two processes on dynamics over time. The scaling functions in Equations 2a–2c demonstrate how to extrapolate all three of these metrics across spatial and ecological scales. Similarly, the relationships among these variables that we summarise in Equation 6, and the correction introduced in Equation 7, demonstrate how to overcome measurement biases that arise as a function of temporal scale.

Jointly, these methods allow empirical estimates of resilience, resistance, and invariance measured at one set of scales to be extrapolated out to other, unobserved scales. These results are therefore important both because they highlight the fact that stability measurements are scale dependent, and because they demonstrate how to correct for this scale-dependence when comparing and interpreting empirical results. The general implications and limitations of these methods are discussed below, and a guide for practitioners showing how to apply our methods using the ecostatscale R package is presented in Box 1.

### General scaling relationships

Because they arise from basic mathematical definitions of variance and covariance, there are several important cases for which the scaling relationships in Equations 2a–2c hold exactly, regardless of the underlying processes governing system dynamics. Most obviously, stability measurements from smaller scales can always be extrapolated to larger scales provided that both include the same total extent (e.g. using the full set of small plots in Figure 1d to estimate stability of the single large plot in Figure 1f). However, even if only partial sampling is available, the scaling relationships provide unbiased extrapolations to larger scales so long as smaller scale sampling is representative (i.e. the average conditions at the
# BOX 1 Recommendations for practitioners

## Harmonising data

Differences in observational scales greatly influence stability measurements, which presents a problem for synthesis and metanalysis because available data are often measured at many different scales. For example, consider the effects of ecological scales (Figure 5): because all three metrics of stability vary as a function of the number of species included in samples, studies with different sampling designs could reach wildly different conclusions about stability, even if they were all conducted in the same system. A potential solution is to “harmonise” records across studies, by extrapolating estimates to a unified set of scales. Doing so removes the statistical effects of scale differences, meaning that remaining differences should be ecologically meaningful. To improve interpretability, we suggest that harmonisation should generally proceed by extrapolating smaller-scale studies up to larger scales (see General scaling relationships). For example, if studies are available with plot sizes of 1, 5, and 10 m², then stability estimates for all three studies should be harmonised to a scale of 10 m².

Unbiased estimates of resilience, resistance, and invariance are independent of temporal scale, meaning that harmonisation is not necessary to account for differential sampling rates (however, a bias correction such as that in Equation 7, and implemented in the xt2fun function, may still be necessary if sampling is slow relative to system dynamics). In contrast, stability measurements do vary with spatial and ecological scale, meaning that harmonisation is required to account for differences in observational scale. Harmonisation is accomplished by taking available estimates of invariance, resistance, or resilience, plugging them into Equations 2a–2c, respectively, and substituting *b* with the scales of observations, and *B* with the larger scale to which observations are to be harmonised. These computations can be carried out using the var_scale, sd_scale, and res_scale functions, respectively. See the help documentation in the ecostatscale R package for more details.

## Extrapolating estimates

Equations 2a–2c and the corresponding R functions also allow extrapolation of expected stability beyond the scale of observations. These extrapolations can be useful for, for example, estimating stability at the landscape level based on plot-level observations, or at the community level given observations of a subset of species. As with harmonisation, extrapolation is accomplished by substituting *b* with the scale of observations, and *B* with the desired scale at which extrapolations are to be made.

Importantly, these extrapolations only describe “average” conditions across observations. Thus, in order for extrapolations to be accurate, observations must be representative of the extrapolated region or species—that is, average conditions across the samples must be the same as the average conditions at the extrapolated scale. In general, random samples of plots or species should achieve this criterion. However, if heterogeneity is structured across space or time, then this must be accounted for in the sampling design. For example, if a landscape is known to be 20% mesic and 80% xeric, then sampling should be stratified to include this ratio of site types, even if it differs from local conditions near survey plots.

## Detecting heterogeneity

Although unobserved heterogeneity can bias extrapolations of stability, data measured across several scales can be leveraged to identify underlying sources of heterogeneity and to validate extrapolations. To do so, we suggest parameterising the statistical scaling relationships in Equations 2a–2c, or the corresponding R functions, with data from the smallest observational scale available, and comparing these to independent estimates of stability observed at larger scales (Pimm et al., 2019; Wang et al., 2017). If empirically observed relationships diverge from statistical expectations, then the associated scaling relationship is likely to be biased, and should not be used to extrapolate stability to other, unobserved scales.

Even if observations are only available from a single scale, it may still be possible to identify heterogeneity by aggregating specific combinations of plots or species through “non-random cross-validation” (Wenger & Olden, 2012). For example, if stability varies among plots in different regions of a site (e.g. northwest vs. southeast corner of a field), or among different species groups (e.g. legumes vs. grasses), then these deviations indicate spatial or ecological heterogeneity, respectively. Similarly, if stability varies between estimates derived from the full dataset vs. for example, a subset of the dataset including every second sampling event, then this indicates that temporal sampling was potentially too infrequent to produce unbiased estimates of resilience and resistance, and a correction such as that in Equation 7, or the accompanying xt2fun R function, should be applied. See function help documentation for example implementations.
smaller scale match those at the larger scale). Thus, the scaling definitions should be applicable for most studies that randomly sample species from a larger community, or plots from a larger landscape—which includes a substantial fraction of existing ecological timeseries datasets. For similar reasons, our methods can also be used to project stability measurements from larger scales to smaller scales, although these estimates summarise *average conditions* at smaller scales—not conditions within individual patches or for individual species.

Another general result is that if sampling is slow relative to system dynamics, raw estimates of resilience and resistance will be biased (leading to under-estimates of *r* and *σ*). In theory, the bias correction that we present in Equation 7 holds only for simple linear systems—however, in practice, it will often hold for more complex systems so long as they are undergoing bounded fluctuations around a static mean (Arnoldi et al 2019; Zelnik et al., 2018). For example, estimates from Equation 7 converged to the correct parameter values for both our multi-patch and multi-species models.

Our study also highlights the similar effects of spatial vs. ecological scale on stability. For invariance, as predicted in Equation 2a, positive average covariance among species-level abundances or site-level abundances leads to a concave-up scaling relationship for var(x) (e.g. Figure 4c, as might be driven by species in nearby plots all responding similarly to a local disturbance). Negative average covariance leads to a concave-down relationship (e.g. Figure 5c, as might be driven by strong effects of competition). Recall, however, that for the simulation results presented here, we assume that species responses to disturbances are independent (i.e. cov(ξᵢ, ξⱼ) = 0), which allows competition effects to dominate. In real-world systems, these responses are generally thought to covary positively over time (e.g. due to similar effects of droughts or floods across species), leading to net positive average covariance in species abundances, and thus a concave-up relationship for var(x) as a function of ecological scale (Houlahan et al., 2007; Loreau & de Mazancourt, 2013). Examples of concave-up, concave-down, and linear scaling relationships derived from empirical data are shown in Figure S2 in the supplement.

The spatial and ecological scaling relationships for resistance in Equation 2b are similar to those for invariance. However, because we define resistance as the instantaneous effect of disturbances on abundances, it is largely independent of internal system processes, such as inter-patch dispersal or species interactions. Scaling relationships are therefore the same for systems with and without species interactions (Figure 5b; n.b. dashed and dotted lines overlap perfectly). This property means that extrapolating resistance across spatial and ecological scales should be relatively straightforward even in complex systems.

The scaling relationship for resilience is the most variable. For temporal scaling, average per-capita growth

rates in our model are constant across space and time, leading to no change in resilience with scale (Figures 3a and 4a). For multi-patch and multi-species communities, the median scaling relationship for resilience can be approximated following Equation 2c, although for multi-species communities, this relationship only describes the instantaneous recovery rate (i.e. as would be observed immediately after a disturbance). This is because resilience also varies as a function of temporal scale in these systems (see Figure S3)—that is, different return rates occur over short vs. long timespans, because interactions among patches or species can lead to complex, non-monotonic dynamics (see Neubert & Caswell, 1997). In theory, this temporal scaling relationship can be expressed analytically (Eq. 6b in Arnoldi et al., 2018), but in practice, doing so requires substantial *a priori* information, suggesting that empiricists should generally be careful to choose sampling intervals that match the scale of the phenomena being tested when studying highly interconnected or diverse systems.

## Interpreting deviations from scaling rules

An important limitation of the statistical scaling relationships in Equations 2a–2c is that they are largely phenomenological. Although they predict *how* patterns change across scales, they do not explain *why* these patterns exist. Indeed, there are usually many potential mechanisms that could explain patterns equally well. For example, although a concave-up relationship between var(x) and ecological scale implies negative average correlations in interspecific abundances, these correlations could result from many different mechanisms, including competition, species differences in their response to disturbances, or observation error (Clark et al., 2019; Lee et al., 2020).

Critically, cases where observed scaling relationships diverge from statistical expectations can be especially informative. Because Equations 2a–2c hold only if samples are representative of the full extent under consideration, divergence from these expectations implies that there is unobserved heterogeneity in observations across time, space, or species. For example, if scaling relationships predicted from patch-level or species-level data do not match larger-scale observations (e.g. from remote sensing data, or at the community level), then this might indicate that samples are not representative of conditions at the larger-scale (König et al., 2017; Wang et al., 2017). Similarly, if observed invariance changes as a function of temporal scale, then this divergence may be indicative of regime shifts (Scheffer et al., 2009), or that equilibrium abundances are not static (Pimm et al., 2019).

Note, however, that without additional information, the statistical scaling relationships cannot predict how unobserved heterogeneity outside of the sampled region will influence predictions. In diverse communities, for
example, stability varies with community structure, meaning that adding or removing species, or changes in how species interact with one another, will alter the scaling relationships (Arnoldi et al., 2018; Hillebrand et al., 2018; Jacquet et al., 2016). Similarly, if the structure of heterogeneity changes across space or time (e.g. clustered disturbances in one region, vs. disaggregated disturbances in another), then resulting biases will be, for all practical purposes, impossible to correct without a priori knowledge of conditions across the entire landscape (Arnoldi et al., 2019; Hallett et al., 2019; Lee et al., 2020).

## Nonstationary systems

In this study, we largely neglect effects of non-stationary behaviour on stability (e.g. dynamics with no fixed equilibrium point). We do so because these kinds of dynamics are generally difficult to summarise through simple statistical relationships. However, it is possible to incorporate some kinds of non-stationary behaviour into our methods. If mean abundances change slowly over time, then our scaling functions can still be applied within subsets of the data. For example, in systems undergoing regime shifts, our scaling functions could be fitted separately for the time periods before vs. after the shift (Carpenter et al., 2001; Scheffer et al., 2009). Similarly, even for systems undergoing rapid changes due to anthropogenic pressures, our methods likely still apply for historical data preceding current disturbance regimes (Coulson, 2021). More broadly, if trends in average abundances are predictable, then we can redefine system dynamics around this moving average. We demonstrate this property in Figure 5c, where losses due to dispersal out of the system lead to declines in average abundance below the expected equilibrium at $N_k(t) = K_k$. Without any additional corrections, this shift in equilibrium causes stability scaling relationships to diverge greatly from statistical expectations. However, if we redefine *x* around the new mean abundance value (i.e. $x_{k,\text{new}}(t) = N_k(t) - \langle N_k \rangle$) then we find that the statistical scaling relationships still hold. Using more sophisticated approaches for predicting changes in mean abundance values over time (e.g. Cenci & Saavedra, 2019; Chesson, 2017; Deyle et al., 2016; Hamilton et al., 2017; Karakoç et al., 2020), it may be possible to similarly “detrend” even very complex nonstationary data to make them suitable for our methods.

## Conclusions

The strong dependence of ecological stability on scale and context has represented a major challenge for attempts to synthesise and extrapolate measurements from

real-world systems. The methods that we present here greatly alleviate these challenges, and provide an opportunity to harmonise measurements from across different sites and systems, and to extrapolate measurements out to scales which are difficult to measure. We are therefore hopeful that these results will be especially helpful in facilitating cross-system comparisons of stability, providing estimates at scales that are more relevant for conservation and management, and identifying underlying drivers of heterogeneity across time, space, and species.

## ACKNOWLEDGEMENTS

This paper was conceived as part of the “Filling in gaps in global understanding of ecological stability and coexistence” (FIGS) workshop, funded by a UFZ Program Synthesis Grant, UFZ IP-11 Project Integration Funds, and the “TULIP” French Laboratory of Excellence (ANR-10-LABX-41; ANR-11-IDEX-0002-02). The contributions of ATC, SH, JMC, and AM, were also supported, in part, by the German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig funded by the German Research Foundation (FZT 118). ML, CdM, and YRZ were supported by the BIOSTASES Advanced Grant funded by the European Research Council under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 666971). GB was funded by the Swedish Research Council (grant VR 2017-05245); LGS was funded by an NSF EPSCoR Track 1 RII grant (NSF award #EPS-1655726). The authors also thank K. Thompson, the Harpole and Chase lab groups at iDiv and UFZ, and the Hillebrand lab group at HIFMB, for helpful feedback on earlier drafts of this manuscript. Open Access funding enabled and organized by Projekt DEAL.

## AUTHORSHIP

All authors contributed to the conception and planning of the paper as part of the FIGS working group at UFZ, and/or contributed significantly to revisions of the manuscript; ATC and YRZ planned and led the FIGS workshop; ATC wrote the first draft of the paper and conducted the analyses; ATC, JFA, and YRZ developed the mathematical framework and scaling relationships.

## PEER REVIEW

The peer review history for this article is available at https://publons.com/publon/10.1111/ele.13760.

## DATA AVAILABILITY STATEMENT

All data used in this study are already available through published, publicly available databases, cited by DOI in this manuscript. Full source code and data used in analyses are available via Zenodo (https://doi.org/10.5281/zenodo.4626668), and the ecostatscale R package is available via CRAN, GitHub (github.com/adamtclark/ecostatscale), and Zenodo (https://doi.org/10.5281/zenodo.4626672).
ORCID

Adam Thomas Clark https://orcid.

org/0000-0002-8843-3278

Jean-François Arnoldi https://orcid.

org/0000-0001-9373-4051

Viktoriia Radchuk https://orcid.

org/0000-0003-3072-0095

Ian Donohue https://orcid.org/0000-0002-4698-6448

Claire Jacquet https://orcid.org/0000-0002-7913-0689

Lauren G. Shoemaker https://orcid.

org/0000-0002-4465-8432

Shaopeng Wang https://orcid.

org/0000-0002-9430-8879

Jonathan M. Chase https://orcid.

org/0000-0001-5580-4303

REFERENCES

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., Loreau, M. & Haegeman, B. (2019) The inherent multidimensionality of temporal variability: How common and rare species shape stability patterns. *Ecology Letters*, 22, 1557–1567.

Black, A.J. & McKane, A.J. (2012) Stochastic formulation of ecological models and their applications. *Trends in Ecology & Evolution*, 27, 337–345.

Carpenter, S., Walker, B., Anderies, J.M. & Abel, N. (2001) From metaphor to measurement: Resilience of what to what? *Ecosystems*, 4, 765–781.

Cenci, S. & Saavedra, S. (2019) Non-parametric estimation of the structural stability of non-equilibrium community dynamics. *Nature Ecology & Evolution*, 3, 912–918.

Chesson, P. (2000) Mechanisms of maintenance of species diversity. *Annual Review of Ecology and Systematics*, 31, 343–366.

Chesson, P. (2017) AEDT: A new concept for ecological dynamics in the ever-changing world. *PLOS Biology*, 15, e2002634.

Clark, A.T., Hillebrand, H. & Harpole, W.S. (2019) Scale both confounds and informs characterization of species coexistence in empirical systems. *The American Naturalist*, 194, 794–806.

Coulson, T. (2021) We live in a changing world, but that shouldn't mean we abandon the concept of equilibrium. *Ecology Letters*, 24(1), 3–5.

Csillag, F., Fortin, M.J. & Dungan, J.L. (2000) On the limits and extensions of the definition of scale. *Bulletin of the Ecological Society of America*, 81(3), 230–232.

Deyle, E.R., May, R.M., Munch, S.B. & Sugihara, G. (2016) Tracking and forecasting ecosystem interactions in real time. *Proceedings of the Royal Society B: Biological Sciences*, 283, 20152258.

Domínguez-García, V., Dakos, V. & Kéfi, S. (2019) Unveiling dimensions of stability in complex ecological networks. *Proceedings of the National Academy of Sciences of the United States of America*, 116, 25714–25720.

Donohue, I., Hillebrand, H., Montoya, J.M., Petchey, O.L., Pimm, S.L., Fowler, M.S. et al. (2016) Navigating the complexity of ecological stability. *Ecology Letters*, 19, 1172–1185.

Donohue, I., Petchey, O.L., Montoya, J.M., Jackson, A.L., McNally, L., Viana, M. et al. (2013) On the dimensionality of ecological stability. *Ecology Letters*, 16, 421–429.

Gillespie, D.T. (1976) A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. *Journal of Computational Physics*, 22, 403–434.

Gillespie, D.T. (1992) *Markov processes: An introduction for physical scientists*. London, UK: Academic Press.

Grimm, V. & Wissel, C. (1997) Babel, or the ecological stability discussions: An inventory and analysis of terminology and a guide for avoiding confusion. *Oecologia*, 109, 323–334.

Haegeman, B., Arnoldi, J.-F., Wang, S., de Mazancourt, C., Montoya, J.M. & Loreau, M. (2016) Resilience, invariability, and ecological stability across levels of organization. *bioRxiv*, 085852.

Hallett, L.M., Shoemaker, L.G., White, C.T. & Suding, K.N. (2019) Rainfall variability maintains grass-forb species coexistence. *Ecology Letters*, 22, 1658–1667.

Hamilton, F., Berry, T. & Sauer, T. (2017) Kalman-Takens filtering in the presence of dynamical noise. *The European Physical Journal Special Topics*, 226, 3239–3250.

Hillebrand, H., Langenheder, S., Lebret, K., Lindström, E., Östman, Ö. & Striebel, M. (2018) Decomposing multiple dimensions of stability in global change experiments. *Ecology Letters*, 21, 21–30.

Holling, C.S. (1973) Resilience and stability of ecological systems. *Annual Review of Ecology and Systematics*, 4, 1–23.

Houlahan, J.E., Currie, D.J., Cottenie, K., Cumming, G.S., Ernest, S.K.M., Findlay, C.S. et al. (2007) Compensatory dynamics are rare in natural ecological communities. *Proceedings of the National Academy of Sciences of the United States of America*, 104, 3273–3277.

Ibell, F., Cowles, J., Dee, L.E., Loreau, M., Reich, P.B., Gonzalez, A. et al. (2018) Quantifying effects of biodiversity on ecosystem functioning across times and places. *Ecology Letters*, 21(6), 763–778. https://doi.org/10.1111/ele.12928

Ives, A.R. & Carpenter, S.R. (2007) Stability and diversity of ecosystems. *Science*, 317, 58–62.

Jacquet, C., Moritz, C., Morissette, L., Legagneux, P., Massol, F., Archambault, P. et al. (2016) No complexity-stability relationship in empirical ecosystems. *Nature Communications*, 7, 12573.

Karakoç, C., Clark, A.T. & Chatzinotas, A. (2020) Diversity and coexistence are influenced by time-dependent species interactions in a predator-prey system: Coexistence in a predator-prey system. *Ecology Letters*, 23(6), 983–993. https://doi.org/10.1111/ele.13500

Kéfi, S., Domínguez-García, V., Donohue, I., Fontaine, C., Thébault, E. & Dakos, V. (2019) Advancing our understanding of ecological stability. *Ecology Letters*, 22, 1349–1356.

König, S., Worrich, A., Centler, F., Wick, L.Y., Miltner, A., Kästner, M. et al. (2017) Modelling functional resilience of microbial ecosystems: Analysis of governing processes. *Environmental Modelling & Software*, 89, 31–39.

Lande, R., Engen, S. & Sæther, B. (1999) Spatial scale of population synchrony: Environmental correlation versus dispersal and density regulation. *The American Naturalist*, 154, 271–281.

Lee, A.M., Sæther, B. & Engen, S. (2020) Spatial covariation of competing species in a fluctuating environment. *Ecology*, 101, e02901.

Leibold, M.A. & Chase, J.M. (2018) *Metacommunity ecology* (Vol. 59). Princeton, NJ: Princeton University Press.

Levin, S.A. (1992) The problem of pattern and scale in ecology: The Robert H. MacArthur Award Lecture. *Ecology*, 73, 1943.

Lewontin, R.C. (1969) The meaning of stability. *Brookhaven Symposium in Biology*, 22, 13–24.

Loreau, M. & de Mazancourt, C. (2013) Biodiversity and ecosystem stability: A synthesis of underlying mechanisms. *Ecology Letters*, 16, 106–115.

May, R.M. (1973) *Stability and complexity in model ecosystems*. Monographs in population biology. Princeton, NJ: Princeton University Press.

Neubert, M.G. & Caswell, H. (1997) Alternatives to resilience for measuring the responses of ecological systems to perturbations. *Ecology*, 78, 653–665.

Pimm, S.L. (1984) The complexity and stability of ecosystems. *Nature*, 307, 321–326.

Pimm, S.L., Donohue, I., Montoya, J.M. & Loreau, M. (2019) Measuring resilience is essential to understand it. *Nature Sustainability*, 2, 895–897.
Pinheiro, J.C. & Bates, D.M. (2000) *Mixed-effects models in S and S-PLUS*. New York, NY: Springer.

R. Development Core Team. (2019) *R: A language and environment for statistical computing*. Vienna, Austria: R Foundation for Statistical Computing.

Scheffer, M., Bascompte, J., Brock, W.A., Brovkin, V., Carpenter, S.R., Dakos, V. et al. (2009) Early-warning signals for critical transitions. *Nature*, *461*, 53–59.

Stommel, H. (1963) Varieties of oceanographic experience: The ocean can be investigated as a hydrodynamical phenomenon as well as explored geographically. *Science*, *139*, 572–576.

Wang, S., Lamy, T., Hallett, L.M. & Loreau, M. (2019) Stability and synchrony across ecological hierarchies in heterogeneous meta-communities: Linking theory to data. *Ecography*, *42*(6), 1200–1211. https://doi.org/10.1111/ecog.04290

Wang, S., Loreau, M., Arnoldi, J.-F., Fang, J., Rahman, K.A., Tao, S. et al. (2017) An invariability-area relationship sheds new light on the spatial scaling of ecological stability. *Nature Communications*, *8*, 15211.

Wenger, S.J. & Olden, J.D. (2012) Assessing transferability of ecological models: An underappreciated aspect of statistical validation: Model transferability. *Methods in Ecology and Evolution*, *3*, 260–267.

Wilkinson, D.J. (2011) *Stochastic modelling for systems biology*. Mathematical and computational biology, 2nd edition. London, UK: Chapman & Hall.

Wilkinson, D.J. (2018) *smfsb: Stochastic modelling for systems biology*. R package version 1.3.

Zelnik, Y.R., Arnoldi, J.-F. & Loreau, M. (2018) The impact of spatial and temporal dimensions of disturbances on ecosystem stability. *Frontiers in Ecology and Evolution*, *6*, 224.

## SUPPORTING INFORMATION

Additional supporting information may be found online in the Supporting Information section.

> How to cite this article: Clark AT, Arnoldi J-F, Zelnik YR, et al. General statistical scaling laws for stability in ecological systems. *Ecology Letters*. 2021;24:1474–1486. https://doi.org/10.1111/ele.13760