Ecological resilience: what to measure and how

Intuitively, resilience is the ability of a system to cope with disturbances, bounce back, and maintain its state and functionality. In the ongoing context of global change, understanding resilience is of utmost importance to achieve sustainable interactions between humans and ecosystems (Cañizares et al 2021). However, moving from intuition to practically measuring resilience has been a real challenge in ecology (Carpenter et al 2001, Kéfi et al 2019, Pimm et al 2019, Capdevila et al 2021).

Measuring resilience in practice has been challenged by the fact that the definition of resilience has lost clarity through time. Resilience has been used across multiple scientific disciplines (Baggio et al 2015), each with a different understanding of what resilience means (Angeler and Allen 2016, Walker 2020). Within the ecological literature, resilience has been defined in multiple ways (Grimm et al 1997), and sometimes the same definition has even been applied in different ways across different ecosystems (Yi and Jackson 2021). Now—at least in ecology—two definitions of resilience dominate in the literature: engineering resilience, that is the rate with which a system returns to a reference state after a disturbance (Pimm 1984), and ecological resilience, that is the magnitude of disturbance that can be absorbed before a system flips to another state (Holling 1973).

There have been numerous efforts aiming at bridging the gap between the clear intuitive concept of resilience and an operational and measurable quantity. These efforts can be summarized in identifying properties that guarantee resilience (e.g. Folke et al 2004, Thrush et al 2009, Biggs et al 2012), suggesting surrogates that could indirectly reflect resilience (Carpenter et al 2005, Cumming et al 2005), or developing qualitative and quantitative approaches for evaluating specific aspects of resilience (Sundstrom et al 2014, Quinlan et al 2016) in ecological and socio-ecological systems. For instance, loss of redundancy and response diversity (Folke et al 2004) or functionality of keystone species (Thrush et al 2009) in an ecosystem are seen as properties that could jeopardise resilience. Identifying and assessing the potential of a socio-ecological system to change under the impact of specific drivers and perturbations can be used to produce surrogates of the system's resilience (Cumming et al 2005). Discontinuities in the distribution of ecosystem functions across spatial scales (Sundstrom et al 2014), or mapping system boundaries, dynamics, cross-scale interactions to other systems, and the system's adaptive capacity are used for assessing the resilience of a system (Quinlan et al 2016). These efforts highlight that there are different approaches in assessing resilience that depend on the system in question and the type of stress or disturbance a system experiences (Carpenter et al 2001). Yet, it is one thing to assess resilience based on system properties that promote resilience and another to measure actual resilience based on system responses. As a result, quantifying resilience in practice remains for some a struggle and for others an abstract metaphor.

Here, we revisit the concept of resilience as defined by Holling in 1973, we list metrics used to quantify it, and explore how these metrics relate to each other in an attempt to address the question of how resilience can be estimated in real systems. We are not disillusioned that there is a definite answer to this challenge nor that we are the first to do this. Mathematical descriptions of Holling's resilience have appeared almost as early as the concept was suggested (Grumm 1976). Also, we do not claim that there is a single notion of resilience. In the same way as its parent term 'stability', resilience is a multifaceted concept that can be quantified in different ways. In particular, we do not consider here all possible ways for assessing resilience as they have been developed across ecological and social sciences (Gunderson 2000, Carpenter et al 2001, Cumming et al 2005, Thrush et al 2009, Folke 2016). Instead, we focus on system responses—not properties—that we can use to quantify resilience following Holling's definition based on the concept of the stability landscape and the geometrical properties of its basins of attraction. Our aim is to synthesize what is known about resilience metrics that are related to the existence of a system's stability landscape by clarifying the connections between ecological and engineering resilience metrics and by reviewing current approaches for measuring resilience in models and data.

The paper is structured as follows. We first introduce the two basic definitions of ecological and engineering resilience. We then present, classify, and identify the relationships between ecological and engineering resilience metrics in simple, widely used ecological models. We list tools and approaches for measuring resilience in models and data and review their use in the ecological literature. We conclude by outlining challenges that, if addressed, could help make progress toward more reliably quantifying resilience in practice.

Two main definitions of resilience pervade the ecological literature. The first definition describes the stability of a system near an equilibrium state. It is typically measured as the ability of an ecosystem to recover to its original state after a small perturbation and the speed at which it does so (Pimm 1984). It is a classical measure of local stability. This has also been referred to as 'engineering resilience' in the ecological literature (Holling 1996). The second definition focuses on cases where disturbances are not small and thus the system may not return to its current equilibrium but it may flip to another state. Resilience can then be evaluated as the magnitude of disturbance that can be absorbed before the system changes state (in terms of structure and/or functioning). This has been termed 'ecological resilience' (Holling 1973) and it is typically referred to as 'Holling's resilience'.

The major difference between these two definitions of resilience is the underlying assumption that systems can either have a single (globally) stable state or multiple (locally) stable states (figure 1). When proposed by Holling in 1973, the concept of ecological resilience was introducing a new dimension of ecological stability contrasting with the more traditional concept of local stability, which assumes that systems have a single attractor to which they return to following any type of perturbation. Holling's concept reflected an extension of Odum's paradigm of a homeostatic mechanism that acts as an ecosystem regulator at equilibrium (Odum 1969), while at the same time formalising earlier concepts of threshold changes and buffers of population dynamics (Errington 1945). Over time, the two definitions of resilience became almost two different 'world-views' of ecological stability (Van Meerbeek et al 2021), occasionally leading to different, sometimes confusing (e.g. Hodgson et al 2015), or even conflicting (e.g. Pimm et al 2019) interpretations in the ecological literature.

Zoom In Zoom Out Reset image size

When we study ecological resilience, we typically assume the existence of a stability landscape with valleys and hilltops (figure 2). Such ball-and-cup stability landscapes are good approximations for understanding resilience concepts (Beisner et al 2003). Valley bottoms correspond to stable states (where a ball, representing the current system state, will eventually end up). A landscape with several valleys reflects a system with several alternative stable states. Hilltops mark the unstable states (or saddle points) that set the borders between the different valleys.

Zoom In Zoom Out Reset image size

Shifts from one valley to another can take place in two ways. First, when a disturbance pushes the ball beyond the neighboring hilltop (figure 2). A disturbance reflects changes in the state variables of the system, such as a sudden loss of a species biomass due to a mass mortality event like a fire, or an increase of a population abundance due to migration. Second, when changes in conditions alter the shape of the stability landscape itself by shrinking the current valley to disappearance so that the alternative valley becomes the only possible state. In that case, changes in conditions reflect changes in the environment or interactions that affect processes between state variables, such as an increase in competition strength between species, or an increase in habitat fragmentation. Shifts between valleys are referred to as 'regime shifts', 'catastrophic shifts', or 'tipping' events (Scheffer et al 2001, Lenton 2013, Petraitis 2013).

The balls-and-cups representation is a useful concept to apprehend the idea of the stability landscape and its link to ecological resilience (Beisner et al 2003). It is clear that the geometrical properties of the stability landscape are defining the basic properties of resilience (figure 2): (a) the size (area or volume), (b) the depth, and (c) the slopes of the each basin of attraction. Altogether these properties define the ability of the system to withstand disturbances without shifting to a different state (i.e. its ecological resilience).

The problem is that, for most systems, we do not know the shape of the stability landscape. In principle, we could derive their shapes mathematically. A useful summary of different mathematical descriptions of a stability landscape can be found in Pawlowski (2006). The most used mathematical description of stability landscapes is given by the potential function. The potential function describes the potential energy of the system for different parameter values (Strogatz 1994), where minima and maxima respectively correspond to stable and unstable equilibria, while the slopes of the potential surface are proportional to the rates of change of the system. In other words, the potential function describes a gradient surface that determines the direction and the speed at which the system (the ball) is moving across and within basins of attraction.

However, deriving the potential of a given system requires knowledge about the actual underlying dynamical model that describes the study system. Things get even more challenging, because even if 'a' model is assumed, it does not necessarily mean that a potential function exists for that model. Think about the well-known Lotka–Volterra two species predator–prey model that leads to a limit cycle behavior. The potential function of such a model should reflect a stability landscape where the system orbits but at the same time converges to the valley bottom of the stability landscape—such landscape would correspond to an impossible object (Rodríguez-Sánchez et al 2020). In general, it will be difficult to find a potential function for systems with more than one dimensions (Rodríguez-Sánchez et al 2020), although there have been ways to approximate stability landscapes for high dimensional systems (Zhou et al 2012) and see Nolting and Abbott (2016) for a very instructive discussion for ecologists.

In what follows, we assume that a generic potential function exists for our system of interest, and we explore how the properties of the stability landscape can be (mathematically) linked to metrics of ecological and engineering resilience. More precise mathematical descriptions of most of the properties we present can be also found in Meyer (2016), Krakovská et al (2021).

4.1. A (generic) potential function

Let's imagine a dynamical system described by the function f(x). The potential function U(x) of such dynamical system is defined as (Strogatz 1994):

$\begin{align} U(x) = -\int_{x_0}^{x}f(x)dx \end{align} \tag{ 1 }$

where f(x) is a set of ordinary differential equations that describe our ecological system:

$\begin{align} \frac{dx}{dt} = f(x). \end{align} \tag{ 2 }$

For simplicity, we will use one-dimensional systems for which a potential function normally exists (it can be either estimated in a closed form or integrated numerically as long it is smooth) (figure 3(B)). In higher-dimensional systems (i.e. with more than one dimensions), the conditions for a gradient potential to exist become extremely difficult (Rodríguez-Sánchez et al 2020).

Zoom In Zoom Out Reset image size

The potential functions U(x) of well-known low-dimensional ecological models used in the ecological literature to study resilience can be derived analytically, are smooth, and many of them can be approximated by a 4th degree polynomial function (Saade et al in prep; appendix A.1). In fact, a 4th-degree polynomial is one of the simplest function that can describe the potential stability landscape of a one-dimensional system with two possible wells (Strogatz 1994). With this in mind, we write a generic 4th-order polynomial function to describe the potential U(x) of our hypothetical ecological system as:

$\begin{align} U(x) = -\left(\frac{1}{4} \alpha_3 x^4 + \frac{1}{3} \alpha_2 x^3 + \frac{1}{2} \alpha_1 x^{\,2} + \alpha_0 x\right) \end{align} \tag{ 3 }$

where α_i ($i = 0,1,2,3$) are the polynomial coefficients that determine the shape of the potential. By tuning these coefficients, we can describe a stability landscape with two basins of attraction. For example, when $a_0 = 0, a_1 = -5, a_2 = 0, a_3 = 1$, the two alternative stable equilibria (attractors) are x₁ and x₂ whereas the hilltop x_un corresponds to the unstable equilibrium (saddle) that marks the border between the two basins of attraction (figure 3(A)).

In this theoretical potential, we can explore how the size and shape of the two basins of the potential landscape are related to commonly used resilience metrics.

4.2. Resilience metrics

From the previous section, it is clear that metrics of ecological and engineering resilience are strongly linked as some are functions of the others. For example, mean exit time depends on the potential depth, whereas return time, variance and autocorrelation are all functions of the potential's curvature. At the same time, some relationships commonly thought as given are not necessarily true. For instance, the probability of escaping from a given valley to the alternative one (which is approximated by the mean exit time) depends on the potential depth and not on the distance to the basin threshold as it is commonly thought.

To confirm these expected relationships (and explore unexpected ones), we estimated correlations between ecological and engineering resilience metrics in the generic one-dimensional potentials described by the fourth degree polynomial (equation (3)). In addition, we estimated correlations for a potential described by a more versatile exponential function (equation (B.2)), as well as the potentials of three classical ecological models exhibiting bistability (appendix A.1). A detailed description of the functions and how we performed the correlations can be found in appendices A.2 and B.

We found strong correlations between all metrics for the polynomial function (figure 4(A)). All engineering resilience metrics (i.e. curvature, variance, autocorrelation) were perfectly correlated with each other (=1) as they are all functions of the curvature (equations (7) and (8)). A strong correlation was also found between the ecological resilience metrics of potential depth and distance to threshold. Similarly, the potential depth was perfectly correlated to the mean exit time. The most interesting result was the strong correlation between the potential depth and the engineering resilience metrics. Although these relationships might be a consequence of the smoothness and symmetry of the polynomial function, we found similarly positive (albeit weaker) correlations for the exponential function especially in the correlations between potential depth and engineering metrics (figure 4(B)), while the correlations between distance to threshold and engineering resilience metrics break down in that case.

Zoom In Zoom Out Reset image size

When we estimated these correlations in three classical bistable ecological models, we found that correlations between engineering and ecological resilience metrics were positive especially between potential depth and engineering resilience, except for the desertification model in which these correlations disappear (figures 4(C)–(E)). Investigating this more in detail revealed that the correlations were different in strength depending on which of the two alternative equilibria was considered. When we looked only at the high equilibrium (the usually desired state), correlations were mostly positive between engineering and ecological resilience metrics, also in the desertification model (figure A.3). Interestingly, in all three models, we found opposite results to the polynomial and exponential functions for correlations between mean exit time and engineering resilience metrics, which seems to be mainly due to what happens in the low equilibrium (figure A.3). Altogether, these results suggest that engineering resilience (local stability) metrics can be used to approximate ecological resilience in some specific simple settings, yet, future work should explore in more details when and why these correlations break down. For example, Nolting and Abbott (2016) have shown that, in a two-dimensional potential described by a double exponential function, resilience metrics such as potential depth and curvature can be decoupled.

We revisited a bibliographic analysis of stability metrics in the ecological literature between 1950 and 2018 (Kéfi et al 2019) to check which metrics were used to quantify either engineering or ecological resilience. We classified metrics such as dominant eigenvalue, whether the system recovers or not, return time, coefficient of variation, variance, standard deviation, and reactivity as metrics of engineering resilience (local stability) (figure 5(C)). Metrics such as type of attractor, distance to bifurcation, skewness, detection of a regime shift or the quasi-potential of the system were categorized as measures of ecological resilience (non-local stability) (figure 5(D)). For more details about the bibliographic analysis refer to Kéfi et al (2019).

Zoom In Zoom Out Reset image size

We found that, out of the 459 papers studied, 223 estimated only engineering resilience metrics (49%), 51 estimated only ecological resilience metrics (11%), and only 18 evaluated both (4%) (figure 5(A)). More strikingly, out of the 805 measures of resilience metrics across the 459 papers of the database, engineering resilience was measured over four times more often than ecological resilience (366 vs 80, figure 5(B)). This difference is even higher if we only consider metrics measured in the field: engineering resilience metrics have been measured empirically 14 times more than ecological resilience metrics. This literature analysis highlights that there is a well-defined set of metrics that allows evaluating engineering resilience, whereas ecological resilience seems to be much more difficult to estimate, especially in real systems (Schroder et al 2005).

Below, we focus on examples of a number of methods that have been used theoretically and practically to estimate resilience metrics based on the properties of the stability landscape (table 1).

Most of these methods follow two approaches: (a) neglect the complexity and approximate an ecological system with a (low-dimensional) model, parameterized in a way that quantitatively matches the empirical system as much as possible, and then use the model to estimate ecological resilience; (b) Estimate resilience metrics that can be directly estimated from empirical data.

6.1. Based on a model

6.2. Based on data

Engineering and ecological resilience have largely been seen as alternative 'worldviews' or 'paradigms' in ecology since the 1970s. While engineering resilience refers to local stability properties, ecological resilience relates more to non-local stability. Engineering resilience is generally easier to measure since it requires information around the present equilibrium state. Instead, ecological resilience is much more difficult to quantify since it requires non-local information way beyond the present equilibrium state. For this reason, these two types of resilience have been mostly studied separately from each other, and only few studies have combined them (figure 5). Still, in the present context of increasing pressures and multiplicity of disturbances on natural systems, the need to operationalize ecological resilience for assessing and managing ecological systems is greater than ever. However, operationalizing ecological resilience has been a difficult, if not chimeric, task. In our view, asking whether such operationalization is possible requires going back to the basic definitions of these two types of resilience and looking at the current state of the theory in order to find out if and how we can go further.

7.1. Ecological and engineering resilience metrics can be strongly correlated, at least in one-dimensional systems

7.2. Stability landscapes are not just metaphors

7.3. Measuring ecological resilience will always be context dependent

7.4. Measuring ecological resilience in relative terms may be easier than measuring it in absolute terms

7.5. Conclusion

We would like to thank Lucie Bourreau for her valuable work on the analysis of the ecological models.

The data that support the findings of this study are available upon reasonable request from the authors.

A.1. Potentials of the models

We here present the potentials of classical ecological models exhibiting bistability (Van Nes and Scheffer 2005, Guttal and Jayaprakash 2008).

A.1.1. Herbivory model

Noy-Meir's model (Noy-Meir 1975) describes the dynamics of grazed vegetation:

$\begin{align} \frac{dx}{dt} & = G(x)-c(x)B \nonumber\\ & = rx\left(1-\frac{x}{K}\right)-\frac{Bx}{a+x} \end{align} \tag{ A.1 }$

with x the vegetation density, G(x) the vegetation growth (which is logistic), c(x) the vegetation consumption by herbivores and B the herbivore population density. K is the carrying capacity, r the vegetation growth rate, and a the half saturation density (May 1977).

The potential of this model, obtained by integration, is:

$\begin{align} U(x) = \frac{r}{3K}x^3-\frac{r}{2}x^{\,2}+Bx-Ba\ln(a+x)+Ba\ln(a).\end{align} \tag{ A.2 }$

Note that the equation can be rewritten as follows (see also Saade et al in prep):

$\begin{align} \frac{dx}{dt} & = rx\left(1-\frac{x}{K}\right) - \frac{Bx}{a+x} \nonumber\\ & = \frac{1}{a+x} \left(-\frac{rx^3}{K} + r\left(1-\frac{A}{K}\right)x^{\,2} +(Ar-B)x \right) \nonumber\\ & = \frac{p(x)}{g(x)}. \end{align} \tag{ A.3 }$

The numerator can be rewritten as a third degree polynomial:

$\begin{align} p(x) = -\alpha_3 x^3 + \alpha_2 x^{\,2} + \alpha_1 x \end{align} \tag{ A.4 }$

with $\alpha_3\gt0$. An illustration of the potential can be seen in figure A.1 panel A.

Zoom In Zoom Out Reset image size

A.1.2. Desertification model

The nonspatial version of Klausmeier's model (Klausmeier 1999) describes the dynamics of vegetation and water in a dryland ecosystem as follows:

$\begin{align} \left\{ \begin{array}{ll} \displaystyle{\frac{dW}{dt} = A-LW-RWx^{\,2}} \\[5pt] \displaystyle{\frac{dx}{dt} = RJWx^{\,2}-Mx} \end{array} \right. \end{align} \tag{ A.5 }$

with W the water, x the plant biomass, A the water supply rate, L the water evaporation rate, R the plant uptake rate of water, J the yield of plant biomass per unit of water consumed, M the vegetation mortality rate.

Because of time scale separation (assuming water dynamics is much faster than vegetation dynamics), we can solve $\frac{dW}{dt} = 0$ and get an expression for W: $W = \frac{A}{L+RN^2}$. We thereby get a single equation:

$\begin{align} \frac{dx}{dt} = RJ\left(\frac{A}{L+Rx^{\,2}}\right)x^{\,2}-Mx. \end{align} \tag{ A.6 }$

The potential function of this model is:

$\begin{align} U(x) = -JA\left(x-\frac{\tan^{-1}\left(\sqrt{\frac{R}{L}}x\right)}{\sqrt{\frac{R}{L}}}\right)+\frac{M}{2}x^{\,2}. \end{align} \tag{ A.7 }$

Note that, in the same vein as for Noy-Meir's model, the equations of this model can be rewritten as:

$\begin{align} \frac{dx}{dt} & = RJ\left(\frac{A}{L+Rx^{\,2}}\right)x^{\,2} - Mx \nonumber \\ & = \frac{-MRx^3+RJAx^{\,2}-MxL}{Rx^{\,2}+L} \nonumber \\ & = \frac{p(x)}{g(x)} \end{align} \tag{ A.8 }$

where the numerator is a third degree polynomial:

$\begin{align} p(x) = -\alpha_3 x^3 + \alpha_2 x^{\,2} - \alpha_1 x \end{align} \tag{ A.9 }$

with $\alpha_i\gt0$. An illustration of the potential can be seen in figure A.1 panel B.

A.1.3. Lake eutrophication model

Carpenter's model (Carpenter et al 1999) is a lake eutrophication model:

$\begin{align} \frac{dx}{dt} = a-bx+r\frac{x^k}{x^k+h^k} \end{align} \tag{ A.10 }$

with x the amount of phosphorus (P) in the water column, a the rate of P input from the watershed, b the rate of P loss, r the maximum rate of recycling of P. The overall recycling rate is assumed to be a sigmoid function of P for which the exponent k controls the steepness of the curve at the inflexion point and h is the recycling half-saturation rate.

The potential function of this model is (Pawlowski 2006) (for k = 2):

$\begin{align} U(x) = \frac{b}{2}x^{\,2}+hr\tan^{-1}\left(\frac{x}{h}\right)-(a+r)x.\end{align} \tag{ A.11 }$

An illustration of the potential can be seen in figure A.1 panel C.

A.2. Correlations between resilience metrics in the three ecological models

We looked at the correlations between the six resilience metrics in the three ecological models (figures 4(C)–(E)). We did this by systematic combinations of each model parameter. Specifically, for the Herbivory model we used $\alpha = [0.1,5]$ at 0.1 intervals, $r = [0.1,3]$ at 0.05 intervals, $B = [0.05,1.5]$ at 0.1 intervals, $K = [0.5,10]$ at 0.5 intervals. For the Desertification model we used $R = [0.1,2]$ at 0.1 intervals, $J = [0.5,10]$ at 0.5 intervals, $A = [0.1,2]$ at 0.1 intervals, $L = [0.1,2]$ at 0.1 intervals, $M = [0.05,1]$ at 0.05 intervals. For the Eutrophication model we used $\alpha = [0.05,1.5]$ at 0.05 intervals, $b = [0.05,2]$ at 0.05 intervals, $r = [0.5,8]$ at 0.5 intervals, $h = [0.05,2]$ at 0.05 intervals.

The three engineering metrics are very strongly correlated to each other in the three models. In the same vein, distance to threshold and potential depth (two ecological resilience metrics) are strongly correlated to each other in the three models. They are also both positively (but not as strongly) correlated to the engineering metrics in the herbivory and in the eutrophication model but very weakly, even negatively correlated to them in the desertification model. One metrics behaves differently, exit time, which is negatively correlated to the engineering metrics in the three models and positively correlated to distance to threshold and potential depth. In sum, in this model, the engineering metrics are positively and strongly correlated to each other; the ecological resilience metrics as well. Correlations between engineering and ecological resilience metrics vary in strength and sign depending on the model.

We also looked at the resilience metrics obtained in these three ecological models for the low and high equilibria (figure A.2). For the three models, the results found for the different metrics are not always consistent. Some of the metrics indicate that one of the equilibrium is more stable than the other and other metrics suggesting the opposite. For other metrics, there is a high variability (e.g. for curvature in the herbivory model; figure A.2(A)) and comparative values can be positive but also strongly negative, suggesting that the high equilibrium is sometimes more stable and sometimes less stable than the low one. Moreover, the differences among metrics for the two equilibria are not always similar across models. For this reason, we displayed correlations among metrics by looking separately at the high and low equilibria (figure A.2).

Zoom In Zoom Out Reset image size

We found that the three engineering resilience metrics (i.e. variance, autocorrelation and curvature) are always strongly positively correlated to each other. This is also true for potential depth and distance to threshold, and this, independent of the equilibrium considered (figure A.3(E)). These two groups of metrics are moreover positively correlated with each other for all models and for both equilibria, except for the high equilibrium of the desertification model (figure A.3(E)). In that latter case, the engineering resilience metrics are only weakly correlated with the ecological resilience metrics. As already noticed earlier, one metric behaves differently than the others for the low equilibrium: the exit time (figures A.3(A)–(C)).

Zoom In Zoom Out Reset image size

We estimated correlations between ecological and engineering resilience metrics based on one-dimensional potentials described by a fourth degree polynomial (equation (B.1)) and an an exponential function (equation (B.2)) (figure B.1).

$\begin{align} & polynomial: U_p(x)\nonumber\\[2pt] & \;\; = -\left(\frac{1}{4} \alpha_3(sx)^4 + \frac{1}{3} \alpha_2(sx)^3 +\frac{1}{2} \alpha_1(sx)^2 + \alpha_0(sx)\right)\end{align} \tag{ B.1 }$

$\begin{align} exponential: U_e(x)& = 1-\alpha exp(-\beta x^{\,2})\nonumber\\ & \quad -exp(-\delta(x-\gamma)^2).\end{align} \tag{ B.2 }$

Zoom In Zoom Out Reset image size

In the polynomial function we introduced the parameter s that shrinks and expands the potential without changing the other properties of its shape. As the polynomial function is not that flexible in changing shapes (and both valleys are completely symmetrical), we also used the exponential function (equation (B.2)), where one can change the position of the two valley bottoms independently with β and delta respectively for each valley bottom. γ changes the distance between the two valleys, and alpha changes the relative depth between the two valleys.

For both functions, we systematically combined parameter values and we estimated Spearmann ρ rank correlations between metrics measured only when two basins of attraction existed (otherwise some ecological resilience metrics, like potential depth, are non-existent). In particular, we estimated correlations between potential depth (Pot_d), distance to basin threshold (Dist_thres), curvature of potential (λ_attr ) (as the dominant eigenvalue to measure recovery rate) and mean exit time ($\tau_x^{x_{un}}$), variance (V ar), and autocorrelation (Corr). We did not measure the size of the basin of attraction because as we assumed x to be unbounded, our two valley are both ending to infinity and are thus equal in size (i.e. infinite size). Also we did not measure reactivity (as it is defined only for more than one-dimensional systems). We present results for each basin of attraction separately (figure B.2) and combined (figure 4).

Zoom In Zoom Out Reset image size

In particular, we used combinations of the five parameters of the polynomial function (equation (B.1)) for the following ranges: s [1,5], $\alpha{_0} = [-5,5]$, $\alpha{_1} = [-5,5]$, $\alpha{_2} = [-5,5]$, $\alpha{_3} = [1,5]$. Each range was sampled at ten equal intervals. For variance and autocorrelation we used noise of σ = 0.25. All estimations were performed for $x = [-5.5,5.5]$ at 1000 interval points. These combinations resulted in 10 000 potential landscapes, out of which 18 074 corresponded to potential landscapes with two basins of attraction. We used these 18 074 potentials to estimate the six resilience metrics described above.

We used combinations of the four parameters of the exponential function (equation (B.1)) for the following ranges: α = [0.01,5], $\beta = [0.01,3]$, $\gamma = [0.1,4]$, $\delta = [0.01,3]$. Each range was sampled at 20 equal intervals. For variance and autocorrelation we used noise of σ = 0.25. All estimations were performed for $x = [-2,6]$ at 1000 interval points. These combinations resulted in 160 000 potential landscapes, out of which 88 201 corresponded to potential landscapes with two basins of attraction. We used these 88 201 potentials to estimate six resilience metrics described above.

All calculations were performed in R.