An invariability-area relationship sheds new light on the spatial scaling of ecological stability

The spatial scaling of stability is key to understanding ecological sustainability across scales and the sensitivity of ecosystems to habitat destruction. Here we propose the invariability–area relationship (IAR) as a novel approach to investigate the spatial scaling of stability. The shape and slope of IAR are largely determined by patterns of spatial synchrony across scales. When synchrony decays exponentially with distance, IARs exhibit three phases, characterized by steeper increases in invariability at both small and large scales. Such triphasic IARs are observed for primary productivity from plot to continental scales. When synchrony decays as a power law with distance, IARs are quasilinear on a log–log scale. Such quasilinear IARs are observed for North American bird biomass at both species and community levels. The IAR provides a quantitative tool to predict the effects of habitat loss on population and ecosystem stability and to detect regime shifts in spatial ecological systems, which are goals of relevance to conservation and policy.

The five inserted plots shows the resulting 500 IARs for respective continents, and the colored lines and shade represent the median and the 25% and 75% quantiles, respectively.

Deriving the asymptotic slope of IAR
Here we derive the asymptotic slope (zasym) of IAR by calculating the log-log slope of IAR as the area A goes to infinity: Eq. (S1) shows that, in order to obtain zasym, we need to derive ̅ , i.e. the average pairwise correlation. Consider a square area A = N 2 (i.e. a group of grids on {1, 2, …, N} × {1, 2, …, N}), we have: where ρ(.) represents the correlation-distance relationship. Below we derive ̅ and zasym under two correlation functions.
(ii) Power function: Under the power-decay function, we have: From the above equation, we have: we have: c 2 A −α/2 < ̅ < c 1 −α/2 . Substitute into Eq. S1, we have: As A goes to infinity, both sides of Eq. S7 converge to the below limit: which is independent of c. Therefore, we have:

The triphasic curve of IAR
In the main text, we show that under the exponential correlation-distance function, IARs exhibit triphasic curves. As we have demonstrated above, IARs have asymptotic slopes of 1 under this scenario. This explains the third-stage steep increase of stability with area. Below we make some explanations on the first two stages of such triphasic curves.
We first consider a special scenario, in which between-patch correlation always equal ρ1 regardless of distance (where 0 < ρ1 < 1). Following Eq. 1 in the main text, IAR can be expressed as: As shown in the Methods, the initial slope is given by: zini = log2(2/(1+ ρ1)), which lies between 0 and 1. As the area A goes to infinity, we have: I(∞) = I1/ρ1. In other words, invariability increases with area at the beginning and converges to a constant when A is very large, i.e. the asymptotic slope of IAR is 0 (Supplementary Note Fig. 1a). This result can be understood more intuitively in an ecological context. Consider a landscape with regularly distributed square patches. The biomass dynamics of each local patch i (Xi) is governed by both environmental and demographic stochasticity: Here, n0 is the temporal mean biomass; E and Di are random variables (mean: 0; variance: ve and vd), which represent, respectively, landscape-level environmental stochasticity and patch-level demographic stochasticity. We assume E and Di are independent of each other, and Di and Dj are also independent. Intuitively, we can expect that the effect of demographic stochasticity decreases with the area considered. So, as the area becomes large, the invariability of biomass should converge at some point that simply reflects the effect of environmental stochasticity. In mathematical language, this writes: the invariability of one patch is: S1 = n0 2 /(ve + vd); invariability of two patches is: S1 = n0 2 /(ve + vd/2); … invariability of A patches is: n0 2 /(ve + vd/A). Thus, IARs exhibit a faster increase with area at the beginning; the increase slows down, and the invariability converges to n0 2 /ve when A is large.
Now we come back to the scenario of exponential correlation-distance functions. Within the area L 2 , the between-patch correlation declines with distance very slowly (e.g. at a distance L, the between-patch correlation is 0.37×ρ1, still in the same magnitude of ρ1). Thus, if L is large, we can expect that IARs exhibit a fast increase at the beginning and followed by a "flat phase" near L 2 . Note that because the between-patch correlation is often lower than ρ1 (though always at the same magnitude with ρ1), invariability at this "flat phase" is larger than 1/ρ1 (Supplementary Note Fig. 1b).

Supplementary Note 2: Sampling issues in IAR
In empirical analyses, the sampling scheme related to both spatial and temporal scales may affect the calculation of invariability and hence IAR. The main focus of this study is on spatial scale.
Spatial scale can be represented by three aspects: extent, resolution, and sampling intensity.
Since our study explores space continuously, the extent and resolution are represented by the largest and smallest area, respectively, on the x-axis of IAR. For the third aspect, i.e. sampling intensity, our paper has examined IARs on both continuous landscapes (i.e. flush grids), in our model and in the primary productivity data, and non-continuous ones (i.e. spatially separated grids due to incomplete sampling), in the bird data. In this Supplementary Note, we will explore the influence of sampling intensity on IAR with our theoretical model. As we will show, incomplete spatial sampling could potentially increase the slope of IAR.
Temporal sampling scheme might also influence the calculation of invariability and its scaling patterns, as previous studies suggested 1-3 . Temporal scale can similarly be represented by three aspects: observation length, sampling resolution, and sampling intensity. In this Supplementary Note, we will investigate the effect of observation length on IARs of bird biomass. As we will show, observation length can slightly alter the intercept and slope of IAR.
However, we have not investigated the influences of sampling resolution and intensity in our study. As for the sampling resolution, our study has fixed it to one year due to both data limitation (e.g. bird biomass data is collected once per year) and research interest (e.g. we are interested in the interannual dynamics of NPP, not seasonal oscillations). As for the sampling intensity, we have fixed it to be annually. This is simply because we want to use most of available information in the data, given the relatively short time series (i.e. NPP data: 15 year; bird data: 21 years). This said, investigations into these aspects of temporal scales may be useful for future studies on IAR.

The effect of incomplete spatial sampling on IAR
In our models, we construct IAR based on a "full observation" of the landscape. However, in reality, field investigations may cover only a small proportion of the full landscape, for instance the North American Breeding Bird Survey in our empirical analyses. Here we explore the impact of incomplete sampling on IAR by constructing IAR based on a randomly sampled proportion of the landscape.
As in the main text, we consider a two-dimensional landscape (e.g. 128×128 grids), in which local ecosystem dynamics have identical temporal mean and variability and the between-patch correlation decays with distance following a power function: ρ(d) = ρ1×d -α . We randomly sample a proportion of grids (e.g. 1/4, 1/16, 1/64, etc.) from the whole landscape, which can be regarded as the sampled area in empirical studies. We construct IARs based on these sampled grids, following similar procedures as in our empirical analyses of the bird data. Specifically, starting from one (sampled) grid, we increase the number of grids by including the closest neighbor (sampled) grids. We calculated the temporal invariability for each respective "area", thus generating an Invariability-Area Relationship. Note that just as for our analysis of BBS data, the "area" here is the sampled area, or the number of grids.
We found that incomplete sampling tended to increase the slopes of IARs ( Supplementary   Fig. 2). We rescale the sampled area (or number of grids) by multiplying it with the average "represented area" of each sample (i.e. 4, 16, 64, respectively), which represents the spatial extent of the sampling efforts ( Supplementary Fig. 2). Note, however, this rescaling does not alter the slope between invariability and area on a log-log scale (but it does alter the intercept).

The effect of observation length on IAR
The observation length (i.e. the length of the census) may also affect temporal invariability and  Fig. 3a). Also, under the shorter observation period, correlation decays with distance more slowly (-0.40 vs. -0.55) (Supplementary Fig. 3b). Specieslevel IARs generally exhibit higher intercepts under the shorter observation period, but no systematic increase or decrease for the slope (Supplementary Fig. 3c,d).
Therefore, a longer observation period can decrease the intercept of IAR, but its effect on the slope of IAR is context dependent. To better understand these patterns, future models should incorporate autocorrelation in environmental fluctuations (e.g. red noise) and investigate how it interacts with local population dynamics and dispersal to regulate IARs.

The effect of grain size on correlation-distance relationship
In our model, we have assumed two correlation-distance functions to calculate IARs on a twodimensional landscape. These correlation-distance functions describe how the correlation between two unit-size patches changes with distance. Here we show that the between-patch correlation also depends on the grain size.
To illustrate this, we consider four patches. Two of them (X1 and X2) are neighbors and located around point X (biomass are denoted as Nx1 and Nx2, respectively), and the other two (Y1 and Y2) are neighbors and located around point Y (biomass: Ny1, Ny2) (see Supplementary Note Figure 2). The distance between X and Y (dxy) are much larger than unit length (which is the distance between X1 and X2 or Y1 and Y2). The correlation between patches X1 and X2 or between patches Y1 and Y2 is ρ1, and the correlation between patches at location X and Y (e.g. X1 and Y1) is ρdxy. In other words, at a grain size of 1, the between-patch correlation at a distance d is ρdxy.
Now we calculate the between-patch correlation at a distance d when the grain size is 2.
Supplementary Note Figure 2: Four patches on a two-dimensional landscape. The patches X1 and X2 are neighbors and located around point X, and the patches Y1 and Y2 are neighbors and located around point Y.

Supplementary Note 3: Ecological determinants of IAR
Our theoretical and empirical analyses show that the shape and slope of IAR are mainly determined by patterns of spatial synchrony of ecosystem dynamics. As spatial synchrony increases (e.g. a higher local correlation or a slower correlation decay with distance), invariability increases with area more slowly (Supplementary Fig. 1). In our model, two correlation-distance functions ("short-tail" and "long-tail") are investigated. Although very simple, they seem to capture the spatial synchrony patterns of plant and bird communities at continental scales, which in turn explain the patterns of IARs in these two systems (Figs. 2 and   3).
Previous studies have highlighted several factors that regulated spatial synchrony. In general, spatial synchrony in populations or ecosystem dynamics increases with environment correlation, dispersal, and community similarity [4][5][6][7] . Trophic interactions may also increase spatial synchrony through phase locking 8 . In our model, these factors are incorporated implicitly into scenarios of spatial synchrony (e.g. correlation-distance functions). However, a further understanding of IAR will require mechanistic models that incorporate lower-level abiotic and biotic processes (e.g. environment, dispersal, etc.) in a dynamical landscape, and clarify how they interact and regulate spatial synchrony and consequently IAR. In particular, the relative importance of these factors is likely scale dependent, and one important but challenging task is to tease apart their effects across scales. Such theoretical insights will advance our understanding on empirical patterns of IAR, e.g. the triphasic curve in Fig. 2 see main text. However, despite ignoring spatial heterogeneity, our model predictions were well supported by the data (e.g. exponential correlation-distance relationships generate triphasic IAR (comparing Fig. 1 vs. Fig. 2), and IAR slope decreases with local correlation but increase with the exponent of power-law decay (comparing Supplementary Fig. 1 vs. Fig. 3)). This said, it would be useful for future theoretical work to clarify how spatial heterogeneity may affect the scaling of invariability.
The overharvesting model is a minimal model that can have alternative stable states 10 , and its spatial extension (i.e. equation S11) also exhibits alternative stable states 9 . In the absence of dispersal and environmental perturbations, a local population shifts from an underexploited state to an overexploited one as the harvesting rate (c) increases (Supplementary Figure 4a). Due to dispersal between patches and despite the fact that not all patches reach their own tipping point at the same time, a global tipping point exists at which all populations collapse onto the overexploited state. Prior to the shift, a slower speed of recovery following perturbations together with a spatial coupling allows perturbation to transmit to distant patches, thus increasing the system's spatial synchrony and altering the shape of IAR.
With this model, we explore how IAR changes as the harvesting rate (c) increases and approaches its critical value at the global tipping point (around c* = 2.0615 for the set of parameters used in our simulations). Away from the regime shift, IAR exhibits a triphasic curve, with steeper increases in invariability at both small and large scales ( Supplementary Fig. 4b).
Approaching the global tipping point, local population dynamics exhibit larger variability, as previous studies showed 11 . This causes the intercept of IAR (i.e. local invariability) to decrease prior to the regime shift ( Supplementary Fig. 4b, d). Moreover, a slowing down of the dynamics in each patch allows perturbations to propagate across space, increasing the correlation between neighboring and distant patches 12 . As a consequence, the initial and final slopes of IAR both decrease ( Supplementary Fig. 4c,e,f) so that the triphasic shape gradually diminishes ( Supplementary Fig. 4c). Remarkably, because spatial correlations do not increase in the same way at all scales but gradually propagate through space, the decrease in the initial slope of IAR occurs first, followed by, close to the shift, a decrease in the final slope ( Supplementary Fig.   4e,f). Thus, monitoring IAR could add to the arsenal of early warning indicators of regime shifts.