Abstract
There are many scales at which to quantify stability in spatial and ecological networks. Local-scale analyses focus on specific nodes of the spatial network, while regional-scale analyses consider the whole network. Similarly, species- and community-level analyses either account for single species or for the whole community. Furthermore, stability itself can be defined in multiple ways, including resistance (the inverse of the relative displacement caused by a perturbation), initial resilience (the rate of return after a perturbation), and invariability (the inverse of the relative amplitude of the population fluctuations). Here, we analyze the scale-dependence of these stability properties. More specifically, we ask how spatial scale (local vs. regional) and ecological scale (species vs. community) influence these stability properties. We find that regional initial resilience is the weighted arithmetic mean of the local initial resiliences. The regional resistance is the harmonic mean of local resistances, which makes regional resistance particularly vulnerable to nodes with low stability, unlike regional initial resilience. Analogous results hold for the relationship between community- and species-level initial resilience and resistance. Both resistance and initial resilience are “scale-free” properties: regional and community values are simply the biomass-weighted means of the local and species values, respectively. Thus, one can easily estimate both stability metrics of whole networks from partial sampling. In contrast, invariability generally is greater at the regional and community-level than at the local and species-level, respectively. Hence, estimating the invariability of spatial or ecological networks from measurements at the local or species level is more complicated, requiring an unbiased estimate of the network (i.e., region or community) size. In conclusion, we find that scaling of stability depends on the metric considered, and we present a reliable framework to estimate these metrics.
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 (; ; ). 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 (; Wang et al., 2017; ; ; ; Figure 1B). Thus, the ecological and spatial scale at which one studies an ecological system can be hypothesized to influence stability assessments.
FIGURE 1
Previous studies have found that species diversity increases the invariability of communities (; ; ), usually as a consequence of the asynchrony of the population dynamics (Yachi and Loreau, 1999; ): variability decreases with species diversity because of the statistical averaging of the fluctuations in species’ abundances for not perfectly synchronous dynamics (). 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 (). This result can create differences between species- and community-scale stability assessments (; ; ).
Analogously, studies of meta-communities (defined as sets of local communities that are linked by dispersal of multiple potentially interacting species; ) have also proven that temporal invariability increases from local species to meta-communities (Wang and Loreau, 2014), again as a consequence of decreasing synchrony (). 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: ; ; , ). Heterogeneous environmental conditions () and the dispersal ability of the species () 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 (; ; ; ; ).
Finally, ecological and spatial scales may interact. For instance, spatial scale affects stability more in communities than in populations (), depending on the position of the focal species within the community (). Furthermore, decreasing the size of spatial networks reduces species richness more than one would expect from spatial samplings ().
To understand the scaling of stability in meta-communities, which will allow the comparison of stability of systems analyzed at different scales, 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 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 (). 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) (). 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 t0,
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 (; ),
Instead of working with this measure, sometimes its inverse Ω−1 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(t0 + δt)/N(t0)) (), whose absolute value will also decrease as systems become more resistant. Actually, in absolute value this logarithmic definition is at first order equivalent to Ω−1 for small perturbations (as can be proven by applying a Taylor expansion on |N(t0 + δt)−N(t0)|/N(t0)≪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
where N* stands for the equilibrium biomass, which is assumed to be equal to the biomass just before the perturbation (N(t0)). This definition stands for the short term recovery rate after a perturbation (), and has been sometimes referred to as reactivity (). This initial resilience ρ can be expressed in terms of the already considered growth rate R (Eq. 1) and resistance Ω (Eq. 2),
Invariability
We define invariability I as the ratio of the square temporal mean of the biomass and its temporal variance ():
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 (; ), 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, Rx,i, as the normalized time derivative of the species local biomass, Nx,i,
Defining the regional biomass of one species as the sum of all local biomasses of that species across the spatial network, Ni ≡ ∑xNx,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
meaning that the regional growth rate of the species is the weighted arithmetic mean of local species growth rates, Rx,i, with weights equal to the local species biomasses at the moment of the perturbation, Nx,i(t0) (Figure 2). Analogously, the local community growth rate , or the growth rate of the sum of biomasses across all the species of the community at a specific location Nx ≡ ∑iNx,i, is
i.e., the local community growth rate is the weighed arithmetic mean of local species growth rates of each of the species, Rx,i, with weights equals to the local species biomasses at the moment of the perturbation, Nx,i(t0) (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, NT ≡ ∑x ∑iNx,i, which is given by
FIGURE 2
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 Rnet 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). Rnet 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 Rnet can be also expressed as
where μN ≡ mean(N) and μR ≡ mean(R) denote unweighted population arithmetic means of biomasses and growth rates, respectively, and are their standard deviations computed as the square root of the variances, and cN,R the normalized correlation between biomass and growth rate (see Table 1 for all the mathematical definitions). Given that −1 ≤ cN,R ≤ 1, the network growth rate Rnet 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
| SYMBOL | DESCRIPTION |
| N | Biomass. The equilibrium biomass is denoted as N*. Moreover, the biomass of species i at location x and time t is denoted as Nx,i(t). |
| R | Growth rate of the biomass after a perturbation, . Measures of growth rate at regional (R(ℛ)), community (R(𝒞)), and regional community (R(ℛ𝒞)) scales can be computed or estimated with weighted arithmetic means, Eqs. (7–9), or with Eq. (10). The error of the estimate is given by Eq. (11). |
| Ω | Resistance of the biomass to a perturbation, Ω≡N(t0)/(N(t0)−N(t0+δt)). Measures of resistance at regional (Ω(ℛ)), community (Ω(𝒞)), and regional community (Ω(ℛ𝒞)) scales can be computed or estimated with harmonic means (Eq. 13), or with Eq. (14). The error or the estimate is given in Eq. (A5) in Supplementary Material. |
| ρ | Initial resilience of the biomass after a perturbation, . Measures of initial resilience at regional (ρ(ℛ)), community (ρ(𝒞)), and regional community (ρ(ℛ𝒞)) scales can be computed or estimated with the estimates of growth rate R and resistance Ω (Eq. 16). The error or the estimate can be also determined with the estimates of R and Ω (Supplementary Appendix A.3). |
| I | Invariability of the biomass, defined as the inverse of the squared temporal coefficient of variation of the biomass, . Measures of invariability at regional (I(ℛ)), community (I(𝒞)), and regional community (I(ℛ𝒞)) scales can be computed or estimated with Eq. (22) (and Eq. D14 of the Supplementary Material), and generally depends on the number of locations and species. The error of the estimate should be determined with bootstrapping techniques. |
| μXandMX | Population and sample unweighted means of variable X. The population mean is computed when all n nodes of the network were measured, , while the sample mean is computed when just a sample of nodes are measured, . |
| Population and sample variances of variable X: and . | |
| cX.YandCX,Y | Population and sample Pearson correlation coefficients of variables X and Y: and . Even though they are usually denoted as ρX,Y and rX,Y, the employed notation was preferred to avoid possible confusions with initial resilience. |
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, . We can then estimate the network growth rate using Eq. (10) with the sampled nodes, replacing the population means (μN and μR) by the sample means (MN and MR), the population standard deviations (σN and σR) by the sample standard deviations (SN and SR), and the population correlation (cN,R) by the sample correlation (CN,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 and validated by (see Supplementary Appendix A)
where is the standard error of the network, corresponding to a 95% confidence level; is the sample size (i.e., the number of sampled nodes); and tñ–1 is the Student’s t distribution with 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 , as expected, decreases when sampling more localities or species (larger ), and when the growth rates vary less across localities or species (smaller SR). Note that, generally, SN≤MN, and so also biomasses that vary less relative to their average biomass will lead to less variable estimates of Rnet. Another implication is that positive or negative correlations between biomass and growth rate will mostly decrease , because SN/MN≤1 and so (SN/MN)4≪(SN/MN)2. 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
Resistance
The resistance of species i at location x, defined as in Eq. (2), is
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),
FIGURE 4
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
where again μ represents unweighted arithmetic means, σ the population standard deviations, and cN,Ω–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 (Ω−1)net is the weighted arithmetic mean of the sub-units estimates of Ω−1. 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
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
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
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Ω−1 (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:
FIGURE 5
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 . And assuming no correlations between resistance and growth rates, the standard error committed with such sampling would be . 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
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)
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
For the asynchronous-space case, the regional invariability is proportional to the number of locations nL (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 and mean 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
The more general case of not perfectly synchronous dynamics can be expressed as
where is the typical correlation between different locations (Eq. D10 of the Supplementary Appendix D). For the case , and since by definition , would simply be the weighted harmonic mean of and , with weights equal to and , respectively. And since, even though does not depend on the number of locations nL (Eq. 20), increases with nL (Eq. 21), the resulting regional species invariability would increase as well with nL (except for the special case ). For the case , since , the regional species invariability would be larger than . Since 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 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 (; ).
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 (; ; ; ; ), we complement these results with results for realistic riverine dendritic networks (Supplementary Figure 1B) generated in R () with the OCNet package (; Supplementary Figures 2–5). All simulations of community dynamics were done in Python 3.7 ().
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 (; ).
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 (; ; ) and spatial scale (Wang and Loreau, 2014; ), as a consequence of not perfectly synchronous dynamics among species and among locations (). 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 (). To solve this issue, 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 (), network resistance and short-term or initial growth rate and resilience [independent of asymptotic resilience, and a better proxy of resilience in experiments (, )] 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 (), depending on the characteristic of the environmental fluctuations affecting the systems (; ), 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 ().
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 (). 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 (). 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 () is expected to have a stronger effect on the network resistance than in the network initial resilience.
FIGURE 7
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 (; ). 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). 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 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 on the scale dependence of resistance can be also explained from biomass normalization (Supplementary Appendix G).
For initial resilience, 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). 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., Ω−1. More resistant systems present smaller values of Ω−1, and Ω−1 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 Ω−1 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 Ω−1 can be easily compensated by nodes with below-average values of Ω−1, having a limited effect on the network-level estimate of Ω−1. This is a clear indication of how the heterogeneous distribution of local species estimates (particularly its skewness ) 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 , and the change in biomass due to a perturbation N(t0)−N(t0δ + 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 Ω−1 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 Ω−1 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 vart(N(t)) is not extensive, because vart(N) = Et[(N−Et[N])2] = Et[N2]−(Et[N])2 is not additive for subsystems. Neither 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 (; ; ) 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 (; ) originated from the disagreement between experimental observations often finding more complex systems to be more stable (), and theoretical analyses finding more complex systems to be less stable (; ). To solve this disagreement, some authors have proposed generalizations of the original work of that account for non-random among-species interactions (Yodzis, 1981; ), or the stabilizing role of dispersal and spatial heterogeneity (; ). Other approaches have suggested that the disagreement is caused by a focus on asymptotic resilience in theoretical studies (; ; ), which is biased by rare species (). 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 (; ). Assessing stability from a multi-dimensional perspective (; ; ) will provide a more comprehensive picture and can reconcile apparent contradictions between and among theoretical and empirical studies.
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.
Statements
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.
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.
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.
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
1
AllesinaS.TangS. (2015). The stability–complexity relationship at age 40: a random matrix perspective.Popul. Ecol.5763–75. 10.1007/s10144-014-0471-0
2
AltermattF. (2013). Diversity in riverine metacommunities: a network perspective.Aquat. Ecol.47365–377. 10.1007/s10452-013-9450-3
3
AmarasekareP. (2008). Spatial dynamics of foodwebs.Annu. Rev. Ecol. Evol. Syst.39479–500. 10.1146/annurev.ecolsys.39.110707.173434
4
ArnoldiJ.-F.BideaultA.LoreauM.HaegemanB. (2018). How ecosystems recover from pulse perturbations: a theory of short- to long-term responses.J. Theor. Biol.43679–92. 10.1016/j.jtbi.2017.10.003
5
ArnoldiJ. F.LoreauM.HaegemanB. (2016). Resilience, reactivity and variability: a mathematical comparison of ecological stability measures.J. Theor. Biol.38947–59. 10.1016/j.jtbi.2015.10.012
6
ArnoldiJ. F.LoreauM.HaegemanB. (2019). The inherent multidimensionality of temporal variability: how common and rare species shape stability patterns.Ecol. Lett.221557–1567. 10.1111/ele.13345
7
BaertJ. M.De LaenderF.SabbeK.JanssenC. R. (2016). Biodiversity increases functional and compositional resistance, but decreases resilience in phytoplankton communities.Ecology973433–3440. 10.1002/ecy.1601
8
CarraraF.AltermattF.Rodriguez-IturbeI.RinaldoA. (2012). Dendritic connectivity controls biodiversity patterns in experimental metacommunities.Proc. Natl. Acad. Sci. U.S.A.1095761–5766. 10.1073/pnas.1119651109
9
CarraroL.BertuzzoE.FronhoferE. A.FurrerR.GounandI.RinaldoA.et al (2020). Generation and application of river network analogues for use in ecology and evolution.Ecol. Evol.107537–7550. 10.1002/ece3.6479
10
ChaseJ. M.BlowesS. A.KnightT. M.GerstnerK.MayF. (2020). Ecosystem decay exacerbates biodiversity loss with habitat loss.Nature584238–243. 10.1038/s41586-020-2531-2
11
ChessonP. (2000). General theory of competitive coexistence in spatially-varying environments.Theor. Popul. Biol.58211–237. 10.1006/tpbi.2000.1486
12
ClarkA. T.ArnoldiJ. F.ZelnikY. R.BarabasG.HodappD.KarakoçC.et al (2021). General statistical scaling laws for stability in ecological systems.Ecol. Lett.241474–1486. 10.1111/ele.13760
13
CochranW. G. (1977). Sampling Techniques.3rd Edn. New York, NY: John Wiley & Sons.
14
De RaedtJ.BaertJ. M.JanssenC. R.De LaenderF. (2019). Stressor fluxes alter the relationship between beta-diversity and regional productivity.Oikos1281015–1026. 10.1111/oik.05191
15
DoakD. F.BiggerD.HardingE. K.MarvierM. A.O’MalleyR. E.ThomsonD. (1998). The statistical inevitability of stability-diversity relationships in community ecology.Am. Nat.151264–276. 10.1086/286117
16
Domínguez-GarcíaV.DakosV.KéfiS. (2019). Unveiling dimensions of stability in complex ecological networks.Proc. Natl. Acad. Sci. U.S.A.11625714–25720. 10.1073/pnas.1904470116
17
DonohueI.PetcheyO. L.MontoyaJ. M.JacksonA. L.McnallyL.VianaM.et al (2013). On the dimensionality of ecological stability.Ecol. Lett.16421–429. 10.1111/ele.12086
18
DowningA. L.BrownB. L.LeiboldM. A. (2014). Multiple diversity-stability mechanisms enhance population and community stability in aquatic food webs.Ecology95173–184. 10.1890/12-1406.1
19
EfronB.TibshiraniR. (1985). The bootstrap method for assessing statistical accuracy.Behaviormetrika121–35. 10.2333/bhmk.12.17_1
20
FaganW. F. (2002). Connectivity, fragmentation, and extinction risk in dendritic metapopulations.Ecology833243–3249.
21
FergerW. F. (1931). The Nature and Use of the Harmonic Mean.J. Am. Stat. Assoc.2636–40. 10.1080/01621459.1931.10503148
22
FlöderS.HillebrandH. (2012). Species traits and species diversity affect community stability in a multiple stressor framework.Aquat. Biol.17197–209. 10.3354/ab00479
23
GatzD. F.SmithL. (1995). The standard error of a weighted mean concentration-I. Bootstrapping vs other methods.Atmos. Environ.291185–1193. 10.1016/1352-2310(94)00210-C
24
GravelD.MassolF.LeiboldM. A. (2016). Stability and complexity in model meta-ecosystems.Nat. Commun.7:12457. 10.1038/ncomms12457
25
GreigH. S.McHughP. A.ThompsonR. M.WarburtonH. J.McIntoshA. R. (2022). Habitat size influences community stability.Ecology103:e03545. 10.1002/ecy.3545
26
GrimmV.WisselC. (1997). Babel, or the ecological stability discussions: an inventory and analysis of terminology and a guide for avoiding confusion.Oecologia109323–334. 10.1007/s004420050090
27
GrossK.CardinaleB. J.FoxJ. W.GonzalezA.LoreauM.Wayne PolleyH.et al (2014). Species richness and the temporal stability of biomass production: a new analysis of recent biodiversity experiments.Am. Nat.1831–12. 10.1086/673915
28
HaegemanB.ArnoldiJ.-F.WangS.de MazancourtC.MontoyaJ.LoreauM. (2016). Resilience, invariability, and ecological stability across levels of organization.bioRxiv[Preprint].10.1101/085852
29
HesterbergT. (2011). Bootstrap.Wiley Interdiscip. Rev. Comput. Stat.3497–526. 10.1002/wics.182
30
HillebrandH.LangenhederS.LebretK.LindströmE.ÖstmanÖStriebelM. (2018). Decomposing multiple dimensions of stability in global change experiments.Ecol. Lett.2121–30. 10.1111/ele.12867
31
IsbellF.CravenD.ConnollyJ.LoreauM.SchmidB.BeierkuhnleinC.et al (2015). Biodiversity increases the resistance of ecosystem productivity to climate extremes.Nature526574–577. 10.1038/nature15374
32
IUPAC (2019). The IUPAC Compendium of Chemical Terminology – The Gold Book, 2nd Edn, ed.GoldV. (Research Triangle Park, NC: International Union of Pure and Applied Chemistry (IUPAC)). 10.1351/goldbook
33
IvesA. R.CarpenterS. R. (2007). Stability and diversity of ecosystems.Science31758–62. 10.1126/science.1133258
34
IvesA. R.KlugJ. L.GrossK. (2000). Stability and species richness in complex communities.Ecol. Lett.3399–411. 10.1046/j.1461-0248.2000.00144.x
35
JarilloJ.SætherB.-E.EngenS.CaoF. J. (2018). Spatial scales of population synchrony of two competing species: effects of harvesting and strength of competition.Oikos1271459–1470. 10.1111/oik.05069
36
JarilloJ.SætherB.-E.EngenS.Cao-GarcíaF. J. (2020). Spatial scales of population synchrony in predator-prey systems.Am. Nat.195216–230. 10.1086/706913
37
KarakoçC.ClarkA. T.ChatzinotasA. (2020). Diversity and coexistence are influenced by time-dependent species interactions in a predator–prey system.Ecol. Lett.23983–993. 10.1111/ele.13500
38
KéfiS.Domínguez-GarcíaV.DonohueI.FontaineC.ThébaultE.DakosV. (2019). Advancing our understanding of ecological stability.Ecol. Lett.221349–1356. 10.1111/ele.13340
39
LandeR.EngenS.SætherB.-E. (1999). Spatial scale of population synchrony: environmental correlation versus dispersal and density regulation.Am. Nat.154271–281. 10.1086/303240
40
LeiboldM. A.HolyoakM.MouquetN.AmarasekareP.ChaseJ. M.HoopesM. F.et al (2004). The metacommunity concept: a framework for multi-scale community ecology.Ecol. Lett.7601–613. 10.1111/j.1461-0248.2004.00608.x
41
LemoineN. P. (2020). Unifying ecosystem responses to disturbance into a single statistical framework.Oikos130, 1–14. 10.1111/oik.07752
42
LevinS. A. (1992). The problem of pattern and scale in ecology: the Robert H. Macarthur award lecture.Ecology731943–1967. 10.2307/1941447
43
LimbergerR.PittA.HahnM. W.WickhamS. A. (2019). Spatial insurance in multi-trophic metacommunities.Ecol. Lett.221828–1837. 10.1111/ele.13365
44
LiuJ.SoininenJ.HanB. P.DeclerckS. A. J. (2013). Effects of connectivity, dispersal directionality and functional traits on the metacommunity structure of river benthic diatoms.J. Biogeogr.402238–2248. 10.1111/jbi.12160
45
LoreauM.De MazancourtC. (2008). Species synchrony and its drivers: neutral and nonneutral community dynamics in fluctuating environments.Am. Nat.17248–66. 10.1086/589746
46
MayR. M. (1972). Will a large complex system be stable?Nature238413–414. 10.1038/238413a0
47
McCannK. S. (2000). The diversity-stability debate.Nature405228–233. 10.1038/35012234
48
McCluneyK. E.PoffN. L.PalmerM. A.ThorpJ. H.PooleG. C.WilliamsB. S.et al (2014). Riverine macrosystems ecology: sensitivity, resistance, and resilience of whole river basins with human alterations.Front. Ecol. Environ.12:48–58. 10.1890/120367
49
MoranP. A. P. (1953). The statistical analysis of the Canadian lynx cycle. II. Synchronization and meteorology.Aust. J. Zool.1291–298. 10.1071/ZO9530291
50
MougiA.KondohM. (2012). Diversity of interaction types and ecological community stability.Science337349–351. 10.1126/science.1220529
51
MougiA.KondohM. (2016). Food-web complexity, meta-community complexity and community stability.Sci. Rep.61–5. 10.1038/srep24478
52
NeubertM. G.CaswellH. (1997). Alternatives to resilience for measuring the responses of ecological systems to perturbations.Ecology78653–665. 10.1111/gcb.12845
53
PennekampF.PontarpM.TabiA.AltermattF.AltherR.ChoffatY.et al (2018). Biodiversity increases and decreases ecosystem stability.Nature563109–112. 10.1038/s41586-018-0627-8
54
PetersonE. E.Ver HoefJ. M.IsaakD. J.FalkeJ. A.FortinM. J.JordanC. E.et al (2013). Modelling dendritic ecological networks in space: an integrated network perspective.Ecol. Lett.16707–719. 10.1111/ele.12084
55
PimmS. L. (1984). The complexity and stability of ecosystems.Nature307321–326. 10.1038/307321a0
56
PlitzkoS. J.DrosselB. (2015). The effect of dispersal between patches on the stability of large trophic food webs.Theor. Ecol.8233–244. 10.1007/s12080-014-0247-3
57
Python Core Team (2019). Python: A Dynamic, Open Source Programming Language. Available online at: https://www.python.org/(accessed December 2021).
58
R Core Team (2020). R: A Language and Environment for Statistical Computing.Vienna: R foundation for statistical computing.
59
RadchukV.LaenderF.De, CabralJ. S.BoulangeatI.CrawfordM.et al (2019). The dimensionality of stability depends on disturbance type.Ecol. Lett.22674–684. 10.1111/ele.13226
60
RooneyN.McCannK.GellnerG.MooreJ. C. (2006). Structural asymmetry and the stability of diverse food webs.Nature442265–269. 10.1038/nature04887
61
SaadeC.KéfiS.Gougat-BarberaC.RosenbaumB.FronhoferE. A. (2020). Spatial distribution of local patch extinctions drives recovery dynamics in metacommunities.bioRxiv[Preprint].10.1098/rspb.2022.0543.2020.12.03.409524
62
SaeedianM.PiganiE.MaritanA.SuweisS.AzaeleS. (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). 10.1098/rsta.2021.0245
63
StevensS. S. (1955). On the averaging of data.Science121113–116. 10.1126/science.121.3135.113
64
SuppS. R.ErnestS. K. M. (2014). Species-level and community-level responses to disturbance: a cross-community analysis.Ecology951717–1723. 10.1890/13-2250.1
65
ThébaultE.LoreauM. (2005). Trophic interactions and the relationship between species diversity and ecosystem stability.Am. Nat.166E95–E114. 10.1086/444403
66
ThibautL. M.ConnollyS. R. (2013). Understanding diversity-stability relationships: towards a unified model of portfolio effects.Ecol. Lett.16140–150. 10.1111/ele.12019
67
TilmanD.ReichP. B.KnopsJ. M. H. (2006). Biodiversity and ecosystem stability in a decade-long grassland experiment.Nature441629–632. 10.1038/nature04742
68
WangS.LamyT.HallettL. M.LoreauM. (2019). Stability and synchrony across ecological hierarchies in heterogeneous metacommunities: linking theory to data.Ecography421200–1211. 10.1111/ecog.04290
69
WangS.LoreauM. (2014). Ecosystem stability in space: α, β and γ variability.Ecol. Lett.17891–901. 10.1111/ele.12292
70
WangS.LoreauM.ArnoldiJ. F.FangJ.RahmanK. A.TaoS.et al (2017). An invariability-area relationship sheds new light on the spatial scaling of ecological stability.Nat. Commun.81–8. 10.1038/ncomms15211
71
YachiS.LoreauM. (1999). Biodiversity and ecosystem productivity in a fluctuating environment: the insurance hypothesis.Proc. Natl. Acad. Sci. U.S.A.961463–1468. 10.1073/pnas.96.4.1463
72
YangQ.FowlerM. S.JacksonA. L.DonohueI. (2019). The predictability of ecological stability in a noisy world.Nat. Ecol. Evol.3251–259. 10.1038/s41559-018-0794-x
73
YodzisP. (1981). The stability of real ecosystems.Nature289674–676. 10.1038/289674a0
Summary
Keywords
scale, stability, resistance, invariability, regional, community, initial resilience
Citation
Jarillo J, Cao-García FJ and De Laender F (2022) Spatial and Ecological Scaling of Stability in Spatial Community Networks. Front. Ecol. Evol. 10:861537. doi: 10.3389/fevo.2022.861537
Received
24 January 2022
Accepted
09 June 2022
Published
30 June 2022
Volume
10 - 2022
Edited by
Adam Clark, University of Graz, Austria
Reviewed by
Azenor Bideault, Université de Sherbrooke, Canada; Pierre Quévreux, UMR 5321 Station d’Ecologie Théorique et Expérimentale (SETE), France
Updates
Copyright
© 2022 Jarillo, Cao-García and De Laender.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Javier Jarillo, jjarillo@ucm.es
†Present address: Javier Jarillo, Departamento de Estadística e Investigación Operativa, Universidad Complutense de Madrid, Madrid, Spain
This article was submitted to Models in Ecology and Evolution, a section of the journal Frontiers in Ecology and Evolution
Disclaimer
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.