Here we develop an approach for predicting temporal stability and resilience from their resistance and recovery components. We began by considering a common measure of temporal stability, the inverse of the coefficient of variation, and defined similar measures of resistance, recovery and resilience that are also the inverse of a deviation. Given that the predicted relationships between integrated measures of stability (that is, temporal stability and resilience) and their component measures of stability (that is, resistance and recovery) will depend on the measures used, we also defined a second set of stability measures as the complement, rather than the inverse, of a deviation. We then showed how temporal stability and resilience can be predicted from their resistance and recovery components, and how resistance can be forecasted from monitoring temporal stability. After generating these new theoretical predictions, we empirically tested them with data from the world’s longest-running biodiversity experiment, allowing us to consider stability at both the ecosystem and the species levels.
Defining two sets of stability measures
We defined temporal stability (also known as invariability) with a common measure, the inverse of the coefficient of variation (CV; D1 in Table 1):
$$I_1\equiv \frac\mu \sigma =\frac1\mathrmCV,$$
and we similarly defined resistance (D2 in Table 1):
$$\Omega _1\equiv \frac\barY_n,$$
and recovery (D3 in Table 1):
$$\Delta _1\equiv |\frac\barY_n-Y_e\barY_n-Y_e+1|,$$
with measures from our previous related work3. We similarly defined resilience (sensu proximity to unperturbed levels after a perturbation5,6,7,8,9,10) as (D4 in Table 1):
$$\Phi _1,x\equiv \frac\barY_n.$$
Here \(\mu \) and \(\sigma \) are the mean and standard deviation, respectively, of the time series and \(\barY_n\), Ye, Ye + 1 and Ye + x are respectively values (for example, ecosystem productivity or other measures of a system level) during normal times (mean across all non-perturbation times), during a perturbation, and during one (or x) time unit (or units) after a perturbation. We refer to these four (D1 to D4) as inverse measures of stability, given that they are all the inverse of a deviation. We used a subscript of one to denote inverse measures of stability.
Given that the predicted relationships between an integrated measure of stability and its components will probably depend on the measures used, we additionally considered a second set of stability measures that define temporal stability (D5 in Table 1):
$$I_2\equiv 1-\frac\sigma \mu =1-\mathrmCV,$$
resistance (D6 in Table 1):
$$\Omega _2\equiv 1-\frac\barY_n,$$
recovery (D7 in Table 1):
$$\Delta _2\equiv 1-|\frac\barY_n-Y_e+1\barY_n-Y_e|,$$
and resilience (D8 in Table 1):
$$\Phi _2,x\equiv 1-\frac\barY_n.$$
We refer to these as complement measures of stability, given that they are the complement (rather than the inverse) of a deviation. We used a subscript of two to denote complement measures of stability.
All of these stability measures (definitions D1 to D8) are dimensionless, and thus can be directly compared among studies and communities with different normal levels, such as different levels of normal productivity, or among data with different units from different dynamical systems. They are also symmetric, and thus can be directly compared between positive and negative perturbations (Fig. 1a), such as wet and dry climate events. Furthermore, these stability measures can be applied to dynamical systems that exhibit either monotonic recovery or damped oscillations after a perturbation.
For both sets of measures, resistance indicates a lack of displacement during a perturbation, quantified relative to the normal level. Inverse resistance (\(\Omega _1\)) is quantified as the inverse of a proportional displacement from normal levels during a perturbation. For example, if the system is reduced by one-quarter of its normal level during a perturbation, then \(\Omega _1\) = 4 (that is, the inverse of one-quarter; Fig. 1a, high-resistance case). Complement resistance (\(\Omega _2\)) is quantified as the complement of a proportional displacement from normal levels during a perturbation. For example, if the system is reduced by one-quarter of its normal level during a perturbation, then \(\Omega _2\) = 0.75 (that is, the complement of one-quarter; Fig. 1a, high-resistance case).
For both sets of measures, recovery indicates the rate of return to normal levels just after a perturbation, relative to the magnitude of displacement. Inverse recovery (\(\Delta _1\)) is quantified as the inverse of a proportional lack of recovery just after a perturbation. For example, if the system recovers all but one-quarter of the way from perturbed to normal levels during the time step following a perturbation, then \(\Delta _1\) = 4 (that is, the inverse of the one-quarter proportion of displacement not yet recovered; Fig. 1a, high-recovery case). Complement recovery (\(\Delta _2\)) is quantified as the complement of a proportional lack of recovery after a perturbation. For example, if the system recovers all but one-quarter of the way from perturbed to normal levels during the time step following a perturbation, then \(\Delta _2\) = 0.75 (that is, the complement of the one-quarter proportion of displacement not yet recovered; Fig. 1a, high-recovery case). Here we focused on the case in which the system returns to normal levels following a perturbation (\(\Delta _1\) > 1, \(\Delta _2\) > 0); however, future studies could additionally consider cases in which a system fails to recover and instead moves away from a repeller or has multiple attractors (for example, \(\Delta _1\) < 1, \(\Delta _2\) < 0).
For both sets of measures, resilience indicates the proximity to normal levels after a perturbation, relative to the normal level. We recognize that the concept of resilience is multifaceted and has multiple definitions45,46. In fact, we were unable to find any definition of resilience that all co-authors of this paper agreed with. In previous work3, some of us have used a narrower definition of resilience that is mathematically independent of resistance and equivalent to how recovery is defined here (D3). Here we instead used the term resilience consistent with several other previous studies5,6,7,8,9,10 where resilience includes both resistance and recovery components. Inverse resilience (\(\Phi _1,x\)) is quantified as the inverse of a proportional deviation from normal levels at some time, x, after a perturbation. For example, if the system is reduced by one-quarter of its normal level during a perturbation (\(\Omega _1\) = 4, \(\Omega _2\) = 0.75, \(Y_e=75\)) and if the system then recovers all but one-quarter of the way from perturbed to normal levels during the first time step following a perturbation (\(\Delta _1\) = 4, \(\Delta _2\) = 0.75, \(Y_e+1=93.75\)), and then recovers all but one-quarter of the remaining distance to normal levels during the second time step following a perturbation (\(Y_e+2=98.44\)), then \(\Phi _1,2=64\) (Fig. 1a, high-resistance and high-recovery case). Complement resilience (\(\Phi _2,x\)) is quantified as the complement of a proportional deviation from normal levels at some time, x, after a perturbation. For example, in the case of high-resistance and high-recovery shown in Fig. 1a, \(\Phi _2,2=0.984\). Our consideration of resilience is a simple discrete-time extension of previous theoretical work in ecology, which considered continuous-time and small perturbations from an equilibrium value11,12. Although less novel, we included resilience in our framework to relate our novel findings for temporal stability and its resistance and recovery components to this earlier work considering resilience.
Complement measures of stability have two advantages over inverse measures: (1) all four complement measures will often be in the range from zero to one with intuitive interpretations (see below), and (2) three of the four complement measures (that is, \(I_2,\Omega _2,\Phi _2,x\)) avoid running off to infinity by avoiding having a difference that can approach zero in their denominator (Extended Data Fig. 1).
Complement measures of stability have intuitive interpretations. The complement measure of resistance indicates the proportion of the normal level retained during the perturbation. For example, \(\Omega _2\) = 0.75 indicates that the system retained 75% of its normal level during the perturbation. The complement measure of recovery indicates that the proportion of the displacement recovered in the time step after the perturbation. For example, \(\Delta _2\) = 0.75 indicates the system recovered 75% of the way back to normal levels just after the perturbation. The complement measure of resilience indicates that the proportion of the normal level at some specified time after the perturbation. For example, \(\Phi _2,2=0.984\) indicates that the system is at 98.4% of its normal level two time steps after the perturbation.
One potential weakness of complement measures of stability is that they can be negative. For complement measures of resilience and recovery, this simply means that the system is moving farther from, rather than closer to, unperturbed levels, or possibly that observation error is overwhelming any recovery. Complement measures of resistance will be negative when the system is displaced by more than 100%, and temporal stability will be negative when the standard deviation is larger than the mean. In the data considered here, only recovery was observed to be negative in some cases (Fig. 2).
Predicting temporal stability from resistance and recovery
Consider a discrete time series of system levels (for example, ecosystem productivities), periodically perturbed (for example, by climate events), with recovery back towards an attractor or away from a repeller after perturbations. If we assume that, during a perturbation, the system level depends only on resistance, then solving for Ye in the resistance formula gives the system level during a negatively perturbed time (for example, productivity during drought), \(\barY_n\) – \(\barY_n\)/\(\Omega _1\), or a positively perturbed time, \(\barY_n\) + \(\barY_n\)/\(\Omega _1\) (Fig. 1a). If we assume recovery occurs at a constant proportional rate, then solving for Ye + 1 in the recovery formula gives the system level during recovery times t either for monotonic recovery, Yt = \(\barY_n\) − (\(\barY_n\) − Yt − 1)/\(\Delta _1\) (Fig. 1a), or for damped oscillations, Yt = \(\barY_n\) + (\(\barY_n\) − Yt − 1)/\(\Delta _1\). Without loss of generality, consider the case in which there is a negative perturbation during time one and recovery via damped oscillations during subsequent times. The first four values of the time series would be given by:
$$Y_1=\barY_n-\frac\barY_n\Omega _1,$$
(M1)
$$Y_2=\barY_n+\frac\barY_n-Y_1\Delta _1=\barY_n+\frac\barY_n-(\barY_n-\barY_n/\Omega _1)\Delta _1=\barY_n+\frac\barY_n\Omega _1\Delta _1,$$
(M2)
$$Y_3=\barY_n+\frac\barY_n-Y_2\Delta _1=\barY_n+\frac\barY_n-\left(\barY_n+\frac\barY_n\Omega _1\Delta _1\right)\Delta _1=\barY_n-\frac\barY_n\Omega _1\triangle _1^2,$$
(M3)
$$Y_4=\barY_n+\frac\barY_n-Y_3\Delta _1=\barY_n+\frac\barY_n-\left(\barY_n-\frac\barY_n\Omega _1\triangle _1^2\right)\Delta _1=\barY_n+\frac\barY_n\Omega _1\triangle _1^3,$$
(M4)
and therefore, by induction, values at time t in this time series can be given by:
$$Y_t=\barY_n+s_t\frac\barY_n\Omega _1\triangle _1^t-1,$$
(M5)
where \(s_t\) is a sign factor indicating whether the system is below or above \(\barY_n\) at time t:
$$s_t=\left\{\beginarraycc-1, & Y_t < \barY_n\\ +1, & Y_t > \barY_n\endarray\right..$$
For monotonic recovery, we take \(s_t=s\) for all t, with \(s\in \-1,+1\\) fixed for a given trajectory (\(s=-1\) for a negative perturbation and \(s=+1\) for a positive perturbation). For damped oscillations, we instead set \(s_t=s(-1)^t-1\) so that the deviation from \(\barY_n\) alternates in sign over time. For the case considered above, which has a negative perturbation followed by recovery via damped oscillations, this gives \(s_t=-(-1)^t-1=(-1)^t\). M5 also applies to positive (rather than negative) perturbations and to monotonic (rather than oscillatory) recovery.
Next, if we further assumed that positive and negative perturbations are symmetric such that the time series mean equals the unperturbed mean, μ = \(\barY_n\), where μ is the mean value across all times, and that perturbations are equally spaced over time with T/P time units between perturbations, where T is the length of the time series and P is the number of perturbations, then the variance of the time series can be given as a function of resistance and recovery:
$$\sigma ^2=\frac1T-1\mathop\sum \limits_t=1^T(Y_t-\barY_n)^2=\fracPT-1\mathop\sum \limits_t=1^T/P(Y_t-\barY_n)^2=\fracPT-1\mathop\sum \limits_t=1^T/P\left(\frac\barY_n\Omega _1\triangle _1^t-1\right)^2,$$
(M6)
where the second form groups time steps by perturbation cycles, and temporal stability can then be specified as a function of its resistance and recovery components:
$$I_1=\frac\mu \sigma =\frac\barY_n{\sqrt{\fracPT-1\sum _t=1^T/P\left(\frac\barY_n\Omega _1\triangle _1^t-1\right)^2}}.$$
(M7)
Examples of dynamics that meet all these assumptions for contrasting levels of resistance and recovery are shown in Fig. 1a. For cases where \(\Omega _1\) > 0, \(\Delta _1\) > 0 and \(\barY_n\) > 0, and given that for any geometric series \(\sum _k=1^n\rmc^k=c\fracc^n-1c-1\), M7 can be further simplified to (E1 in Table 1):
$$I_1=\frac\Omega _1\Delta _1\sqrt\fracP(1-\triangle _1^-2T/P)(T-1)(\triangle _1^2-1).$$
E1 in Table 1 specifies how increasing either resistance or recovery could increase temporal stability, as shown in Extended Data Fig. 1b.
E1 in Table 1 also shows that temporal stability may often depend much more on resistance than on recovery. First, note that recovery contributes very little to temporal stability as the frequency of perturbations increases. That is, as P approaches T (with 0 ≤ P ≤ T), E1 in Table 1 simplifies to the following approximation (A1 in Table 1):
$$I_1\approx \frac\Omega _1\sqrt\fracP(T-1),$$
and temporal stability depends only on resistance. To see this simplification, recall that \(1-\Delta _1^-2\) can be rewritten as \((\Delta _1^2-1)/\Delta _1^2\) and thus \(\Delta _1\sqrt((\Delta _1^2-1)/\Delta _1^2)/(\Delta _1^2-1)=\Delta _1\sqrt1/\Delta _1^2=1\). Furthermore, for frequent perturbations and sufficiently large T (that is, in the limit as \(P\to T\) and \(T\to \infty \)), A1 becomes \(I_1\approx \Omega _1\), and thus temporal stability becomes nearly equivalent in magnitude to resistance. When perturbations are frequent, recovery contributes little to temporal stability because the period of recovery between perturbations is brief. Second, note that recovery also contributes very little to temporal stability when recovery is relatively rapid, even if perturbations are infrequent. When recovery is rapid or perturbations are infrequent, \(1-\triangle _1^-2T/P\) approaches one; and when recovery is rapid, \(\triangle _1\sqrt1/(\triangle _1^2-1)\) approaches one. In these cases, E1 in Table 1 again simplifies to approximation A1. When recovery is rapid, even if perturbations are infrequent, recovery contributes little to temporal stability because the system is essentially fully recovered during much of the period between perturbations. Therefore, recovery will only substantially contribute to temporal stability when the system spends much of its time recovering, which happens when two conditions are met: (1) perturbations are infrequent, and (2) recovery is slow.
Consequently, temporal stability can often be approximated by the resistance of the system to perturbations and frequency of perturbations, irrespective of its rate of recovery following perturbations. For relatively long time series, \(T-1\,\approx T\), and approximation A1 can be further simplified to:
$$I_1\approx \Omega _1/\sqrtp,$$
where p is the proportion of times in which there is a perturbation. This simple approximation may be useful for many dynamical systems.
Similar predictions can be made for complement measures of stability. Using the relationship between the two sets of measures, \(I_2=1-\frac1\rmI_1\), M7 can be rewritten as:
$$I_2=1-\frac1I_1=1-\frac{\sqrt{\fracPT-1\sum _t=1^T/P\left(\frac\barY_n\Omega _1\triangle _1^t-1\right)^2}}\barY_n.$$
(M8)
Substituting \(\Omega _2=1-\frac1\Omega _1\) and \(\Delta _2=1-\frac1\Delta _1\) gives the complement measure of temporal stability as a function of its resistance and recovery components:
$$I_2=1-\frac{\sqrt\fracPT-1\sum _t=1^T/P(\barY_n(1-\Omega _2)(1-\Delta _2)^t-1)^2}\barY_n=1-(1-\Omega _2)\sqrt\fracPT-1\sum _t=1^T/P(1-\Delta _2)^2t-2.$$
(M9)
For cases where \(\Omega _2\) > 0, \(\Delta _2\) > 0 and \(\barY_n\) > 0, and given that for any geometric series \(\sum _k=1^n\rmc^k=c\fracc^n-1c-1\), M9 can be further simplified to (E2 in Table 1):
$$I_2=1+(\Omega _2-1)\sqrt\fracP[(1-\Delta _2)^2T/P-1](T-1)\Delta _2(\Delta _2-2).$$
E2 in Table 1 also shows that temporal stability may often depend much more on resistance than on recovery. First, note that recovery contributes very little to temporal stability as the frequency of perturbations increases. For example, as P approaches T, E2 in Table 1 simplifies to (A2 in Table 1):
$$I_2\approx 1+(\Omega _2-1)\sqrt\fracPT-1$$
and for sufficiently large T, this becomes \(I_2\approx \Omega _2\). Second, note that recovery also contributes very little to temporal stability when recovery is relatively rapid. That is, as \(\Delta _2\) approaches one, E2 in Table 1 again simplifies to approximation A2.
For relatively long time series, \(T-1\,\approx T\), and approximation A2 can be further simplified to:
$$I_2\approx 1+(\Omega _2-1)\sqrtp,$$
These equations and approximations (E1, E2, A1 and A2 in Table 1) can be applied to both pulse and press perturbations because: (1) the proportion of times in which a perturbation occurs can range from 0 to 1, and (2) the time steps can be defined over very brief (for example, daily measurements of soil moisture) or long (for example, multiyear droughts) periods. Note that there may be a similar level of temporal stability resulting from low resistance to an infrequent event (for example, p = 0.01 for a 1 in 100 year drought) or high resistance to a more frequent event (for example, p = 0.1 for a 1 in 10 year drought).
Also, note that these simple approximations (A1 and A2) can be used to estimate temporal stability over any time interval at least as long as the time interval over which the system recovers relatively quickly. For example, temporal stability is often insensitive to recovery rates when \(\Delta _1 > 2\) (Extended Data Fig. 1b). Thus, if it takes 3 years for the system to recover halfway back to normal levels, then it may be possible to approximate temporal stability across intervals of 3 or more years.
In addition, to avoid using a long time series for estimating the unperturbed system level, \(\barY_n\), one can use other data (for example, unperturbed controls in an experiment) or the value observed one time step before the perturbation event, \(Y_e-1\), assuming the system was unperturbed at that time.
Predicting resilience from resistance and recovery
Next, we considered how resilience depends on its resistance and recovery components. M5 can be rearranged to solve for the deviation between normal levels and the level at a specified time of interest (that is, where \(Y_t=Y_e+x\)):
$$|\barY_n-Y_e+x|=\frac\barY_n\Omega _1\triangle _1^t-1.$$
(M10)
Then, given definition D4 above and that x = t − 1 because the perturbation occurred at time one and x counts the time steps after the perturbation (E3 in Table 1):
$$\Phi _1,x=\Omega _1\Delta _1^x.$$
This shows that resilience depends on its resistance and recovery components. Initially, one time step after a perturbation (x = 1), resilience depends equally on resistance and recovery. Longer after a perturbation, resilience depends increasingly and eventually almost exclusively on recovery.
Similar predictions can be made for complement measures of stability. Given that \(\Phi _2,x=1-\frac1\Phi _1,x\), E3 in Table 1 can be rewritten as:
$$\Phi _2,x=1-\frac1\Omega _1\triangle _1^x.$$
(M11)
Substituting \(\Omega _2=1-\frac1\Omega _1\) and \(\Delta _2=1-\frac1\Delta _1\) gives the complement measure of resilience as a function of its resistance and recovery components (E4 in Table 1):
$$\Phi _2,x=1-(1-\Delta _2)^x+\Omega _2(1-\Delta _2)^x.$$
Initially, just after a perturbation (when x = 1), E4 simplifies to \(\Phi _2,1=\Omega _2+\Delta _2-\Omega _2\Delta _2\). Thus, as shown above for the inverse measures of stability, resilience initially depends equally on resistance and recovery, but increasingly depends on its recovery component after there has been more time for recovery following a perturbation.
Forecasting resistance by monitoring temporal stability
Next, we show an example of how the above equations and approximations can be rearranged to solve for other variables of interest. It may be useful to forecast resistance to a perturbation, before it occurs, by monitoring temporal stability. To do so, one could simply rearrange the above approximations A1 and A2, solving for resistance, to get (A3 in Table 1):
$$\Omega _1\approx I_1\sqrt\fracPT-1,$$
and (A4 in Table 1):
$$\Omega _2\approx 1+\fracI_2-1\sqrt\fracPT-1,$$
and then monitor temporal stability to forecast resistance, as empirically demonstrated below.
Predictions that differ between the two sets of stability measures
In the main text, we have emphasized and tested general predictions that are consistent between the two sets of stability measures. Here we note some differences between inverse and complement measures of stability and encourage future studies to test these as well.
For inverse measures of stability, temporal stability is predicted to depend on recovery more at high than at low resistance (Extended Data Fig. 1b). For complement measures of stability, temporal stability is predicted to depend on recovery more at low than at high resistance (Extended Data Fig. 1d). For inverse measures of stability, initial resilience is predicted to increase especially when both recovery and resistance are increased (concave up relationship in Extended Data Fig. 1c). By contrast, for complement measures of stability, initial resilience is predicted to increase when either recovery or resistance are increased (concave down relationship in Extended Data Fig. 1e).
Considering an alternative definition of resistance
Here and in our previous work3, we defined resistance relative to the unperturbed system level, \(\barY_n\), the mean value observed across the unperturbed times. Some previous studies have instead quantified resistance relative to the pre-perturbed system level, Ye − 1 (for example, replacing \(\barY_n\) with Ye − 1 in D2 or D6). The two references will be similar (\(Y_e-1\approx \barY_n\)) when the system has had sufficient time to nearly recover between perturbations. By contrast, when there are multiple sequential perturbations, the two references make very different assumptions about how the system responds. Our definitions above consider how far the perturbation displaces the system from the unperturbed level, regardless of its pre-perturbed level. Alternatively, perturbations may displace the system some distance from its pre-perturbed level, regardless of its unperturbed level. For example, consider what might occur if a slowly recovering system has no recovery time between two consecutive climate events. If the system starts at its unperturbed level during year 1, then decreases to 50% below its unperturbed level during a negative perturbation in year 2, then a positive perturbation in year 3 might increase the system to:
-
(1)
150% of its unperturbed level, decreasing and increasing by the same percentage (that is, 50%) above and below the unperturbed system level (as assumed in D2 and D6);
-
(2)
100% of its unperturbed level, decreasing and increasing by the same absolute amount (that is, 50% of the unperturbed level) during the dry and wet year; or
-
(3)
75% of its unperturbed level, decreasing and increasing by the same percentage (50%) during the dry and wet year, relative to the pre-perturbed level, which is 100% for the dry event and 50% for the wet event.
Our definitions above (D2 and D6) assume the first of these three possibilities. Below, we empirically tested this assumption and found some support for it (see ‘Testing assumptions’ below). Next, we also considered mathematically the third possibility: the system changes by a constant percentage relative to the pre-perturbed system level.
To consider how our main results might differ had we used a different definition of resistance, we determined how temporal stability depends on resistance and recovery when resistance is defined in reference to the pre-perturbation, rather than the unperturbed, system level. To do so, we replaced \(\barY_n\) with Ye − 1 in D2, providing an alternative definition of resistance as:
$$\hat\Omega _1\equiv \fracY_e-1Y_e-1-Y_e,$$
M1–M5 remain unchanged, simply with Ye − 1 in place of \(\barY_n\), and thus M5 becomes:
$$Y_t=Y_e-1+s_t\fracY_e-1\Omega _1\triangle _1^t-1,$$
(M12)
When resistance is defined with respect to the pre-perturbation system level, rather than the unperturbed system level, M6 becomes:
$$\sigma ^2=\fracPT-1\mathop\sum \limits_t=1^T/P\left(Y_e-1-\barY_n+s_t\fracY_e-1\hat\Omega _1\triangle _1^t-1\right)^2,$$
(M13)
and M7 becomes:
$$\hatI_1=\frac\mu \sigma =\frac\barY_n{\sqrt{\fracPT-1\mathop\sum \limits_t=1^T/P\left(Y_e-1-\barY_n+s_t\fracY_e-1\hat\Omega _1\triangle _1^t-1\right)^2}}.$$
(M14)
When the pre-perturbation system level is near the unperturbed system level, \(Y_e-1-\barY_n\approx 0\), the two definitions for resistance are equivalent (\(\Omega _1=\hat\Omega _1\)) and the predictions and approximations above hold (that is, M14 simplifies to E1 and A1 in Table 1).
However, when the pre-perturbation system level is far from unperturbed levels, as might be the case when there are multiple sequential perturbations, this difference cannot be ignored. Furthermore, in this case, the sign factor does not immediately cancel and thus, for monotonic recovery, E1 in Table 1 becomes:
$$\beginarrayl\hatI_1,m=\\ \frac\barY_n{\sqrt\fracPT-1\left[\fracTP(Y_e-1-\barY_n)^2+\frac2sY_e-1(Y_e-1-\barY_n)\hat\Omega _1\,\frac\triangle _1(1-\triangle _1^-T/P)\triangle _1-1+\fracY_e-1^2\hat\Omega _1^2\,\frac\triangle _1^2(1-\triangle _1^-2T/P)\triangle _1^2-1\right]}.\endarray$$
(M15)
and, for damped oscillations, E1 in Table 1 becomes:
$$\beginarrayl\hatI_1,d=\\ \frac\barY_n{\sqrt\fracPT-1\left[\fracTP(Y_e-1-\barY_n)^2+\frac2sY_e-1(Y_e-1-\barY_n)\hat\Omega _1\,\frac\triangle _1(1-(-\triangle _1^-1)^T/P)\triangle _1+1+\fracY_e-1^2\hat\Omega _1^2\,\frac\triangle _1^2(1-\triangle _1^-2T/P)\triangle _1^2-1\right]}.\endarray$$
(M16)
assuming in both cases positive parameters and \(\triangle _1\ne 1\). To obtain M15 and M16, we expanded the squared term in M14 and then summed over all time steps t between perturbations, grouping terms for clarity. This produces three components: (1) the deviation of the pre-perturbation system level from the unperturbed mean, (2) a cross-term representing the interaction between this deviation and subsequent recovery dynamics, and (3) a term describing the recovery trajectory. That is, summing each component separately gives:
$$\mathop\sum \limits_t=1^T/P\left(Y_e-1-\barY_n+s_t\fracY_e-1\hat\Omega _1\triangle _1^t-1\right)^2=\fracTP(Y_e-1-\barY_n)^2+2sY_e-1(Y_e-1-\barY_n)\frac1\hat\Omega _1\mathop\sum \limits_t=1^T/Ps_t\triangle _1^-(t-1)+\fracY_e-1^2\hat\Omega _1^2\mathop\sum \limits_t=1^T/P\triangle _1^-2(t-1).$$
For monotonic recovery, \(s_t=s\) for all t, giving the geometric sum:
$$\mathop\sum \limits_t=1^T/P\triangle _1^-(t-1)=\triangle _1(1-\triangle _1^-T/P)/(\triangle _1-1).$$
For damped oscillations, \(s_t=s(-1)^t-1\), giving an alternating geometric series:
$$\mathop\sum \limits_t=1^T/P(-1)^t-1\triangle _1^-(t-1)=\triangle _1(1-(-\triangle _1^-1)^T/P)/(\triangle _1+1).$$
Substituting these into M14 gives the two closed forms M15 and M16, respectively.
If perturbations are symmetric, such that positive and negative deviations are equally frequent and the time-series mean equals the unperturbed mean (μ = \(\barY_n\)), the contributions of the sign factor cancel when summed over all perturbations. We note that, over infinite time, \(\mu \ne \barY_n\) with this definition of resistance (\(\hat\Omega _1\)) and these dynamics (M12). However, for finite time, μ can approximately equal \(\barY_n\), especially when the system nearly recovers between perturbations (that is, rapid recovery, infrequent perturbations, or both). With the symmetry assumption, the cross-term vanishes, and both M15 and M16 simplify to the following general expression:
$$\hatI_1=\frac\barY_n{\sqrt\fracTT-1\,(Y_e-1-\barY_n)^2+\fracPT-1\,\fracY_e-1^2\hat\Omega _1^2\,\frac\triangle _1^2(1-\triangle _1^-2T/P)\triangle _1^2-1\,}.$$
(M17)
Thus, M17 describes both monotonic and oscillatory recovery dynamics under symmetric perturbations with resistance defined by the pre-perturbation, rather than the unperturbed, system level.
This general expression can be further simplified as recovery becomes rapid (\(\triangle _1\to \infty \)), perturbations become frequent (\(P\to T\)), or both:
$$\frac\triangle _1^2(1-\triangle _1^-2T/P)\triangle _1^2-1\to 1,$$
the recovery term vanishes, and M17 simplifies to the following approximation (in place of A1 in Table 1):
$$\hatI_1\approx \frac\barY_n{\sqrt\fracTT-1\,(Y_e-1-\barY_n)^2+\fracPT-1\,\fracY_e-1^2\hat\Omega _1^2\,}.$$
Therefore, considering this alternative definition of resistance provides further support that temporal stability can be approximated by resistance alone when recovery is rapid, perturbations are frequent, or both. Only when perturbations are infrequent and recovery is slow does recovery substantially influence temporal stability. These results parallel those derived above for M7, and E1 and A1 in Table 1, confirming that they hold even when resistance is defined relative to the pre-perturbation system level, rather than the unperturbed level.
Study site
Our field experiment is located at Cedar Creek Ecosystem Science Reserve, near East Bethel, MN, USA. The climate of the study site is continental, with mean January temperature of −11 °C, mean July temperature of 22 °C, and mean annual precipitation of 660 mm per year47. The field site is located on a glacial outwash plain and has sandy and nutrient-poor soils. Before European colonization, the experimental site was an oak savanna, with wolves (Canis lupus), bison (Bison bison), elk (Cervus elaphus) and other animals and plants that later became extirpated. Subsequently, the site was used for various row crops for an unknown number of years48,49. Between 1956 and 1968, fields on this part of the Reserve were abandoned from agriculture and became old fields undergoing succession.
BioDIV experiment
We tested these new theoretical predictions with data from the world’s longest-running biodiversity experiment: the BioDIV experiment (E120 of the Cedar Creek Long-Term Ecological Research program)24,26,27,28,29,30,50. In 1993, in a 7-ha area: (1) glyphosate herbicide was applied; (2) the dead vegetation was burned; (3) the top 6–8 cm of soil was removed by a bulldozer to reduce the seed bank, which also reduced soil carbon and nutrients; and (4) the remaining soil was ploughed and harrowed. This site management homogenized the 7-ha area, in preparation for planting the BioDIV experiment. In 1994 and 1995, different numbers and combinations of grassland plant species were randomized to replicate plots sown with 1, 2, 4, 8 or 16 plant species from a pool of 18 species, including 16 herbaceous species and 2 oak (Quercus) species. The main experiment consists of 154 plots, each 9 × 9 m, that have been burned each spring and hand-weeded each summer to remove non-sown species. Additional experimental design details are available online (see E120 on Cedar Creek’s Long-Term Ecological Research website: https://www.cedarcreek.umn.edu).
Since 1996, peak aboveground biomass has been clipped in two strips, each 10 × 600 cm, in each of the 154 plots that have been maintained. Clipped biomass was dried to a constant mass and weighed. In these temperate grasslands, aboveground peak biomass is a proxy for annual aboveground net primary productivity because all aboveground tissues die in the winter season and are removed by annual spring burning. Plant aboveground biomass data and metadata for this experiment are available from the Environmental Data Initiative51.
Identifying wet and dry climate events
We defined all wet and dry events based on the standardized precipitation–evapotranspiration index (SPEI)31. Drought occurs when water availability remains below normal levels over some period of time31. Several drought indices consider water balance, accounting for both water inputs (precipitation) and water losses (potential evapotranspiration (PET)). Warming temperatures can contribute to drought by increasing PET. Here we used SPEI to identify and quantify the magnitude of wet and dry climate events (Extended Data Fig. 2). SPEI is a standard normal variable for water balances aggregated over a given number of months at a particular location. SPEI values are based on month-by-month variations in climate observed since 1901, based on monthly means of measurements made at more than 4,000 weather stations worldwide, and provided on 0.5 × 0.5 degree grids globally. For example, our study site experienced a value of SPEI-04 = 1.9 for August of 2002 (Extended Data Fig. 2), which corresponds to a 1-in-35-year wet event for the 4 months (thus the 04 in SPEI-04) of the growing season before the biomass harvest (that is, May, June, July and August), relative to the water balances observed for those 4 months at our site since 1901.
We obtained SPEI values from SPEIbase (v2.9)52 for our study site. We used SPEI-04 throughout the article because we found stronger relationships with growing season, rather than annual, water balances at our study site, probably partly because there is a short growing season and sandy, well-drained soils, making the system most sensitive to precipitation that occurs during the summer months while the plants are actively growing.
We classified experiment years as extremely dry, moderately dry, normal, moderately wet and extremely wet, as in our previous work3 (Extended Data Fig. 2). Extreme events (extremely dry or extremely wet) were defined as those that historically occurred less frequently than once per decade. Moderate events were defined as those that historically occurred between once in 4 years and once per decade. Normal years were defined as those within the interquartile range of historical water balances that did not follow an extreme event (Extended Data Fig. 2).
Statistical analyses
All analyses were conducted in R (v4.4.3)53 using the tidyverse (v2.0.0)54, DImodels (v1.3.2)38 and ggpubr (v0.6.0)55 packages.
For each BioDIV plot, we log transformed plot-level productivity and then calculated D1 to D8 in Table 1. That is, to help meet the assumption of the theoretical predictions that positive and negative perturbations are symmetric such that μ = \(\barY_n\), we first log transformed productivity. Temporal stability measures were calculated for two different time intervals: 1996–2022 and 1996–2020, with the latter used for forecasting. For all measures, \(\barY_n\) is the mean plot-level productivity, averaged across all normal years that did not follow an extreme event (Extended Data Fig. 2). Resistance measures were quantified for the 2002 wet event and the 2021 drought (that is, Ye is the plot-level value observed in 2002 or 2021), with the latter used for forecasting. Recovery measures were quantified for 2003, the year after the wet event (that is, Ye + 1 is the plot-level value observed in 2003). Resilience measures were quantified for 2004, 2 years after the wet event (that is, Ye + 2 is the plot-level value observed in 2004).
Measures of stability that have a difference in their denominator tend to run off to infinity as this difference approaches zero, creating outliers. This issue arises for all the inverse measures of stability, but only for the complement measure of recovery (Table 1). Thus, here we largely avoid the issues of outliers by considering complement, rather than inverse, measures of stability. Although there are no perfect solutions, the influence of outliers can be assessed by censoring or transforming the data, excluding outliers, or through other regression approaches (for example, generalized linear models, beta regression and generalized additive models). Here we censored the complement measure of recovery to a minimum value of zero, which assumes that individual plots are not moving farther from normal the year after a perturbation. Our main results were qualitatively similar when we instead skipped this censoring step (although, in this case, parameter estimates are influenced by outliers), excluded these negative outliers (although, in this case, sample size is reduced to n = 99 plots), or used minimum–maximum normalization (although, in this case and with other transformations, there is an altered interpretation of parameter estimates). When considering inverse measures of stability, it might be helpful to assume that the system moves at least 1% from normal levels during a perturbation (\(\Omega _1\le 100\)) and that the system recovers no more than 99% of the way back to normal levels just after a perturbation (\(\Delta _1\le 100\)).
Next, to account for observation error when estimating stability components, and to get species-level estimates for all measures, we fit diversity–interactions (DI) models36,37 using least squares estimation via the autoDI function in the DImodels package38. DI models estimate species-specific effects, known as species identity effects in the modelling framework, and species interaction effects.
DI models have flexibility in the specification of species interactions. For example, the autoDI function compares the fits of models that assume: (1) there are no species interactions; (2) all pairs of species interact with the same strength (for example, interspecific interactions differ from intraspecific interactions); (3) each pair of functional groups has a unique pairwise interaction; or (4) each pair of species has a unique pairwise interaction. The second option was the most parsimonious for most of our measures of stability (Extended Data Table 1). For this case, where interspecific interactions differ from intraspecific interactions and are similar between all pairs of species, the DI model is:
$$y=\mathop\sum \limits_i=1^S\beta _ip_i+\delta \sum _1\le i < j\le S(p_ip_j)^\theta +\epsilon $$
where y is the community-level response (for example, biomass for a monoculture or mixture plot), \(S\) is richness (the number of species), \(p_i\) and \(\beta _i\) are the initial proportion and identity effect of species i, respectively, \(\delta \) is the interaction effect for all pairs of species, \(\theta \) determines the shape of the relationship between biodiversity and ecosystem functioning or stability (for example, linear, decelerating and saturating), and \(\epsilon \) is a normally distributed random error term with constant variance.
Parameter estimates from a DI model can be used to predict the response for any combinations and proportions of species. For monocultures, predictions rely only on the species identity effect in the model (for example, the prediction for species i in monoculture is \(\haty=\hat\beta _ip_i=\hat\beta _i\times 1=\hat\beta _i\)). For mixtures, predictions are determined by both the identity effects of the component species and any interaction effects (for example, for a two-species equi-proportional mixture of species 1 and 2: \(\haty=\hat\beta _1\times 0.5+\hat\beta _2\times 0.5+\hat\delta (0.5\times 0.5)^\hat\theta \)). Here we used these models to predict species-level estimates of stability measures, which are the estimated values for monocultures of each of the 16 study species, as well as the estimated stability measures for each of the 154 experimental plots. A summary of the most parsimonious models is shown in Extended Data Table 1. Species-level predictions are the identity effects estimated in these models. This gives, for each species, an estimated value in monoculture, accounting for observation error. Ecosystem-level predictions are the plot-level estimates from these models.
Using predictions from DI models to estimate stability measures, instead of raw data, can help to reduce the influence of outliers that may arise from observation error. This is particularly helpful when using stability measures that may otherwise run off to infinity for resistance or recovery, which then leads to unrealistic predictions for temporal stability and resilience.
Note that a DI model may give a negative prediction for a community, if identity effect and/or interaction effect estimates are negative. For instance, for recovery, we found that the most parsimonious model produced negative estimates for the identity effects of two species (Extended Data Table 1). However, note that these estimates were not significantly different from zero (compare these estimates and their standard errors in Extended Data Table 1). Thus, this simply means that, in the time step after a perturbation, there was no evidence that these species significantly recovered (or significantly moved away from normal levels).
Note that our stability predictions that follow can be applied to many types of experimental or observational data, and to other types of perturbations, and thus they do not require the data from biodiversity experiments, the consideration of climate extremes, or the DI modelling framework described above. Here we tested our theoretical stability predictions by using data from a long-term biodiversity experiment because this allows us to consider stability at two levels of biological organization: species and ecosystem. We considered extreme climate events because, in grasslands, these perturbations contribute significantly to the temporal dynamics of species and ecosystem productivity. We used DI models because they are designed for the analysis of data from biodiversity experiments. We encourage future studies to apply our stability predictions below to other observational or experimental data, to carefully consider which types of perturbations are important in those systems, and to use, as a first step, a statistical modelling framework that is suited for the data and study design.
Using these plot-level and species-level estimates of resistance and recovery, we then predicted temporal stability and resilience using E1–E4 and A1–A2 in Table 1. Using plot-level and species-level estimates of temporal stability from 1996 to 2020, we used approximations A3–A4 to forecast resistance to the 2021 drought. In all cases, P = 16 non-normal years and T = 27 total years.
We chose P based on previous results, which have shown that productivity responds not only to extreme and infrequent climate events but also to more moderate and frequent wet and dry events3. Note that E1 and A1 in Table 1 can account for the fact that systems may be less resistant to infrequent than frequent perturbations. For example, the same estimate of temporal stability could arise from an observation of either relatively low resistance to an extreme 1-in-100-year drought (P/T = p = 0.01) or relatively high resistance to a more moderate 1-in-4-year drought (P/T = p = 0.25). Thus, we could have chosen different p estimates based on the extremity of the wet event in 2002 and the drought in 2021. However, updating p in this way would assume that resistance scales with the inverse of the square root of the frequency of perturbations, which may not always hold. In our case, updating p in this way would lead to extreme overestimates of ecosystem stability, probably because we know from previous results that productivity responds not only to these extreme and infrequent events but also to more moderate and frequent wet and dry events3. In other words, in our system, productivity is frequently perturbed by climate events, rather than infrequently perturbed with substantial periods of recovery between perturbations.
Next, we tested for correlations between each integrated measure of stability and its resistance and recovery components. We fit linear models (black lines and confidence intervals in all figures) using the lm function in the stats package via the geom_smooth function in the ggplot2 package and added P values and R2 values with the stat_cor function in the ggpubr package55.
Finally, we tested for relationships between the observed and predicted or forecasted values. For comparability and consistency across all relationships, here we used sequential (type I) sums of squares for all relationships. For temporal stability, a case may be made for plotting the observed values on the x axis rather than the y axis, given that we expect much greater uncertainty in the predicted values, which are based on single-year estimates, than the observed values, which are based on a quarter century of data. We did not do this here, however, because the same argument would not apply to resilience and we wanted to compare with other predictors, such as species richness. For resilience, future studies may consider using major axis regression (type II sums of squares) to account for substantial uncertainty in both the predicted and the observed variables, which are based on only 2 or 1 year of data, respectively. We did not do this here, however, because this is only recommended for sample sizes much larger than those that we were able to consider for species-level relationships. Although there is no perfect solution for all the relationships that we considered, for comparability and consistency, here we used sequential (type I) sums of squares and observed versus predicted values for all relationships. We did not focus on the estimated slopes and intercepts for these relationships because they depend on the choices described above (that is, type I versus type II sums of squares and observed versus predicted, or vice versa). Instead, we quantified error in the predictions and forecasts by quantifying the mean absolute percentage error and the R2 values for each relationship.
Testing assumptions
We tested whether our data met several assumptions of our model. First, consistent with our assumption that there are symmetrical responses to perturbation, we found no significant difference in resistance between positive (wet) and negative (dry) perturbations (Extended Data Fig. 5a). Using the measure of resistance defined here and SPEI, we previously found symmetric responses of grassland plant productivity to observed wet and dry events across 46 experiments3. By contrast, a separate meta-analysis of 83 studies found that productivity often responded more to experimental wet than to experimental drought treatments, with the magnitude of wet or dry events defined by the percentage change in precipitation56. These previous results are not necessarily inconsistent, given the many differences between how wet and dry events were defined and whether climate events were observed or experimentally created. Instead, these contrasting results indicate that the symmetry or asymmetry of responses to wet and dry events probably depends on how these events are defined and imposed. We suspect the symmetry assumption of our analysis can often be met by log-transforming productivity and defining wet and dry events based on SPEI, which accounts for water balance, including both water inputs and water losses due to PET3,56.
Second, consistent with our assumption that resistance can be defined with respect to the average unperturbed system level (\(\Omega _1\)), rather than the level at the pre-perturbation time step (\(\hat\Omega _1\)), we found that resistance — quantified with D2 — does not depend on the level of biomass observed in the year before the perturbation (Extended Data Fig. 5b). We were only able to test this for wet events, given the small number of dry events observed during the study (Extended Data Fig. 2a). Furthermore, see above (‘Considering an alternative definition of resistance’) for further evidence that the deviation of the state variable from the mean after a perturbation does not depend on the value of the state variable before the perturbation.
Third, consistent with our assumption of a constant recovery rate, we found that resilience two time steps after the perturbation can be predicted from resistance during the perturbation and recovery the year after it at the ecosystem level (Fig. 3a). This prediction assumes that the recovery rate during the second year after the perturbation is the same as was observed during the first year after the perturbation. We are unable to have a longer-term test of the assumption of a constant recovery rate because there are no extended periods of recovery with many years between climate events (Extended Data Fig. 2). For other empirical dynamical systems with extended periods of recovery, we encourage additional tests of the assumption of a constant recovery rate.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
