Spatial and Ecological Scaling of Stability in Spatial Community Networks

Introduction

Ecological stability is a property that can be broadly defined as the ability of an ecosystem to remain unaltered when challenged by perturbations. However, there exist multiple ways of characterizing stability, which leads to different stability definitions or components (Pimm, 1984; Grimm and Wissel, 1997; McCann, 2000). Different components include resistance (to perturbation), initial resilience (i.e., the ability to recover from a perturbation) or invariability (i.e., the ability to remain unaltered to repeated perturbations) (Figure 1A). Different stability components can also vary with scale, impeding cross-system comparison of stability, or be scale-independent instead (Levin, 1992; Wang et al., 2017; Domínguez-García et al., 2019; Kéfi et al., 2019; Greig et al., 2022; Figure 1B). Thus, the ecological and spatial scale at which one studies an ecological system can be hypothesized to influence stability assessments.

FIGURE 1

Figure 1. (A) Different stability components considered throughout the manuscript. In solid lines, we represent the hypothetical dynamics of the biomass of one species affected by multiple periodic disturbances. Resistance (blue lines) is defined as the inverse of the change of biomass as a consequence of the direct effect of the disturbance. The growth rate and the initial resilience of the system (green lines) both quantify the rate at which the system tends to recover to the equilibrium after the effect of the disturbance. The main difference between growth rate and initial resilience is whether it is computed relatively to the biomass prior the disturbance (growth rate), or relatively to the distance between the biomass after the perturbation and the equilibrium biomass (initial resilience). In both cases, we will consider the short-term response, in contrast to the asymptotic resilience (not considered in this manuscript). Finally, invariability (red dashed lines) represents the ability of the systems to remain close to the equilibrium when multiple disturbances affect (periodically or randomly) the species. Panel inspired by Clark et al. (2021). (B) In meta-communities, multiple species (“sp”) form communities located at multiple locations (“loc”) forming a region. We define the regional stability as the stability of the total biomass of one species across locations, and community stability as the stability of the total biomass across species in one location.

Previous studies have found that species diversity increases the invariability of communities (Thébault and Loreau, 2005; Tilman et al., 2006; Gross et al., 2014), usually as a consequence of the asynchrony of the population dynamics (Yachi and Loreau, 1999; Ives et al., 2000): variability decreases with species diversity because of the statistical averaging of the fluctuations in species’ abundances for not perfectly synchronous dynamics (Doak et al., 1998). Hence, communities are more invariable than their constituent species, and the ratio of the community and species invariabilities has been proposed as an estimate of the across-species synchrony of the population dynamics (Loreau and De Mazancourt, 2008). This result can create differences between species- and community-scale stability assessments (Flöder and Hillebrand, 2012; Mougi and Kondoh, 2012; Downing et al., 2014).

Analogously, studies of meta-communities (defined as sets of local communities that are linked by dispersal of multiple potentially interacting species; Leibold et al., 2004) have also proven that temporal invariability increases from local species to meta-communities (Wang and Loreau, 2014), again as a consequence of decreasing synchrony (Wang et al., 2019). The ratio between regional and local invariabilities could be employed as a proxy for a region-wide synchrony, which would represent the global degree of synchrony in the spatial network (not to be confused with the regional synchrony that might occur at smaller spatial scales: Moran, 1953; Lande et al., 1999; Jarillo et al., 2018, 2020). Heterogeneous environmental conditions (Chesson, 2000) and the dispersal ability of the species (Amarasekare, 2008) might further cause the spatial scale to influence population and community dynamics, and therefore spatial scale-dependence of stability. Similarly, also the precise spatial organization of the network may influence meta-community stability, as has been found when comparing riverine vs. linear networks (Fagan, 2002; Carrara et al., 2012; Altermatt, 2013; Liu et al., 2013; Peterson et al., 2013).

Finally, ecological and spatial scales may interact. For instance, spatial scale affects stability more in communities than in populations (Mougi and Kondoh, 2016), depending on the position of the focal species within the community (Limberger et al., 2019). Furthermore, decreasing the size of spatial networks reduces species richness more than one would expect from spatial samplings (Chase et al., 2020).

To understand the scaling of stability in meta-communities, which will allow the comparison of stability of systems analyzed at different scales, Clark et al. (2021) provided general scaling laws of common stability components: resistance, initial resilience, and invariance. They found that—if these measures are not normalized by biomass—invariance, resistance, and initial resilience decrease with the spatial scale (the size of the considered spatial region). While they found that ecological scale (the number of species) also in most cases reduced invariance and resistance, it increased initial resilience. Because total biomass often changes with scale, we first wanted to revisit the scaling laws established by Clark et al. (2021) by considering normalized stability measures.

Our objective is to investigate whether normalized regional and community stability measures are simple summary statistics of local and species stability, respectively, and do not change with the spatial or ecological scale at which they are studied (Figure 1B). Additive quantities (also known as extensive) are quantities that are simply added when combining several subsystems (IUPAC, 2019). They include the total biomass, its derivative, or its change due to a perturbation. The quotient of two additive quantities gives what is call an intensive quantity, which are scale-free quantities (independent of the number of subsystems) (IUPAC, 2019). Therefore, additive quantities automatically give scale-free quantities when normalized by another additive quantity like the number of subsystems or the total biomass.

We begin our paper (section “Biomass Normalized Stability Measures”) by introducing biomass-normalized stability measures of growth rate, resistance, initial resilience, and invariability. In section “Spatial and Ecological Scaling of Stability Measures,” we then show that resistance, growth rate, and initial resilience are all independent of scale. This is because the introduced normalized definitions of these stability components are intensive magnitudes, which do not scale with the size of the systems. In contrast, invariability is shown to be independent of scale only in the case of perfectly synchronous dynamics of species and local nodes (the different locations forming the spatial network). In the more realistic scenarios with at least some level of asynchrony, invariability is not an intensive quantity, and it increases with the network size. We then discuss how to use these results to statistically estimate regional/community stability from partial information (values on some local nodes or species). Section “Model Simulations” shows that our formulas compare well with model simulations. Finally, Section “Discussion” discusses various ecologically relevant aspects of the results: the influence of low-stability nodes or species on network stability, the relevance of the mathematical definitions, the implications for the empirical measurement of stability, and the implications for the stability-complexity debate.

Biomass Normalized Stability Measures

We introduce biomass-normalized measures for the main stability components: resistance (to a perturbation), growth rate (after a perturbation), initial resilience, and invariance (Figure 1A). These biomass-normalized measures avoid potential scale effects due to the usual scale-dependence of biomass. Biomass, its derivative and the change of biomass due to a perturbation are additive (extensive) quantities. Their quotients are expected to be scale-free measures (intensive magnitudes). We advance that for invariance the study will be more complicated as the variance is not additive.

Growth Rate

We define the growth rate R as the relative instantaneous return rate to the equilibrium of the biomass N after any sudden biomass change caused by any external perturbation at time t₀,

$R\equiv {\frac{1}{N\left(t\right)}\frac{\text{d}N\left(t\right)}{\text{d}t}|}_{t=t_0}.(1)$

Growth rate after a perturbation could be argued not to be a proper stability measure. Growth rate provides the rate of change of the population relative to the remaining population. Instead, initial resilience (section “Initial Resilience”) provides the rate of change of the population relative to the departure of the population to its equilibrium value. This makes the initial resilience a more intuitive stability measure, as it estimates the initial rate of return to equilibrium. However, as we will show below, the growth rate is directly proportional to the scale-dependence of initial resilience, which shows that it contains information on stability and therefore can be considered a component of stability.

Resistance

We define the resistance Ω as the inverse of the relative change of biomass as a consequence of a perturbation (Isbell et al., 2015; Baert et al., 2016),

$\mathrm{\Omega }\equiv \frac{N\left(t_0\right)}{N\left(t_0\right)-N\left(t_0+\mathrm{\delta }t\right)}.(2)$

Instead of working with this measure, sometimes its inverse Ω⁻¹ is referred to as resistance (Yang et al., 2019), with the possible conceptual disadvantage of presenting smaller values for more resistant systems. Other studies define resistance as the logarithm of the ratio of biomasses before and after any disturbance, ln⁡(N(t₀ + δt)/N(t₀)) (Hillebrand et al., 2018), whose absolute value will also decrease as systems become more resistant. Actually, in absolute value this logarithmic definition is at first order equivalent to Ω⁻¹ for small perturbations (as can be proven by applying a Taylor expansion on |N(t₀ + δt)−N(t₀)|/N(t₀)≪1), so its inverse is at first order equivalent to Eq. (2).

Initial Resilience

We define the initial resilience as the initial rate at which a biomass perturbation disappears, normalized by the extent of the perturbation

$\mathrm{\rho }\equiv {\frac{1}{\left|N\left(t\right)-N^{*}\right|}\frac{\text{d}\left(N\left(t\right)-N^{*}\right)}{\text{d}t}|}_{t=t_0},(3)$

where N* stands for the equilibrium biomass, which is assumed to be equal to the biomass just before the perturbation (N(t₀)). This definition stands for the short term recovery rate after a perturbation (Arnoldi et al., 2018), and has been sometimes referred to as reactivity (Neubert and Caswell, 1997). This initial resilience ρ can be expressed in terms of the already considered growth rate R (Eq. 1) and resistance Ω (Eq. 2),

$\begin{array}{c}\mathrm{\rho }\equiv {\frac{1}{\left|N\left(t\right)-N^{*}\right|}\frac{\text{d}\left(N\left(t\right)-N^{*}\right)}{\text{d}t}|}_{t=t_0}\hfill \\ ={\frac{N^{*}}{\left|N\left(t\right)-N^{*}\right|}\frac{1}{N^{*}}\frac{\text{d}N\left(t\right)}{\text{d}t}|}_{t=t_0}=\mathrm{\Omega }R.\hfill \end{array}(4)$

Invariability

We define invariability I as the ratio of the square temporal mean of the biomass and its temporal variance (Thibaut and Connolly, 2013):

$I\equiv \frac{{\left[{\text{mean}}_t\left(N\left(t\right)\right)\right]}^2}{{\text{var}}_t\left(N\left(t\right)\right)}.(5)$

This quantity is the inverse of the squared coefficient of variation of the biomass. Note that if the system is stable enough to stay away from extinction, this invariability will necessary be greater than 1; otherwise, environmental fluctuations might bring the biomass to zero. Other invariability estimates further normalize this invariability by the amplitude of environmental stochasticity (Haegeman et al., 2016; Arnoldi et al., 2019), in order to compare the invariability of systems subject to different environmental variability conditions.

Spatial and Ecological Scaling of Stability Measures

In this section, we address the spatial (local vs. regional) and ecological (species vs. community) scaling behavior (Figure 1B) of the previously introduced stability measures (section “Biomass Normalized Stability Measures,” Figure 1A). Spatial scaling refers to how the measure changes from the local level (e.g., one location) to the regional level (e.g., all locations). Ecological scaling refers to how the measure changes from the species level to the community level (e.g., all species). Knowing the response of stability metrics to scaling is important to build estimators that can be applied to the empirical study of extended ecological networks, which can only be partially sampled. We start with the study of the growth rate, as the simpler case, and follow with resistance, initial resilience, and invariability.

Growth Rate

Using our definition of growth rate (Eq. 1), we can compute the growth rate after a perturbation of a given species i at a given specific location x, R_x,i, as the normalized time derivative of the species local biomass, N_x,i,

$R_{x,i}={\frac{1}{N_{x,i}\left(t\right)}\frac{\text{d}{N}_{x,i}\left(t\right)}{\text{d}t}|}_{t=t_0}.(6)$

Defining the regional biomass of one species as the sum of all local biomasses of that species across the spatial network, N_i ≡ ∑_x N_x,i, and based on the mathematical definition of growth rate (Eq. 1) and on the sum rule of the derivative, we obtain that the regional growth rate R of the species i is

$R_i^{\left(ℛ\right)}\equiv {\frac{1}{N_i\left(t\right)}\frac{\text{d}{N}_i\left(t\right)}{\text{d}t}|}_{t=t_0}=\frac{{\sum }_x{N}_{x,i}\left(t_0\right)R_{x,i}}{{\sum }_x{N}_{x,i}\left(t_0\right)},(7)$

meaning that the regional growth rate of the species is the weighted arithmetic mean of local species growth rates, R_x,i, with weights equal to the local species biomasses at the moment of the perturbation, N_x,i(t₀) (Figure 2). Analogously, the local community growth rate $R_x^{\left(𝒞\right)}$, or the growth rate of the sum of biomasses across all the species of the community at a specific location N_x ≡ ∑_i N_x,i, is

$R_x^{\left(𝒞\right)}\equiv {\frac{1}{N_x\left(t\right)}\frac{\text{d}{N}_x\left(t\right)}{\text{d}t}|}_{t=t_0}=\frac{{\sum }_i{N}_{x,i}\left(t_0\right)R_{x,i}}{{\sum }_i{N}_{x,i}\left(t_0\right)},(8)$

i.e., the local community growth rate is the weighed arithmetic mean of local species growth rates of each of the species, R_x,i, with weights equals to the local species biomasses at the moment of the perturbation, N_x,i(t₀) (Figure 2). Finally, we can define regional community growth rate R^(ℛ𝒞) (equal to the regional growth rate of the community, or to the community growth rate of the spatial network), as the growth rate of the total biomass across species and locations, N_T ≡ ∑_x ∑_i N_x,i, which is given by

$\begin{array}{c}{R}^{\left(ℛ𝒞\right)}\equiv {\frac{1}{N_T\left(t\right)}\frac{\text{d}{N}_T\left(t\right)}{\text{d}t}|}_{t=t_0}=\frac{{\sum }_x{\sum }_i{N}_{x,i}\left(t_0\right)R_{x,i}}{{\sum }_x{\sum }_i{N}_{x,i}\left(t_0\right)}\hfill \\ =\frac{{\sum }_x{N}_x\left(t_0\right)R_x^{\left(𝒞\right)}}{{\sum }_x{N}_x\left(t_0\right)}=\frac{{\sum }_i{N}_i\left(t_0\right)R_i^{\left(ℛ\right)}}{{\sum }_i{N}_i\left(t_0\right)}.\hfill \end{array}(9)$

FIGURE 2

Figure 2. Regional (A) and community (B) growth rates, compared to local and species growth rates and to their unweighted arithmetic mean (AM ) (yellow circles), and biomass weighted arithmetic mean (AM_w) (green circles). AM and AM_w values closer to the identity line (black, dashed) estimate more precisely the regional (A) and community growth rate (B). Black dots are the growth rates of individual localities (A) or species (B) on which the means AM and AM_w were computed. Data were generated for 6,300 random communities of 10 competitors in 10-node random spatial networks (Supplementary Figure 1A and Supplementary Text), after a biomass decrease from the equilibrium affecting all species at all locations. AM_w are found to be better estimators of the growth rate, as expected from the results shown in the text.

Thus, the regional growth rate of a community after a perturbation, R^(ℛ𝒞), is the arithmetic mean of local community growth rates, weighted by the local total biomass; or equivalently the community growth rate of a spatial network is the arithmetic mean of the regional species growth rates, weighted by the species regional biomass.

Generally speaking, the network growth rate R_net of any ecological or spatial network is then given by the biomass-weighted arithmetic mean of the growth rates at the nodes (either representing the species for an ecological network, or the local nodes for a spatial network). R_net can represent either the regional value in a spatial network, the community value in an ecological network, or even the regional community value; computed using Eqs. (7), (8), or (9), respectively. This R_net can be also expressed as

$R_{net}={\mathrm{\mu }}_R\left(1+\frac{{\mathrm{\sigma }}_N}{{\mathrm{\mu }}_N}\frac{{\mathrm{\sigma }}_R}{{\mathrm{\mu }}_R}{c}_{N,R}\right),(10)$

where μ_N ≡ mean(N) and μ_R ≡ mean(R) denote unweighted population arithmetic means of biomasses and growth rates, respectively, ${\mathrm{\sigma }}_N\equiv \sqrt{\text{var}\left(N\right)}$ and ${\mathrm{\sigma }}_R\equiv \sqrt{\text{var}\left(R\right)}$ are their standard deviations computed as the square root of the variances, and c_N,R the normalized correlation between biomass and growth rate (see Table 1 for all the mathematical definitions). Given that −1 ≤ c_N,R ≤ 1, the network growth rate R_net can be greater or smaller than the unweighted mean of growth rates μ_R depending on the positive or negative correlations between the node biomasses and node growth rates.

TABLE 1

Table 1. Summary table: symbols used in the article and their descriptions.

Previous Eqs. (7–9) provide accurate computations of the regional, community or regional community growth rate if we know the growth rates of all involved nodes (either species or locations) (Figure 2). However, in many practical situations, we only can sample a limited number of nodes, $\tilde{n}$. We can then estimate the network growth rate $\widetilde{R_{net}}$ using Eq. (10) with the sampled nodes, replacing the population means (μ_N and μ_R) by the sample means (M_N and M_R), the population standard deviations (σ_N and σ_R) by the sample standard deviations (S_N and S_R), and the population correlation (c_N,R) by the sample correlation (C_N,R) (see Table 1). I.e., computing the biomass-weighted arithmetic mean of the sampled node growth rates. As the network estimate of growth rate corresponds to the weighed arithmetic mean of the node growth rates, we can estimate the standard error that arises from a partial (but representative) sampling using the formulas provided by Cochran (1977) and validated by Gatz and Smith (1995) (see Supplementary Appendix A)

$\begin{array}{c}\text{SE}\left(\widetilde{R_{net}}\right)=\hfill \\ t_{\tilde{n}-1}\frac{S_R}{\sqrt{\tilde{n}-1}}\sqrt{\mathrm{\hspace{0.25em}1}+\left(1-C_{N,R}^2\right){\left(\frac{S_N}{M_N}\right)}^2+C_{N,R}^2{\left(\frac{S_N}{M_N}\right)}^4},\hfill \end{array}(11)$

where $\text{SE}\left(\widetilde{R_{\text{net}}}\right)$ is the standard error of the network, corresponding to a 95% confidence level; $\tilde{n}$ is the sample size (i.e., the number of sampled nodes); and t_ñ–1 is the Student’s t distribution with $\tilde{n}-1$ degrees of freedom associated with a 95% confidence level, whose value is approximately 1.96 for large enough sampling sizes. Eq. (11) shows that the uncertainty in the determination of the regional community growth rate $\text{SE}\left(\widetilde{R_{\text{net}}}\right)$, as expected, decreases when sampling more localities or species (larger $\tilde{n}$), and when the growth rates vary less across localities or species (smaller S_R). Note that, generally, S_N≤M_N, and so also biomasses that vary less relative to their average biomass will lead to less variable estimates of R_net. Another implication is that positive or negative correlations between biomass and growth rate will mostly decrease $\text{SE}\left(\widetilde{R_{\text{net}}}\right)$, because S_N/M_N≤1 and so (S_N/M_N)⁴≪(S_N/M_N)². Nevertheless, since negative correlations decreases the network growth rate (Eq. 10), negative correlations between N and R require larger sample sizes to control the relative error on the estimated network growth rate (Figure 3).

FIGURE 3

Figure 3. Sample size required to estimate the growth rate of an ecological or spatial network from estimates at the nodes such that the standard error of the estimate is smaller than 10%. The required sample size mainly increases with the coefficient of variation of the node estimates of growth rate (S_R/M_R), and to a lesser extent with the coefficient of variation of the node biomass (S_N/M_N). Larger sampling sizes are required for more negative values of the cross-correlation between the biomass and the growth rate C_N,R (Table 1).

Resistance

The resistance of species i at location x, defined as in Eq. (2), is

${\mathrm{\Omega }}_{x,i}\equiv \frac{N_{x,i}\left(t_0\right)}{N_{x,i}\left(t_0\right)-N_{x,i}\left(t_0+\mathrm{\delta }t\right)}.(12)$

Then, from its definition, the regional and community resistances are the harmonic means of local or species resistances, weighted by the local or species biomasses (Supplementary Appendix B and Figure 4),

$\begin{array}{c}\hfill {\mathrm{\Omega }}_i^{\left(ℛ\right)}={\left[\frac{{\sum }_x{N}_{x,i}\left(t_0\right)\frac{1}{{\mathrm{\Omega }}_{x,i}}}{{\sum }_x{N}_{x,i}\left(t_0\right)}\right]}^{-1},\hfill \\ \hfill {\mathrm{\Omega }}_x^{\left(𝒞\right)}={\left[\frac{{\sum }_i{N}_{x,i}\left(t_0\right)\frac{1}{{\mathrm{\Omega }}_{x,i}}}{{\sum }_i{N}_{x,i}\left(t_0\right)}\right]}^{-1}.\hfill \end{array}(13)$

FIGURE 4

Figure 4. Regional (A) and community (B) resistances, compared to local and species resistances and to their biomass weighted arithmetic mean (AM_w) (blue circles), unweighted harmonic mean (HM ) (yellow circles), and biomass weighted harmonic mean (HM_w) (green circles). AM_w, HM, and HM_w values closer to the identity line (black, dashed) estimate more precisely the regional (A) and community resistance (B). Black dots are the resistances of the individual localities (A) or species (B) on which the means were computed. Data were generated for 6,300 random communities of 10 competitors in 10-node random spatial networks (Supplementary Figure 1A and Supplementary Text), from the biomass decrease caused by a reduction of the species growth rates at all locations. The harmonic mean of resistances is always smaller than the arithmetic mean (see Supplementary Appendix C). Moreover, HM_w are found to be better estimators for the resistance, as expected from the results shown in the text.

Hence, for any ecological or spatial network the network resistance is the weighted harmonic mean of node resistances, weighted by the node biomasses (see Supplementary Appendix B). Again, network resistance can be also rewritten in a more general way as

${\mathrm{\Omega }}_{net}=\frac{1}{{\mathrm{\mu }}_{{\mathrm{\Omega }}^{-1}}}\frac{1}{1+\frac{{\mathrm{\sigma }}_N}{{\mathrm{\mu }}_N}\frac{{\mathrm{\sigma }}_{{\mathrm{\Omega }}^{-1}}}{{\mathrm{\mu }}_{{\mathrm{\Omega }}^{-1}}}{c}_{N,{\mathrm{\Omega }}^{-1}}},(14)$

where again μ represents unweighted arithmetic means, σ the population standard deviations, and c_N,Ω^–1 the correlation between biomasses and the inverse of resistance (see Table 1). Hence, it is possible also to estimate the network resistance from a partial sampling of the network, replacing in Eq. (14) the population means, standard deviations and correlations by their sample equivalents.

Even though the resistance of a network Ω_net is the weighted harmonic mean of resistances, the network estimate of its inverse (Ω⁻¹)_net is the weighted arithmetic mean of the sub-units estimates of Ω⁻¹. This result again allows us to estimate its standard error arising from incomplete but representative network sampling (Supplementary Appendix A), which let us to obtain an expression for the standard error of the resistance obtained from a partial sampling of a network (Supplementary Appendix B, Eq. B8). The relative uncertainty of the network resistance will be dominated by the number of samples from the network, and by the variance of the inverse of resistances.

Initial Resilience

In Eq. (4), we show that initial resilience ρ is given by the product of resistance Ω and growth rate R, i.e., ρ = ΩR. Therefore, to estimate ρ for different scales one can use the already obtained scaling of Ω and R (which are scale invariant and have the simple estimators described). For example, defining the local species initial resilience as

${\mathrm{\rho }}_{x,i}\equiv {\frac{1}{\left|N_{x,i}\left(t\right)-N_{x,i}^{*}\right|}\frac{\text{d}\left(N_{x,i}\left(t\right)-N_{x,i}^{*}\right)}{\text{d}t}|}_{t=t_0},(15)$

and by defining the regional species initial resilience as the initial resilience of the regional biomass of a given species, it is easy to see that it coincides with the product of the regional resistance and the regional growth rates

${\mathrm{\rho }}_i^{\left(ℛ\right)}\equiv {\frac{1}{\left|N_i\left(t\right)-N_i^{*}\right|}\frac{\text{d}\left(N_i\left(t\right)-N_i^{*}\right)}{\text{d}t}|}_{t=t_0}={\mathrm{\Omega }}_i^{\left(ℛ\right)}{R}_i^{\left(ℛ\right)}.(16)$

That is, we can express the regional initial resilience as the product of the regional estimates of resistance (the weighted harmonic mean of local resistances) and growth rate (the weighted arithmetic mean of local growth rates). An analogous result links the community initial resilience to the community estimates of resistance and growth rate. In particular, since both resistance and growth rates were scale invariant (the network estimates of R and Ω were, respectively, the harmonic and arithmetic mean of the node estimates), also the initial resilience would be scale invariant. Actually, the regional initial resilience can be also expressed as

${\mathrm{\rho }}_i^{\left(ℛ\right)}=\frac{{\sum }_x{N}_{x,i}^{*}{\mathrm{\Omega }}_{x,i}^{-1}{\mathrm{\rho }}_{x,i}}{{\sum }_x{N}_{x,i}^{*}{\mathrm{\Omega }}_{x,i}^{-1}}.(17)$

Thus, the network estimate of initial resilience is the arithmetic mean of the node estimates of initial resilience, weighted by the change of the node biomass caused by a perturbation, NΩ⁻¹ (Figure 5). This result reinforces that initial resilience is another scale invariant stability property of ecosystems, and that less resistant nodes with higher biomass will disproportionally influence the initial resilience of the total network. Again, for any spatial or ecological network, we can express the network estimate of the initial resilience in a general way, as the product of the network resistance (the weighted harmonic mean of node resistances, highly influenced by nodes with higher biomasses and lower resistances) and the network growth rate:

$\begin{array}{c}{\mathrm{\rho }}_{net}={\mathrm{\Omega }}_{net}{R}_{net}\hfill \\ ={\mathrm{\mu }}_{\mathrm{\rho }}\left(1+\frac{\frac{{\mathrm{\sigma }}_N}{{\mathrm{\mu }}_N}\frac{{\mathrm{\sigma }}_{\mathrm{\rho }}}{{\mathrm{\mu }}_{\mathrm{\rho }}}{c}_{N,\mathrm{\rho }}+\frac{{\mathrm{\sigma }}_{\mathrm{\rho }}}{{\mathrm{\mu }}_{\mathrm{\rho }}}\frac{{\mathrm{\sigma }}_{{\mathrm{\Omega }}^{-1}}}{{\mathrm{\mu }}_{{\mathrm{\Omega }}^{-1}}}{c}_{\mathrm{\rho },{\mathrm{\Omega }}^{-1}}}{1+\frac{{\mathrm{\sigma }}_N}{{\mathrm{\mu }}_N}\frac{{\mathrm{\sigma }}_{{\mathrm{\Omega }}^{-1}}}{{\mathrm{\mu }}_{{\mathrm{\Omega }}^{-1}}}{c}_{N,{\mathrm{\Omega }}^{-1}}}\right).\hfill \end{array}(18)$

FIGURE 5

Figure 5. Regional (A) and community (B) initial resiliences, compared to local and species initial resiliences and to their unweighted arithmetic mean (AM ) (yellow circles) and biomass weighted arithmetic mean (AM_w) (green circles) (with weights given by the difference between the actual and the equilibrium biomasses, Eq. 17). AM and AM_w values closer to the identity line (black, dashed) estimate more precisely the regional (A) and community initial resilience (B). Black dots are the initial resiliences of the individual localities (A) or species (B) on which the means were computed. Data were generated for 6,300 random communities of 10 competitors in 10-node random spatial networks (Supplementary Figure 1A and Supplementary Text), after a biomass decrease from the equilibrium affecting all species at all locations. AM_w are found to be better estimators of the initial resilience, as expected from the results shown in the text.

This result indicates that correlations between the node biomass, resistances and initial resiliences can make the network initial resilience higher or smaller than the unweighted mean of node initial resilience estimates.

As for growth rate and resistance, we can also estimate network initial resilience from an incomplete sampling of the network as $\widetilde{{\mathrm{\rho }}_{net}}=\widetilde{{\mathrm{\Omega }}_{net}}\widetilde{R_{net}}$. And assuming no correlations between resistance and growth rates, the standard error committed with such sampling would be $\text{SE}\left(\widetilde{{\mathrm{\rho }}_{net}}\right)=\widetilde{{\mathrm{\Omega }}_{net}}\text{SE}\left(\widetilde{R_{net}}\right)+\widetilde{R_{net}}\text{SE}\left(\widetilde{{\mathrm{\Omega }}_{net}}\right)$. In such a case, the relative uncertainty of the estimated initial resilience would be first given by the sample size, and then by the variances of the growth rates and reciprocal resistances.

Invariability

We consider the invariability definition of Eq. (5), so the local species invariability reads

$I_{x,i}\equiv \frac{{\left[{\text{mean}}_t\left(N_{x,i}\left(t\right)\right)\right]}^2}{{\text{var}}_t\left(N_{x,i}\left(t\right)\right)}.(19)$

The regional invariability can be defined as the invariability of the total biomass of one species across all locations in a spatial network. For the synchronous space (“ss”) case, for which the local biomass dynamics are perfectly positively correlated, the regional invariability of species i reads (see Supplementary Appendix D)

$I_i^{(ℛ;ss)}={\left[\frac{{\sum }_x{N}_{x,i}^{*}}{{\sum }_x{N}_{x,i}^{*}\frac{1}{\sqrt{I_{x,i}}}}\right]}^2.(20)$

Then, for perfectly synchronous local dynamics, the regional species invariability is the square of the harmonic mean of the square root of the local species invariabilities, weighted by the equilibrium local biomass densities. Conversely, if the local species biomass dynamics is spatially asynchronous (asynchronous space, “as”), the regional species invariability of the whole spatial network is

$\begin{array}{c}{I}_i^{(ℛ;as)}=\frac{{\left({\sum }_x{N}_{x,i}^{*}\right)}^2}{{\sum }_x{\left(N_{x,i}^{*}\right)}^2\frac{1}{I_{x,i}}}\hfill \\ =n_L\frac{1}{1+{\left(\frac{{\mathrm{\sigma }}_{N_i^{*}}}{{\mathrm{\mu }}_{N_i^{*}}}\right)}^2}\frac{{\sum }_x{\left(N_{x,i}^{*}\right)}^2}{{\sum }_x{\left(N_{x,i}^{*}\right)}^2\frac{1}{I_{x,i}}}.\hfill \end{array}(21)$

For the asynchronous-space case, the regional invariability is proportional to the number of locations n_L (Figure 6). I.e., for the asynchronous-space case, invariability would be an extensive stability property, that grows linearly with the size of the system. In this asynchronous case, invariability is also proportional to the harmonic mean of local species invariabilities weighted by the squared local species biomasses. In addition, it is modulated by the spatial variance ${\mathrm{\sigma }}_{N_i^{*}}^2$ and mean ${\mathrm{\mu }}_{N_i^{*}}$ of the local equilibrium biomasses of species i. It can be proven that invariability is higher for asynchronous than for synchronous dynamics (see Supplementary Appendix D). Moreover, the number of locations does not modify the local invariability estimates, and there are not significative differences between cases with synchronous or asynchronous dynamics (Figure 6B).

FIGURE 6

Figure 6. Local (A) and regional (B) invariability estimates in random spatial networks of random communities of 10 competitor species, for different sizes of the spatial network; and species (C) and community (D) invariability estimates of random communities of competitors at 10-node random spatial networks, for different number of species forming the communities. In (A,B), we have considered three different scenarios: asynchronous local dynamics ($\overline{c}=0$, yellow box plots), partially synchronous local dynamics ($\overline{c}=0.36$, blue), and perfectly synchronous local dynamics ($\overline{c}=1$, green). In (C,D), we have considered the cases of asynchronous ($\tilde{c}=0$), partially synchronous ($\tilde{c}=0.17$) and perfectly synchronous ($\tilde{c}=1$) species dynamics. Solid lines depict the invariabilities predicted by analytical expressions (Eqs. 20–22), and analogous expressions for community invariability).

The more general case of not perfectly synchronous dynamics can be expressed as

$I_i^{\left(ℛ\right)}={\left[\left(1-{\overline{c}}_i\right)\frac{1}{I_i^{(ℛ;as)}}+{\overline{c}}_i\frac{1}{I_i^{(ℛ;ss)}}\right]}^{-1},(22)$

where ${\overline{c}}_i$ is the typical correlation between different locations (Eq. D10 of the Supplementary Appendix D). For the case ${\overline{c}}_i>0$, and since by definition ${\overline{c}}_i\le 1$, $I_i^{\left(ℛ\right)}$ would simply be the weighted harmonic mean of $I_i^{(ℛ;as)}$ and $I_i^{(ℛ;ss)}$, with weights equal to $1-{\overline{c}}_i$ and ${\overline{c}}_i$, respectively. And since, even though $I_i^{(ℛ;ss)}$ does not depend on the number of locations n_L (Eq. 20), $I_i^{(ℛ;as)}$ increases with n_L (Eq. 21), the resulting regional species invariability would increase as well with n_L (except for the special case ${\overline{c}}_i=1$). For the case ${\overline{c}}_i<0$, since $I_i^{(ℛ;as)}>I_i^{(ℛ;ss)}$, the regional species invariability $I_i^{\left(ℛ\right)}$ would be larger than $I_i^{(ℛ;as)}$. Since $I_i^{(ℛ;as)}$ increases linearly with the number of locations, the regional species invariability would then also increase with the number of locations. In summary, when the local population dynamics are not perfectly synchronized (so the typical spatial correlation of the local biomasses $\overline{c}$ is less than 1), the regional invariability increases with the number of locations of the spatial network (Figure 6).

For community invariability, we can obtain completely analogous expressions to Eqs. (20–22), In particular, this proves that community invariability increases with the number of species forming the community, except for the special case of perfectly synchronous dynamics across species (Figure 6C), while the degree of synchrony and the number of species do not significantly affect the invariabilities at the species level (Figure 6D).

In general, network invariability is not a mean of the invariability estimates at the network nodes, so we cannot estimate its standard error in the same way that we did for resistance, growth rate and initial resilience (Supplementary Appendix A). We did not pursue here the characterization of such network invariability standard error. To estimate the error that arises from incomplete network sampling, general bootstrapping techniques should be applied instead (Efron and Tibshirani, 1985; Hesterberg, 2011).

Model Simulations

In this study, we have investigated how different stability components such as growth rate, initial resilience, resistance, and invariability scale from the local or species level to the regional or community level. We now compare these scaling laws to numerically simulated population dynamics of a community of 10 competitors with the Lotka-Volterra model (see Supplementary Appendix E) in 10-node random spatial networks (Supplementary Figure 1A and Figures 2, 4–6). To ensure that the results do not depend on the chosen network, and motivated by fundamental differences of meta-community stability between linear and riverine networks (Fagan, 2002; Carrara et al., 2012; Altermatt, 2013; Liu et al., 2013; Peterson et al., 2013), we complement these results with results for realistic riverine dendritic networks (Supplementary Figure 1B) generated in R (R Core Team, 2020) with the OCNet package (Carraro et al., 2020; Supplementary Figures 2–5). All simulations of community dynamics were done in Python 3.7 (Python Core Team, 2019).

The simulation results confirm our theoretical prediction that growth rate and initial resilience are scale-free stability properties, where regional and community estimates equal to the weighted arithmetic mean of the estimates at the local or species level (Figures 2, 5 and Supplementary Figures 2, 4). Also, the simulations confirm that resistance is another scale-free property: the regional and community estimates of resistance are the harmonic mean of the local and species resistance estimates, weighted by the local biomasses or the species proportions (Figure 4 and Supplementary Figure 3). The numerical simulations also confirmed that invariability is a scale-free property solely in networks with perfectly synchronous dynamics for which all sub-units effectively act as a unique single unit (species or location). In more realistic networks, with imperfect synchrony across subunits, the invariability is higher than for the perfectly synchronous case (Figure 6 and Supplementary Figure 5), and it increases with the network size, so the regional or community invariability is actually larger than the average of its elements, and this difference is more pronounced in larger networks. Thus, realistic spatial networks are more invariable than their individual locations, and community dynamics are more invariable than the population dynamics of the species forming the community (Loreau and De Mazancourt, 2008; Haegeman et al., 2016).

Discussion

We have shown that resistance and initial resilience (and growth rate) of ecological or spatial networks, unlike invariability, are biomass-weighted means of the estimates of these stability measures at the nodes of the network. In this section, we will discuss the consequences of this fundamental difference between these stability components.

Resistance and Initial Resilience Are Scale-Free Network Properties, While Invariability Is Not

Some stability components, such as invariability, have been found to increase with the ecological (Thébault and Loreau, 2005; Tilman et al., 2006; Gross et al., 2014) and spatial scale (Wang and Loreau, 2014; Wang et al., 2019), as a consequence of not perfectly synchronous dynamics among species and among locations (Doak et al., 1998). Hence, communities and regions are more invariable than their constituent species and locations. On the contrary, other stability components seem to not depend on scale (Haegeman et al., 2016). To solve this issue, Clark et al. (2021) have proposed different scaling laws of three different common stability components: resistance, initial resilience, and invariance.

Our analysis confirms that regional and community invariability is larger than local and species invariability, and generally increases with the size of the studied network (Figure 6). Similar results were obtained by Wang and Loreau (2014), who showed that the regional variability decreases with the species richness and the region size. However, and as is the case for asymptotic resilience (Haegeman et al., 2016), network resistance and short-term or initial growth rate and resilience [independent of asymptotic resilience, and a better proxy of resilience in experiments (Arnoldi et al., 2016, 2018)] is the mean of the local and species values (Figures 2, 4, 5). As this result is a consequence of the mathematical definition of these stability components, it will hold for any spatial and ecological network of any complexity, and for any meta-community dynamics model (Supplementary Figures 2–5).

These results contribute to a better understanding of the multidimensional nature of ecological stability. While stability properties can be correlated (Donohue et al., 2013), depending on the characteristic of the environmental fluctuations affecting the systems (Arnoldi et al., 2019; Radchuk et al., 2019), their scaling laws can introduce another axis of fundamental differentiation between different stability components. Indeed, one could distinguish between network-level stability components (those fundamentally depending on the topology and size of the ecological network) and node-level stability components (those reflecting network averages of the node-level estimates). As a consequence, the analysis of multiple components of stability of the ecosystems might be preferred to the employment of single metrics that aim to reproduce the whole ecosystem stability (Lemoine, 2020).

Resistance Is More Affected Than Initial Resilience by the Presence of Low-Stable Nodes or Species

Although both resistance and initial resilience are scale-free properties, they differ in how the network estimate is averaged from the node measures, which has important ecological consequences. Harmonic means are more affected by the presence of low numbers, and less affected by the presence of high numbers, than arithmetic means (Ferger, 1931). Hence, the presence of less resistant species and locations will affect network-level resistance much more than network-level initial resilience (see Supplementary Appendix C). High resilient nodes can easily compensate low resilient nodes, as such bolstering network initial resilience. This will lower the impact of stressors that only affect a fraction of the nodes (Supp and Ernest, 2014). However, this is not the case for resistance (Figure 7). Low-resistant nodes limit the resistance of a network much more, which makes resistance a stability component that is more difficult to protect in ecological networks. For example, in spatially heterogeneous meta-populations, meta-community resistance will be mainly determined by the resistance of the less stable regions, while the meta-community initial resilience will be mostly given by the average spatial conditions. Thus, the heterogenous presence of stressors in the meta-community (McCluney et al., 2014) is expected to have a stronger effect on the network resistance than in the network initial resilience.

FIGURE 7

Figure 7. Schematic comparison of different stability components at local or species level vs. at regional or community level, assuming normal distributions for the local and species estimates. The regional or community growth rate (A) and initial resilience (B) is the weighted arithmetic mean of the estimates at the local or species levels. On the contrary, regional/community resistance (C) is the weighted harmonic mean of the local/species estimates, so locations or species with low stability will limit the resistance of spatial or ecological networks.

Influence of Mathematical Definitions of Stability

In this study, we have shown how different stability components scale from the local and species level to the regional and community level. Starting from common mathematical definitions, we showed that resistance and initial resilience are scale-free properties, while regions and communities are fundamentally more invariable than local species population dynamics. However, we anticipate that this result will depend on the employed mathematical definition (and then, on the proposed measurements) for these stability properties.

As previously noted, there is evidence that communities and spatial networks are more invariable than local species populations, as a consequence of imperfect synchronization on the local population dynamics (Gross et al., 2014; Wang et al., 2019). In this article, we have obtained the same result: except for the unrealistic case of perfect inter-species and inter-location synchrony, meta-communities are more invariable than local populations, and such meta-community invariability increases with the number of species and the size of the spatial network (Figure 6). Clark et al. (2021) have shown that invariance (the inverse of the biomass variance, not normalized by the average species biomasses) decreases with scale. Here, we show that invariability (i.e., biomass-normalized) increases with scale. Repeating our analyses using invariance recovers Clark et al.’s (2021) results (Supplementary Appendix F). Overall, this shows that normalizing the variance by the mean biomass (which increases with the network size) has a large influence on how we appreciate the scaling of stability. As networks will often be less variable than their nodes, we advocate the use of normalized stability properties related to variability when studying the effect of scale on this kind of stability. A similar difference between our results and those by Clark et al. (2021) on the scale dependence of resistance can be also explained from biomass normalization (Supplementary Appendix G).

For initial resilience, Clark et al. (2021) already employed a normalized definition. Here, with an analogous definition, we showed that the network initial resilience is the arithmetic mean of local species initial resiliences, being independent of the network size (Figure 5). Clark et al. (2021) also found that the median of the initial resiliences was proportional to the ratio of the expected values of the invariance and the resistance for linear models. Since their defined invariance and resistance decrease with the characteristic ecological or spatial scale, their scale dependencies cancel out. An interesting question beyond the purpose of the present work is to what extent this relation between stability measures holds in the non-linear case.

The different means and behavior between resistance and initial resilience, discussed in section “Resistance Is More Affected Than Initial Resilience by the Presence of Low-Stable Nodes or Species,” depend on their mathematical definition and on the distribution of those stability components. For example, instead of resistance Ω (inverse of the relative change in biomass after a perturbation) we can define an alternative stability measure just given by the relative change of the biomass after a perturbation, i.e., Ω⁻¹. More resistant systems present smaller values of Ω⁻¹, and Ω⁻¹ represents the plasticity of the system against perturbations. Since the harmonic mean of a random variable is the inverse of the arithmetic mean of the reciprocals, it is easy to prove that a network estimator of Ω⁻¹ would simply be the weighted arithmetic mean of the estimates at the nodes. For this new defined resistance, the presence of outliers affects the network resistance in the same way than the presence of outliers affected the network initial resilience, so nodes with above-average values of Ω⁻¹ can be easily compensated by nodes with below-average values of Ω⁻¹, having a limited effect on the network-level estimate of Ω⁻¹. This is a clear indication of how the heterogeneous distribution of local species estimates (particularly its skewness Stevens, 1955) affects the network resistance and initial resilience.

The scale-free property found for the growth rate R, the resistance Ω, and the initial resilience ρ is due to their character as intensive quantities. The total biomass N, its derivative $\frac{\text{d}N}{\text{d}t}$, and the change in biomass due to a perturbation N(t₀)−N(t₀δ + t), are are additive for subsystems. Their quotients have allowed us to construct quantities independent of the extent of the system, i.e., scale-free quantities. Namely, growth rate R, resistance, Ω and initial resilience ρ. For the simpler cases, the growth rate R and the inverse of the resistance Ω⁻¹ are given by a biomass-weighted arithmetic mean, which compensates the total biomass increase as the considered scale increases. The expression of Ω as a biomass-weighted harmonic mean is equivalent to the expression of Ω⁻¹ as a biomass-weighted arithmetic mean and conserves the scale-free properties. The scale-free property of initial resilience ρ can then be seen as a consequence of being the product (or quotient) of two intensive (or scale-free) quantities.

This view also shows why, in general, invariability is not scale-free. The temporal variance var_t(N(t)) is not extensive, because var_t(N) = E_t[(N−E_t[N])²] = E_t[N²]−(E_t[N])² is not additive for subsystems. Neither $\sqrt{{\text{var}}_t\left(N\right)}$ is extensive in general. This makes that only for completely synchronous dynamics the invariability is scale-free, as previously shown.

Implications for Measuring Stability Empirically

Resistance and initial resilience of ecological spatial networks are biomass-weighted means of the local species estimates at the nodes of the networks, so they can be easily estimated from partial samples of the network. This property is important for the assessment of stability in large experiments (De Raedt et al., 2019; Karakoç et al., 2020; Saade et al., 2020) or field campaigns. Using our equations for the relative standard error of these stability indices (Eq. 11, and Eq. A5 of Supplementary Appendix A), one is able to estimate the sampling size required to control the error committed in the estimation of the network stability components from partial network samplings (illustrations in Figure 3 and Supplementary Figure 6). The coefficient of variation of the studied stability property (resistance or initial resilience) affects this required sampling effort most (Figure 3). Thus, the coefficient of variation will be higher for networks with more variable stability. Two other factors influence this variation: (1) a greater coefficient of variation of the biomasses; (2) a negative correlation between the biomass and the stability (see Supplementary Figure 6).

With respect to invariability, the standard error associated with an incomplete sampling is more difficult to estimate, since generally the network invariability is not a mean of the nodes’ invariabilities, and depends on the size of the network. Hence, for this stability component the standard error should generally be assessed directly with a bootstrap. Moreover, for controlling the error associated with the estimation of network invariability from node-level invariability, it would be important to have an unbiased estimate of network size.

Implications for the Stability-Complexity Debate

The stability-complexity debate (McCann, 2000; Allesina and Tang, 2015) originated from the disagreement between experimental observations often finding more complex systems to be more stable (Ives and Carpenter, 2007), and theoretical analyses finding more complex systems to be less stable (May, 1972; Pimm, 1984). To solve this disagreement, some authors have proposed generalizations of the original work of May (1972) that account for non-random among-species interactions (Yodzis, 1981; Rooney et al., 2006), or the stabilizing role of dispersal and spatial heterogeneity (Plitzko and Drossel, 2015; Gravel et al., 2016). Other approaches have suggested that the disagreement is caused by a focus on asymptotic resilience in theoretical studies (Pimm, 1984; McCann, 2000; Saeedian et al., 2022), which is biased by rare species (Haegeman et al., 2016). Our results adhere to this point of view, by showing that more complex systems (i.e., ecological and spatial networks, as opposed to single species and locations) are inevitably less variable, if using normalized estimates that correct for inherent effects on system size of system complexity. Consequently, such correction for system size leads to no relationship between complexity and the other stability properties (resistance and initial resilience). Taken together, these findings confirm that using a sole stability component (e.g., as asymptotic resilience) does not fully capture the complex ways in which biological systems deal with environmental changes (Pennekamp et al., 2018; Arnoldi et al., 2019). Assessing stability from a multi-dimensional perspective (Donohue et al., 2013; Arnoldi et al., 2019; Radchuk et al., 2019) will provide a more comprehensive picture and can reconcile apparent contradictions between and among theoretical and empirical studies.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

JJ and FDL conceived the presented idea. JJ performed the numerical simulations and the analytical computations, wrote the first draft, and prepared the figures. FJC-G and FDL verified the analytical derivations. All authors discussed the results and contributed to the final manuscript.

Funding

This work was supported by the CEFIC under LRI project ECO50, and the special research fund (FSR) from UNamur. Computational resources have been provided by the Consortium des Équipements de Calcul Intensif (CÉCI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under (Grant No. 2.5020.11) and by the Walloon Region. FJC-G was funded by the European Regional Development Fund (ERDF) and by the Spanish Ministry of Economy and Competitiveness through (Grant No. RTI2018-095802-B-I00), by European Union’s Horizon 2020 through (grant agreement No. 817578 TRIATLAS), and FDL acknowledges support from his Namur Research College fellowship, granted by the University of Namur.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We acknowledge constructive discussions on these results with Jeff Arnoldi, Bart Haegeman, Michel Loreau, and other scientists of the Ecological Station of Moulis (CNRS, France). FJC-G acknowledges the warm welcome at the Ecological Station of Moulis (CNRS, France). We also acknowledge the reviewers for their valuable feedback and efforts toward improving this manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2022.861537/full#supplementary-material

References

Allesina, S., and Tang, S. (2015). The stability–complexity relationship at age 40: a random matrix perspective. Popul. Ecol. 57, 63–75. doi: 10.1007/s10144-014-0471-0

CrossRef Full Text | Google Scholar

Altermatt, F. (2013). Diversity in riverine metacommunities: a network perspective. Aquat. Ecol. 47, 365–377. doi: 10.1007/s10452-013-9450-3

CrossRef Full Text | Google Scholar

Amarasekare, P. (2008). Spatial dynamics of foodwebs. Annu. Rev. Ecol. Evol. Syst. 39, 479–500. doi: 10.1146/annurev.ecolsys.39.110707.173434

CrossRef Full Text | Google Scholar

Arnoldi, J.-F., Bideault, A., Loreau, M., and Haegeman, B. (2018). How ecosystems recover from pulse perturbations: a theory of short- to long-term responses. J. Theor. Biol. 436, 79–92. doi: 10.1016/j.jtbi.2017.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Arnoldi, J. F., Loreau, M., and Haegeman, B. (2016). Resilience, reactivity and variability: a mathematical comparison of ecological stability measures. J. Theor. Biol. 389, 47–59. doi: 10.1016/j.jtbi.2015.10.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Arnoldi, J. F., Loreau, M., and Haegeman, B. (2019). The inherent multidimensionality of temporal variability: how common and rare species shape stability patterns. Ecol. Lett. 22, 1557–1567. doi: 10.1111/ele.13345

PubMed Abstract | CrossRef Full Text | Google Scholar

Baert, J. M., De Laender, F., Sabbe, K., and Janssen, C. R. (2016). Biodiversity increases functional and compositional resistance, but decreases resilience in phytoplankton communities. Ecology 97, 3433–3440. doi: 10.1002/ecy.1601

PubMed Abstract | CrossRef Full Text | Google Scholar

Carrara, F., Altermatt, F., Rodriguez-Iturbe, I., and Rinaldo, A. (2012). Dendritic connectivity controls biodiversity patterns in experimental metacommunities. Proc. Natl. Acad. Sci. U.S.A. 109, 5761–5766. doi: 10.1073/pnas.1119651109

PubMed Abstract | CrossRef Full Text | Google Scholar

Carraro, L., Bertuzzo, E., Fronhofer, E. A., Furrer, R., Gounand, I., Rinaldo, A., et al. (2020). Generation and application of river network analogues for use in ecology and evolution. Ecol. Evol. 10, 7537–7550. doi: 10.1002/ece3.6479

PubMed Abstract | CrossRef Full Text | Google Scholar

Chase, J. M., Blowes, S. A., Knight, T. M., Gerstner, K., and May, F. (2020). Ecosystem decay exacerbates biodiversity loss with habitat loss. Nature 584, 238–243. doi: 10.1038/s41586-020-2531-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Chesson, P. (2000). General theory of competitive coexistence in spatially-varying environments. Theor. Popul. Biol. 58, 211–237. doi: 10.1006/tpbi.2000.1486

PubMed Abstract | CrossRef Full Text | Google Scholar

Clark, A. T., Arnoldi, J. F., Zelnik, Y. R., Barabas, G., Hodapp, D., Karakoç, C., et al. (2021). General statistical scaling laws for stability in ecological systems. Ecol. Lett. 24, 1474–1486. doi: 10.1111/ele.13760

PubMed Abstract | CrossRef Full Text | Google Scholar

Cochran, W. G. (1977). Sampling Techniques. 3rd Edn. New York, NY: John Wiley & Sons.

Google Scholar

De Raedt, J., Baert, J. M., Janssen, C. R., and De Laender, F. (2019). Stressor fluxes alter the relationship between beta-diversity and regional productivity. Oikos 128, 1015–1026. doi: 10.1111/oik.05191

CrossRef Full Text | Google Scholar

Doak, D. F., Bigger, D., Harding, E. K., Marvier, M. A., O’Malley, R. E., and Thomson, D. (1998). The statistical inevitability of stability-diversity relationships in community ecology. Am. Nat. 151, 264–276. doi: 10.1086/286117

PubMed Abstract | CrossRef Full Text | Google Scholar

Domínguez-García, V., Dakos, V., and Kéfi, S. (2019). Unveiling dimensions of stability in complex ecological networks. Proc. Natl. Acad. Sci. U.S.A. 116, 25714–25720. doi: 10.1073/pnas.1904470116

PubMed Abstract | CrossRef Full Text | Google Scholar

Donohue, I., Petchey, O. L., Montoya, J. M., Jackson, A. L., Mcnally, L., Viana, M., et al. (2013). On the dimensionality of ecological stability. Ecol. Lett. 16, 421–429. doi: 10.1111/ele.12086

PubMed Abstract | CrossRef Full Text | Google Scholar

Downing, A. L., Brown, B. L., and Leibold, M. A. (2014). Multiple diversity-stability mechanisms enhance population and community stability in aquatic food webs. Ecology 95, 173–184. doi: 10.1890/12-1406.1

CrossRef Full Text | Google Scholar

Efron, B., and Tibshirani, R. (1985). The bootstrap method for assessing statistical accuracy. Behaviormetrika 12, 1–35. doi: 10.2333/bhmk.12.17_1

CrossRef Full Text | Google Scholar

Fagan, W. F. (2002). Connectivity, fragmentation, and extinction risk in dendritic metapopulations. Ecology 83, 3243–3249.

Google Scholar

Ferger, W. F. (1931). The Nature and Use of the Harmonic Mean. J. Am. Stat. Assoc. 26, 36–40. doi: 10.1080/01621459.1931.10503148

CrossRef Full Text | Google Scholar

Flöder, S., and Hillebrand, H. (2012). Species traits and species diversity affect community stability in a multiple stressor framework. Aquat. Biol. 17, 197–209. doi: 10.3354/ab00479

CrossRef Full Text | Google Scholar

Gatz, D. F., and Smith, L. (1995). The standard error of a weighted mean concentration-I. Bootstrapping vs other methods. Atmos. Environ. 29, 1185–1193. doi: 10.1016/1352-2310(94)00210-C

CrossRef Full Text | Google Scholar

Gravel, D., Massol, F., and Leibold, M. A. (2016). Stability and complexity in model meta-ecosystems. Nat. Commun. 7:12457. doi: 10.1038/ncomms12457

PubMed Abstract | CrossRef Full Text | Google Scholar

Greig, H. S., McHugh, P. A., Thompson, R. M., Warburton, H. J., and McIntosh, A. R. (2022). Habitat size influences community stability. Ecology 103:e03545. doi: 10.1002/ecy.3545

PubMed Abstract | CrossRef Full Text | Google Scholar

Grimm, V., and 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. doi: 10.1007/s004420050090

PubMed Abstract | CrossRef Full Text | Google Scholar

Gross, K., Cardinale, B. J., Fox, J. W., Gonzalez, A., Loreau, M., Wayne Polley, H., et al. (2014). Species richness and the temporal stability of biomass production: a new analysis of recent biodiversity experiments. Am. Nat. 183, 1–12. doi: 10.1086/673915

PubMed Abstract | CrossRef Full Text | Google Scholar

Haegeman, B., Arnoldi, J.-F., Wang, S., de Mazancourt, C., Montoya, J., and Loreau, M. (2016). Resilience, invariability, and ecological stability across levels of organization. bioRxiv[Preprint]. doi: 10.1101/085852

CrossRef Full Text | Google Scholar

Hesterberg, T. (2011). Bootstrap. Wiley Interdiscip. Rev. Comput. Stat. 3, 497–526. doi: 10.1002/wics.182

CrossRef Full Text | Google Scholar

Hillebrand, H., Langenheder, S., Lebret, K., Lindström, E., and Östman, Ö, and Striebel, M. (2018). Decomposing multiple dimensions of stability in global change experiments. Ecol. Lett. 21, 21–30. doi: 10.1111/ele.12867

PubMed Abstract | CrossRef Full Text | Google Scholar

Isbell, F., Craven, D., Connolly, J., Loreau, M., Schmid, B., Beierkuhnlein, C., et al. (2015). Biodiversity increases the resistance of ecosystem productivity to climate extremes. Nature 526, 574–577. doi: 10.1038/nature15374

PubMed Abstract | CrossRef Full Text | Google Scholar

IUPAC (2019). The IUPAC Compendium of Chemical Terminology – The Gold Book, 2nd Edn, ed. V. Gold (Research Triangle Park, NC: International Union of Pure and Applied Chemistry (IUPAC)). doi: 10.1351/goldbook

CrossRef Full Text | Google Scholar

Ives, A. R., and Carpenter, S. R. (2007). Stability and diversity of ecosystems. Science 317, 58–62. doi: 10.1126/science.1133258

PubMed Abstract | CrossRef Full Text | Google Scholar

Ives, A. R., Klug, J. L., and Gross, K. (2000). Stability and species richness in complex communities. Ecol. Lett. 3, 399–411. doi: 10.1046/j.1461-0248.2000.00144.x

CrossRef Full Text | Google Scholar

Jarillo, J., Sæther, B.-E., Engen, S., and Cao, F. J. (2018). Spatial scales of population synchrony of two competing species: effects of harvesting and strength of competition. Oikos 127, 1459–1470. doi: 10.1111/oik.05069

CrossRef Full Text | Google Scholar

Jarillo, J., Sæther, B.-E., Engen, S., and Cao-García, F. J. (2020). Spatial scales of population synchrony in predator-prey systems. Am. Nat. 195, 216–230. doi: 10.1086/706913

PubMed Abstract | CrossRef Full Text | Google Scholar

Karakoç, C., Clark, A. T., and Chatzinotas, A. (2020). Diversity and coexistence are influenced by time-dependent species interactions in a predator–prey system. Ecol. Lett. 23, 983–993. doi: 10.1111/ele.13500

PubMed Abstract | CrossRef Full Text | Google Scholar

Kéfi, S., Domínguez-García, V., Donohue, I., Fontaine, C., Thébault, E., and Dakos, V. (2019). Advancing our understanding of ecological stability. Ecol. Lett. 22, 1349–1356. doi: 10.1111/ele.13340

PubMed Abstract | CrossRef Full Text | Google Scholar

Lande, R., Engen, S., and Sæther, B.-E. (1999). Spatial scale of population synchrony: environmental correlation versus dispersal and density regulation. Am. Nat. 154, 271–281. doi: 10.1086/303240

PubMed Abstract | CrossRef Full Text | Google Scholar

Leibold, M. A., Holyoak, M., Mouquet, N., Amarasekare, P., Chase, J. M., Hoopes, M. F., et al. (2004). The metacommunity concept: a framework for multi-scale community ecology. Ecol. Lett. 7, 601–613. doi: 10.1111/j.1461-0248.2004.00608.x

CrossRef Full Text | Google Scholar

Lemoine, N. P. (2020). Unifying ecosystem responses to disturbance into a single statistical framework. Oikos 130, 1–14. doi: 10.1111/oik.07752

CrossRef Full Text | Google Scholar

Levin, S. A. (1992). The problem of pattern and scale in ecology: the Robert H. Macarthur award lecture. Ecology 73, 1943–1967. doi: 10.2307/1941447

CrossRef Full Text | Google Scholar

Limberger, R., Pitt, A., Hahn, M. W., and Wickham, S. A. (2019). Spatial insurance in multi-trophic metacommunities. Ecol. Lett. 22, 1828–1837. doi: 10.1111/ele.13365

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Soininen, J., Han, B. P., and Declerck, S. A. J. (2013). Effects of connectivity, dispersal directionality and functional traits on the metacommunity structure of river benthic diatoms. J. Biogeogr. 40, 2238–2248. doi: 10.1111/jbi.12160

CrossRef Full Text | Google Scholar

Loreau, M., and De Mazancourt, C. (2008). Species synchrony and its drivers: neutral and nonneutral community dynamics in fluctuating environments. Am. Nat. 172, 48–66. doi: 10.1086/589746

PubMed Abstract | CrossRef Full Text | Google Scholar

May, R. M. (1972). Will a large complex system be stable? Nature 238, 413–414. doi: 10.1038/238413a0

PubMed Abstract | CrossRef Full Text | Google Scholar

McCann, K. S. (2000). The diversity-stability debate. Nature 405, 228–233. doi: 10.1038/35012234

PubMed Abstract | CrossRef Full Text | Google Scholar

McCluney, K. E., Poff, N. L., Palmer, M. A., Thorp, J. H., Poole, G. C., Williams, B. S., et al. (2014). Riverine macrosystems ecology: sensitivity, resistance, and resilience of whole river basins with human alterations. Front. Ecol. Environ. 12:48–58. doi: 10.1890/120367

CrossRef Full Text | Google Scholar

Moran, P. A. P. (1953). The statistical analysis of the Canadian lynx cycle. II. Synchronization and meteorology. Aust. J. Zool. 1, 291–298. doi: 10.1071/ZO9530291

CrossRef Full Text | Google Scholar

Mougi, A., and Kondoh, M. (2012). Diversity of interaction types and ecological community stability. Science 337, 349–351. doi: 10.1126/science.1220529

PubMed Abstract | CrossRef Full Text | Google Scholar

Mougi, A., and Kondoh, M. (2016). Food-web complexity, meta-community complexity and community stability. Sci. Rep. 6, 1–5. doi: 10.1038/srep24478

PubMed Abstract | CrossRef Full Text | Google Scholar

Neubert, M. G., and Caswell, H. (1997). Alternatives to resilience for measuring the responses of ecological systems to perturbations. Ecology 78, 653–665. doi: 10.1111/gcb.12845

PubMed Abstract | CrossRef Full Text | Google Scholar

Pennekamp, F., Pontarp, M., Tabi, A., Altermatt, F., Alther, R., Choffat, Y., et al. (2018). Biodiversity increases and decreases ecosystem stability. Nature 563, 109–112. doi: 10.1038/s41586-018-0627-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Peterson, E. E., Ver Hoef, J. M., Isaak, D. J., Falke, J. A., Fortin, M. J., Jordan, C. E., et al. (2013). Modelling dendritic ecological networks in space: an integrated network perspective. Ecol. Lett. 16, 707–719. doi: 10.1111/ele.12084

PubMed Abstract | CrossRef Full Text | Google Scholar

Pimm, S. L. (1984). The complexity and stability of ecosystems. Nature 307, 321–326. doi: 10.1038/307321a0

CrossRef Full Text | Google Scholar

Plitzko, S. J., and Drossel, B. (2015). The effect of dispersal between patches on the stability of large trophic food webs. Theor. Ecol. 8, 233–244. doi: 10.1007/s12080-014-0247-3

CrossRef Full Text | Google Scholar

Python Core Team (2019). Python: A Dynamic, Open Source Programming Language. Available online at: https://www.python.org/ (accessed December 2021).

Google Scholar

R Core Team (2020). R: A Language and Environment for Statistical Computing. Vienna: R foundation for statistical computing.

Google Scholar

Radchuk, V., Laender, F., De, Cabral, J. S., Boulangeat, I., Crawford, M., et al. (2019). The dimensionality of stability depends on disturbance type. Ecol. Lett. 22, 674–684. doi: 10.1111/ele.13226

PubMed Abstract | CrossRef Full Text | Google Scholar

Rooney, N., McCann, K., Gellner, G., and Moore, J. C. (2006). Structural asymmetry and the stability of diverse food webs. Nature 442, 265–269. doi: 10.1038/nature04887

PubMed Abstract | CrossRef Full Text | Google Scholar

Saade, C., Kéfi, S., Gougat-Barbera, C., Rosenbaum, B., and Fronhofer, E. A. (2020). Spatial distribution of local patch extinctions drives recovery dynamics in metacommunities. bioRxiv[Preprint]. doi: 10.1098/rspb.2022.0543.2020.12.03.409524

PubMed Abstract | CrossRef Full Text | Google Scholar

Saeedian, M., Pigani, E., Maritan, A., Suweis, S., and Azaele, S. (2022). Effect of delay on the emergent stability patterns in generalized Lotka–Volterra ecological dynamics. Phil. Trans. R. Soc. A, 380:20210245. (accessed December 2021). doi: 10.1098/rsta.2021.0245

PubMed Abstract | CrossRef Full Text | Google Scholar

Stevens, S. S. (1955). On the averaging of data. Science 121, 113–116. doi: 10.1126/science.121.3135.113

PubMed Abstract | CrossRef Full Text | Google Scholar

Supp, S. R., and Ernest, S. K. M. (2014). Species-level and community-level responses to disturbance: a cross-community analysis. Ecology 95, 1717–1723. doi: 10.1890/13-2250.1

CrossRef Full Text | Google Scholar

Thébault, E., and Loreau, M. (2005). Trophic interactions and the relationship between species diversity and ecosystem stability. Am. Nat. 166, E95–E114. doi: 10.1086/444403

PubMed Abstract | CrossRef Full Text | Google Scholar

Thibaut, L. M., and Connolly, S. R. (2013). Understanding diversity-stability relationships: towards a unified model of portfolio effects. Ecol. Lett. 16, 140–150. doi: 10.1111/ele.12019

PubMed Abstract | CrossRef Full Text | Google Scholar

Tilman, D., Reich, P. B., and Knops, J. M. H. (2006). Biodiversity and ecosystem stability in a decade-long grassland experiment. Nature 441, 629–632. doi: 10.1038/nature04742

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, S., Lamy, T., Hallett, L. M., and Loreau, M. (2019). Stability and synchrony across ecological hierarchies in heterogeneous metacommunities: linking theory to data. Ecography 42, 1200–1211. doi: 10.1111/ecog.04290

CrossRef Full Text | Google Scholar

Wang, S., and Loreau, M. (2014). Ecosystem stability in space: α, β and γ variability. Ecol. Lett. 17, 891–901. doi: 10.1111/ele.12292

PubMed Abstract | CrossRef Full Text | Google Scholar

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. Nat. Commun. 8, 1–8. doi: 10.1038/ncomms15211

PubMed Abstract | CrossRef Full Text | Google Scholar

Yachi, S., and Loreau, M. (1999). Biodiversity and ecosystem productivity in a fluctuating environment: the insurance hypothesis. Proc. Natl. Acad. Sci. U.S.A. 96, 1463–1468. doi: 10.1073/pnas.96.4.1463

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Q., Fowler, M. S., Jackson, A. L., and Donohue, I. (2019). The predictability of ecological stability in a noisy world. Nat. Ecol. Evol. 3, 251–259. doi: 10.1038/s41559-018-0794-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Yodzis, P. (1981). The stability of real ecosystems. Nature 289, 674–676. doi: 10.1038/289674a0

CrossRef Full Text | Google Scholar