Spatial and Ecological Scaling of Stability in Spatial Community Networks

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.


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) (Fig. 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) (Fig. 1B). Thus, the ecological and spatial scale at which one studies an ecological system can be hypothesized to influence stability assessments.

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 metacommunities, 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., 2018Jarillo et al., , 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 (Fig. 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 2) by introducing biomass-normalized stability measures of growth rate, resistance, initial resilience, and invariability. In Section 3, 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 4 shows that our formulas compare well with model simulations. Finally, Section 5 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 (Fig. 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 as the relative instantaneous return rate to the equilibrium of the biomass after any sudden biomass change caused by any external perturbation at time 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 2.3) 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), 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� ( 0 + )/ ( 0 )� (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 Ω −1 for small perturbations (as can be proven by applying a Taylor expansion on , 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 * stands for the equilibrium biomass, which is assumed to be equal to the biomass just before the perturbation ( ( 0 )). 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 (Eq. (1)) and resistance Ω (Eq. (2)),

Invariability
We define invariability as the ratio of the square temporal mean of the biomass and its temporal variance (Thibaut and Connolly, 2013): 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 (Fig. 1B) of the previously introduced stability measures (Section 2, Fig. 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 at a given specific location , , , as the normalized time derivative of the species local biomass, , , Defining the regional biomass of one species as the sum of all local biomasses of that species across the spatial network, ≡ ∑ , , 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 of the species is meaning that the regional growth rate of the species is the weighted arithmetic mean of local species growth rates, , , with weights equal to the local species biomasses at the moment of the perturbation, , ( 0 ) (Fig. 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 ≡ ∑ , , is i.e., the local community growth rate is the weighed arithmetic mean of local species growth rates of each of the species, , , with weights equals to the local species biomasses at the moment of the perturbation, , ( 0 ) (Fig. 2). Finally, we can define regional community growth rate (ℛ ) (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, ≡ ∑ ∑ , , which is given by Thus, the regional growth rate of a community after a perturbation, (ℛ ) , 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 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).
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 Eq. (7), (8) or (9) respectively. This can be also expressed as where ≡ mean( ) and ≡ mean( ) denote unweighted population arithmetic means of biomasses and growth rates, respectively, ≡ �var( ) and ≡ �var( ) are their standard deviations computed as the square root of the variances, and , the normalized correlation between biomass and growth rate. (See Table 1 for all the mathematical definitions). Given that −1 ≤ , ≤ 1, the network growth rate can be greater or smaller than the unweighted mean of growth rates depending on the positive or negative correlations between the node biomasses and node growth rates.
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 (panel A) and community growth rate (panel 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 6300 random communities of 10 competitors in 10-node random spatial networks ( Fig. S1A and text on Supplementary Material), 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.

SYMBOL DESCRIPTION
Biomass. The equilibrium biomass is denoted as * . Moreover, the biomass of species at location and time is denoted as , ( ).
Ω Resistance of the biomass to a perturbation, 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 the 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 and resistance Ω (Eq. (16)). The error or the estimate can be also determined with the estimates of and Ω (see Appendix A3 of the Supplementary Material).
. and , Population and sample Pearson correlation coefficients of variables and : . ≡ ). Even though they are usually denoted as , and , , the employed notation was preferred to avoid possible confusions with initial resilience. 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) (Fig. 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 ( and ) by the sample means ( and ), the population standard deviations ( and ) by the sample standard deviations ( and ), and the population correlation ( , ) by the sample correlation ( , ) (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 where SE� net � � 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 �−1 is the Student's t distribution with � − 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 SE� net � �, as expected, decreases when sampling more localities or species (larger �), and when the growth rates vary less across localities or species (smaller ). Note that, generally, < , and so also biomasses that vary less relative to their average biomass will lead to less variable estimates of . Another implication is that positive or negative correlations between biomass and growth rate will mostly decrease SE� net � �, because / < 1 and so ( / ) 4 ≪ ( / ) 2 . Nevertheless, since negative correlations decreases the network growth rate (Eq. (10)), negative correlations between and require larger sample sizes to control the relative error on the estimated network growth rate (Fig. 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 ( / ), and to a lesser extent with the coefficient of variation of the node biomass ( / ). Larger sampling sizes are required for more negative values of the cross-correlation between the biomass and the growth rate , (Table 1).

Resistance
The resistance of species at location , 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 (Appendix B, Fig. 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 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 ,Ω −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.
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 (panel A) and community resistance (panel B). Black dots are the resistances of the individual localities (A) or species (B) on which the means were computed. Data were generated for 6300 random communities of 10 competitors in 10-node random spatial networks ( Fig. S1A and text on Supplementary Material), 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 Appendix C of the Supplementary Material). Moreover, HM w are found to be better estimators for the resistance, as expected from the results shown in the text.
Even though the resistance of a network Ω is the weighted harmonic mean of resistances, the network estimate of its inverse (Ω −1 ) 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 (Appendix A of the Supplementary Material), which let us to obtain an expression for the standard error of the resistance obtained from a partial sampling of a network (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 , i.e., = Ω . Therefore, to estimate for different scales one can use the already obtained scaling of Ω and (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 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, Ω −1 (Fig. 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: 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.
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 (panel A) and community initial resilience (panel 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 6300 random communities of 10 competitors in 10-node random spatial networks ( Fig. S1A and text on Supplementary Material), 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.
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 SE( � ) = Ω � SE� � � + � SE�Ω � �. 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 (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 reads (see 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 (Fig. 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 * 2 and mean * of the local equilibrium biomasses of species . It can be proven that invariability is higher for asynchronous than for synchronous dynamics (see 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 (Fig. 6B) 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 panels A and B, we have considered three different scenarios: asynchronous local dynamics ( ̅ = 0, yellow box plots), partiallysynchronous local dynamics ( ̅ = 0.36, blue), and perfectly-synchronous local dynamics ( ̅ = 1, green). In panels C and D, we have considered the cases of asynchronous (̃= 0), partiallysynchronous (̃= 0.17) and perfectly-synchronous (̃= 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 where ̅ is the typical correlation between different locations (Eq. (D10) of the appendix D). For the case ̅ > 0, and since by definition ̅ ≤ 1, (ℛ) would simply be the weighted harmonic mean of (ℛ; ) and (ℛ; ) , with weights equal to 1 − ̅ and ̅ , respectively. And since, even though does not depend on the number of locations (Eq. (20)), (ℛ; ) increases with (Eq. (21)), the resulting regional species invariability would increase as well with (except for the special case ̅ = 1). For the case ̅ < 0, 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 (Fig. 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 (Fig. 6C), while the degree of synchrony and the number of species do not significantly affect the invariabilities at the species level (Fig. 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 (Appendix A of the Supplementary Material). 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 Appendix E of the Supplementary Material) in 10-node random spatial networks (Fig. S1A of the Supplementary Material) (Figs. 2,(4)(5)(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 ( 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 (Figs. 2 and 5, Figs. S2 and S4 of the Supplementary Material). 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 (Fig. 4, Fig. S3 of the Supplementary Material). 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 (Fig. 6, Fig. S5 of the Supplementary Material), 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 (Fig. 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(Arnoldi et al., , 2018) is the mean of the local and species values (Figs. 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 (Figs. S2-S5 of the Supplementary Material).
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 Appendix C of the Supplementary Material). 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 (Fig.  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.
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 (Fig. 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 results (Appendix F of the Supplementary Material). 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 (Appendix G of the Supplementary Material).
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 (Fig. 5). Clark et al. 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 5.2, 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 (Stevens, 1955)) affects the network resistance and initial resilience.
The scale-free property found for the growth rate , the resistance Ω, and the initial resilience is due to their character as intensive quantities. The total biomass , its derivative d d , and the change in biomass due to a perturbation ( 0 ) − ( 0 + ), 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 , resistance Ω, and initial resilience . For the simpler cases, the growth rate 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 var � ( )� is not extensive, because var ( ) = [( − [ ]) 2 ] = [ 2 ] − ( [ ]) 2 is not additive for subsystems. Neither �var ( ) 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 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 Fig. 3, and Fig. S6 of the Supplementary Material). The coefficient of variation of the studied stability property (resistance or initial resilience) affects this required sampling effort most (Fig. 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 Fig. S6 of the Supplementary Material.) 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., 2021), 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.

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.

Author Contributions
JJ and FDL conceived the presented idea. JJ performed the numerical simulations and the analytical computations. FJCG and FDL verified the analytical derivations. JJ wrote the first draft and prepared the figures. All authors discussed the results and contributed to the final manuscript.  Section 3 of the main text shows that local population biomasses and local population-level growth rates suffice to estimate the growth rate of the spatial networks of communities (Eqs. (7-9)). The network-level growth rate is the arithmetic mean of node-level growth, weighted by the node biomasses. However, in real networks, measuring the growth rate at each of the nodes of the network might be not feasible. The level of uncertainty that we will have in our regional or community estimate � coming from an incomplete number of measurements will be given by the standard error on this estimate. The standard error of a weighted arithmetic mean (here defined as half the width of the 95% confidence interval of the weighted mean from a sampling of size �) follows the expression obtained by , and validated with bootstrapping techniques , where is the unweighted arithmetic sample mean of the biomasses, ≡ 1 � ∑ �

=1
(see Table 1 of the main text); � ≡ ∑ � =1 / ∑ =1 is the estimate of network growth rate, obtained as the weighted arithmetic mean of growth rates of the sampled species or locations, with weights equal to the local or species biomasses ; and � is the number of sampled locations or species. �−1 is the Student t-distribution with � − 1 degrees of freedom corresponding to the 95% confidence interval. Its value for large enough sampling sizes ( � → ∞) is approximately 1.96. Using the general expression for the covariance of the product of random variables (Bohrnstedt and Goldberger, 1969), we can write the covariance of with the product ⋅ as cov( , = cov( , ) + var( ) (where "cov" denotes the sampling covariance , cov( , ) ≡ ⋅ − , and "var" the sampling variance, var( ) ≡ 2 ≡ 2 − 2 ). Then, this general expression can be further reduced to where = �var( ) is the standard deviation of the sampled values of growth rate at the nodes of the network, / is the coefficient of variation of the node biomasses (the ratio of the sample standard deviation and the sample mean of the node biomasses), and , ≡ cov( , )/( ) is the correlation of the nodes biomasses and growth rates (see Table 1 of the main text). Note that for the special case of no variation in the local or species biomasses, the standard error of the unweighted arithmetic mean is recovered, SE( ) = �−1 √ �−1 . The standard error estimated with Eq. (A2) equals the standard error that can be obtained numerically with bootstrapping techniques , although in our simulations it produces a slight overestimation (Figs. S7-S8). It is important to note that this analytical expression is valid for any possible community or spatial network, regardless of its complexity. We can further normalize by the network growth rate to obtain an expression for the relative standard error Using Eq. (A3), we can estimate the required sample size needed to control the standard error of the growth rate estimate of a sampled network based on the coefficients of variation and on the correlation of the node biomasses and growth rates (Fig. 3 of the main text, Fig. S6). This required sample size is mainly determined by the coefficient of variation of the node growth rates (Fig. 3), having the coefficient of variation of the node biomasses and the correlation between the biomass and the growth rate smaller effects except for the case of high anticorrelation (Fig. S6).

A.2. Resistance
Network resistance has been obtained as the weighted harmonic mean of local or species resistances (Eq. (12) of the main text, and Appendix B). Hence, it is easy to check that the network estimate of the inverse of resistance, Ω −1 , would simply be given as the weighted arithmetic mean of the individual values. The standard error of Ω −1 would then be equivalent to that expressed in Eq. (A2), changing by Ω −1 , Moreover, the standard error of the unweighted harmonic mean of any random variable is just equal to the standard error of the arithmetic mean of −1 , multiplied by the squared harmonic mean of (Norris, 1940). Assuming this approximately holds also for the weighted case, we would obtain that the standard error of the network resistance is the harmonic mean of the resistances of the sampled nodes weighted by the biomasses. Again, this analytical expression for the standard error of the resistance produces results compatible with those obtained numerically with a bootstrap (Figs. S7-S8). Moreover, normalizing by the network resistance Ω � , we recover an expression analogous to Eq. (A3), From Eq. (A6) it is possible to estimate the sample size required to control the error of the resistance estimated of a sampled ecological network. Analogously than for the growth rate, the higher element which will mainly determine this sample size is the coefficient of variation of the reciprocal of the node resistance, Ω −1 / Ω −1 , except for the case of high anticorrelation between biomass and the reciprocal of resistance.

A.3. Initial resilience
With respect to the initial resilience, in the main text we have seen that it is just the product of the growth rate and the resistance, = Ω . Hence, we can write the network estimate of initial resilience as the product of the estimates of the growth rates and resistances, � = � Ω net � . And assuming no correlations between the growth rates and the resistances, the standard error of the estimated network initial resilience would be SE( � ) = Ω � SE� � � + � SE�Ω � �. Hence, with the estimations of the network growth rate after a perturbation and resistance we can obtain the estimation (and the standard error) of the network initial resilience.

A.4. Invariability
For the case of invariability, the results are more complicated. For the particular case of perfectly synchronous dynamics, the network invariability is the square of the weighted harmonic mean of the square roots of the invariabilities (Eq. (20) main text). The standard error of the harmonic mean of � would then given by Eq. (A5), replacing Ω by √ , and the biomasses by their steady states * . Regarding the standard error of the square of such mean, the results can be easily computed via error propagation, achieving the estimate However, the more general cases of asynchronous or partially asynchronous dynamics are more involved, since for them the invariability is not a mean of individual invariabilities, but depend on the size of the network. Hence, to assess the uncertainty of the network invariability, we would generally not be able to employ any of these analytical expression, and we should compute it numerically with bootstrapping techniques .
In general, the resistance of any ecological or spatial network can be rewritten as where Ω represents the resistance of any ecological or spatial network, Ω −1 is the unweighted population mean of the node estimates of resistance (so 1/ Ω −1 represents the unweighted population harmonic mean of node resistances); Ω −1 = �var(Ω −1 ) the standard deviation of the reciprocals of the node resistances, and ,Ω −1 = cov( , Ω −1 )/( Ω −1 ) the normalized correlation (see Table 1 of the main text). This equation also indicates that correlations between the node biomasses and the inverse of the resistance can make the network resistance Ω greater or smaller than the unweighted harmonic mean of resistances.
Eq. (B6) can also be applied to estimate the network resistance Ω when only a limited number of nodes have been sampled from the network, so the sample estimate of resistance would be where again represent sample means, sample standard deviations, and sample correlations (see Table 1).
Even though the resistance of a network Ω is the weighted harmonic mean of resistances, the network estimate of its inverse (Ω −1 ) 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 (Appendix A, Eq. (A5)), while the relative standard error follows a formula analogous to that of the growth (Appendix A, Eq. (A6)). Hence, 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.

Appendix C: Relation Between Weighted Harmonic and Arithmetic Mean of Resistances
The regional or community resistance has been obtained to be the harmonic mean of the local or population resistances, weighted by the biomass. (Eq. (13) of the main text). Here, we want to prove that such weighted harmonic mean is upper-bounded by the weighted arithmetic mean of resistances, i.e., We denote ( 0 ) ≡ ( 0 )/ ∑ ( 0 ) the proportion of the biomass of the network at the node . Such node could represent either a location or a species.Prove the relation (C1) is equivalent to prove that Let's start to simplify the expression in Eq. (C2) From the definition of , ∑ = 1 (the sum of proportions over all nodes of the networks is 1). Hence But also All the terms , Ω , and (Ω − Ω ) 2 are non-negative. Hence, the second term in Eq. (C3) is also non-negative, and in consequence where Ψ = 1 only if Ω = Ω for each pair of nodes of the network with non-zero biomass.
Eq. (C4) proves that Eq (C2) is true, so it proves that the weighted harmonic mean of resistances is upper-bounded by the weighted arithmetic mean of resistances (Eq. (C1)). And both means are equal only if there is no variation in the individual resistances of the nodes with non-zero biomass. This result implies that the network resistance estimate, equal to the weighted harmonic mean of the node resistance estimates, is upper bounded by the weighted arithmetic mean of resistances, being more affected by low-resistant nodes than the arithmetic mean is.
For the asynchronous-space case, the regional invariability is proportional to the number of locations (Fig. 6 of the main text, Fig. S5), and to the harmonic mean of local species invariabilities weighted by the squared local species biomasses; and it is modulated by the spatial variance * 2 and mean * 1 10 �. The diagonal terms of the interaction strengths were randomly sampled from a normal distribution with mean 1/100 and a standard deviation of For the dispersal, we generate random networks with different number of nodes, and connections between some of the nodes (Fig. S1). Individuals of each of the population will be able to diffuse through connected nodes of the spatial network, with the expression where the first term represents the emigration from the local node to each of the connected nodes, and the second term represents immigration to from any connected node. The elements → , are equal to zero for any pair of nodes not connected by an edge; while for connected nodes they are sampled from a normal distribution → , ∼ 1 � = 1, = 1 10 �, being . the spatial distance between the nodes.
Finally, for the environmental stochasticity term we have employed where the amplitudes of the environmental stochasticity , have been sampled from an exponential distribution with mean 1 20 . The terms d , d ( ) can be taken as independent white noises for each population and location, for studying the dynamics in the spatial and population asynchronous case; or as a unique common white noise shared by all populations and locations, for studying the dynamics in a perfectly synchronous scenario.
For estimating the invariability of spatial community networks, we run the community dynamics for a time span of 300 units, with time-steps of 0.25, from initial populations equal to the species carrying capacities (in absence of interspecific competition). Then, we select the last half of the time series, in order to remove the transient dynamics. From these truncated time series, we compute the local population temporal means and variances. All the numerical simulation were done in Python 3.7 (Python Core Team, 2019).
For estimating the growth rate, initial resilience, and resistance of spatial community networks, starting again at populations equal to the local carrying capacities, we run the simulations for 100 time units (with time-steps of 0.25) to reach the network equilibrium state. Then, from time 100 we reduce each of the relative local population growth rates a quantity sampled from a normal distribution, , → , �1 − , �, with , ∼ 1 2 � = 1, = 1 4 �, and we run the community dynamics from time 100 to time 200 (to let the network reach again its new equilibrium state), again with time steps equal to 0.25. Resistance is estimated from the relative decrease in the population abundances at times 100 (equilibrium before the perturbation) and times 200 (equilibrium after the perturbation). Finally, at time 200 we set again the local population growth rates to their pre-perturbated values, and we run the community dynamics until time 300 (to check that we return the pre-perturbed equilibrium state). From these last time series, growth rate is estimated as the initial return rate to the pre-perturbed equilibrium in the first time-step.

Appendix F: Scaling of not-normalized invariance
We define the invariance as the inverse of the temporal variance of the population, This invariance can be considered as another component of the ecosystem stability. It is analogous to the invariability (Eq. (5) of the main text), but without normalizing the variance of the abundances with the temporal mean of the biomasses.
If we want to study the invariance of a (either spatial or ecological) network, we define it as the inverse of the variance in the network total biomass, where holds for the node of the network, which can represent either a species or a location, and ℐ holds for the global (either community or regional) invariance. Eq. (F2) can be rewritten as with cov � ( ), ( )� the temporal covariance between the nodes and . (1/ℐ ) where 〈ℐ 〉 is the harmonic mean of the invariances at the nodes of the network.
This not-normalized resistance is equal to the normalized resistance that we employed throughout the main text (Eq. (2) of the main text), without normalizing it with the biomass before the perturbation.
If we define the local not-normalized resistance as the inverse of the biomass change at a node caused by a perturbation, we can simply define the global not-normalized resistance of the network as the not-normalized resistance of the total biomass which can be rewritten as That is, the global estimate of the not-normalized resistance is equal to the harmonic mean of the local estimates, divided by the number of nodes of the network. Thus, networks would have fundamentally smaller values of this estimate than the nodes forming these networks, so it is not a scale-free property. Figure S1. Spatial random networks employed for the computation of regional estimates of the different stability estimates. (A) Random spatial networks. The nodes and connections between the nodes are constructed as Erdős-Rényi random graphs ( , ), with the number of edges and the probability of presence of each possible edge between the nodes (taken as = 1 10 for all the plots depicted through the text). These graphs were generated using the Networkx package (Hagberg et al., 2008) in Python 3.7 (Python Core Team, 2019). The positions of the nodes (both x and y axes) are assigned from a uniform distribution (0, ), to ensure that bigger random spatial networks cover larger regions. In each of the patches, local dynamics is governed by a Lotka-Volterra Model. Moreover, diffusion is implemented at each time step between each pair of connected nodes. (B) Random spatial dendritic networks, generated in R (R Core Team, 2020) with the OCNet package . These dendritic networks resemble those of riverine structures, and present a higher spatial order than the previous random networks. As before, local dynamics is assumed to be Lotka-Volterra, and diffusion is possible just between connected nodes of the network. Figure S2: Analogous to Fig. 2 of the main text, but for random communities of 10 competitors in 10 node random dendritic networks (Fig. S1B). 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 (panel A) and community growth rate (panel B). Black dots are the growth rates of individual localities (A) or species (B) on which the means AM and AM w were computed. AM w are found to be better estimators of the growth rate, as expected from the results shown in the text Figure S3: Analogous to Fig. 4 of the main text, but for random communities of 10 competitors in 10 node random dendritic networks (Fig. S1B). 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 (panel A) and community resistance (panel B). Black dots are the resistances of the individual localities (A) or species (B) on which the means were computed. HM w are found to be better estimators for the resistance, as expected from the results shown in the text. Figure S4: Analogous to Fig. 5 of the main text, but for random communities of 10 competitors in 10 node random dendritic networks (Fig. S1B). 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). AM and AM w values closer to the identity line (black, dashed) estimate more precisely the regional (panel A) and community initial resilience (panel B). Black dots are the initial resiliences of the individual localities (A) or species (B) on which the means were computed. AM w are found to be better estimators of the initial resilience, as expected from the results shown in the text. Figure S5: Analogous to Fig. 6 of the main text, but for random spatial dendritic networks (Fig.  S1B). Local (A) and regional (B) invariability estimates in random dendritic 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 dendritic networks, for different number of species forming the communities. In panels A and B, we have considered three different scenarios: asynchronous local dynamics ( ̅ = 0, yellow box plots), partially-synchronous local dynamics ( ̅ = 0.44, blue), and perfectly-synchronous local dynamics ( ̅ = 1, green). In panels C and D, we have considered the cases of asynchronous (̃= 0), partiallysynchronous (̃= 0.17) and perfectly-synchronous (̃= 1) species dynamics. Solid lines depict the invariabilities predicted by analytical expressions (Eqs. (20)-(22), and analogous expressions for community invariability) Figure S6: Dependence of the sample size required to control the relative standard error of the network growth rate with the correlation of the nodes' biomasses and growth rates, , , for different coefficients of variation of growth rate ( / ) and biomass ( / ). Confirming results of Fig. 3 of the main text, the required sample size increases mainly with the coefficient of variation of the node growth rates, and to a lesser extent with the coefficient of variation of the node biomasses except for the case of high anticorrelation. Moreover, the required sample size decreases with the correlation of biomasses and growth rates. Figure S7: The analytical approximation for the standard error of regional growth rate (Eq. (11) of the main text) and resistance (Eq. (A5) of Appendix A) compares well to the numerically simulated standard errors obtained with bootstrapping techniques . These results have been obtained for the population dynamics of a unique species in random spatial networks of 20 nodes (Fig. S1A).