Ecological resilience: what to measure and how

The question of what and how to measure ecological resilience has been troubling ecologists since Holling 1973s seminal paper in which he defined resilience as the ability of a system to withstand perturbations without shifting to a different state. This definition moved the focus from studying the local stability of a single attractor to which a system always converges, to the idea that a system may converge to different states when perturbed. These two concepts have later on led to the definitions of engineering (local stability) vs ecological (non-local stability) resilience metrics. While engineering resilience is associated to clear metrics, measuring ecological resilience has remained elusive. As a result, the two notions have been studied largely independently from one another and although several attempts have been devoted to mapping them together in some kind of a coherent framework, the extent to which they overlap or complement each other in quantifying the resilience of a system is not yet fully understood. In this perspective, we focus on metrics that quantify resilience following Holling’s definition based on the concept of the stability landscape. We explore the relationships between different engineering and ecological resilience metrics derived from bistable systems and show that, for low dimensional ecological models, the correlation between engineering and ecological resilience can be high. We also review current approaches for measuring resilience from models and data, and we outline challenges which, if answered, could help us make progress toward a more reliable quantification of resilience in practice.


Introduction
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 , 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 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 ecologytwo 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 socioecological 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 . 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, Cumming et al 2005, Thrush et al 2009, Folke 2016. Instead, we focus on system responsesnot 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.

Defining resilience: the two paradigms of local vs non-local stability
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

The stability landscape for mapping resilience
When we study ecological resilience, we typically assume the existence of a stability landscape with valleys and hilltops (figure 2). Such ball-andcup 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 Figure 1. Local vs global stability: a locally stable equilibrium can also be globally stable when the system always converges to the same equilibrium (illustrated with the black point) regardless from its starting point (the black line reflects an example of a trajectory in state-space) (a). In contrast, an equilibrium can be locally stable but globally unstable if only the trajectories starting from a subset of starting points converge to it, while other trajectories converge toward another state (b). In that case, there are more than one possible stable equilibria (two in this example). Local stability thus concerns what happens in the close neighborhood of the equilibrium. In contrast, global stability is concerned with what happens in the whole state-space. Here, we use non-local stability to refer to what happens beyond the close neighborhood of the equilibrium (while not necessarily investigating the whole state-space). 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.
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, Figure 2. A (hypothetical) stability landscape of a two-dimensional system with hilltops and valleys, also known as a marble-in-a-cup or balls-and-cups landscape. Black balls are found in the bottom of the valley and represent stable states (if disturbed the ball will 'gravitate' toward the bottom of the valley). White balls are found on the hilltops and represent unstable states (if disturbed they will move away from their hilltop position). See text for more details. 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)

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): where f (x) is a set of ordinary differential equations that describe our ecological system: 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).
The potential functions U(x) of well-known lowdimensional 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: 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 1 and x 2 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.

Resilience metrics 4.2.1. Resilience to perturbations on state variables
The size and shape of the potential can be related to properties of both engineering and ecological resilience. Engineering resilience provides information about the response of the system around equilibrium to small perturbations on state variables (figure 2). Ecological resilience refers to how the system would respond to non-small perturbations, so that if the stability landscape has several wells one would need to evaluate the overall risk of shifting from one attractor to the other. Below we classify resilience metrics based on this distinction and explain their relationships to the properties of the potential.

Metrics of engineering resilience (local stability)
• Curvature of potential. In general, the absolute value of the potential curvature measures the speed of convergence or divergence after a small perturbation. This can be quantified at any point along the potential. Typically, though, it is the curvature measured at equilibrium (figure 3(B)) quantified by the dominant eigenvalue λ of the Jacobian matrix of the ecological model (equation (2)) that determines the local stability of the equilibrium: Different information can be derived from the quantification of the potential curvature. First, the sign of the eigenvalue λ qualifies whether the state is stable (λ < 0) or not (λ > 0). Second, the magnitude of the eigenvalue λ determines how fast the system returns back to the disturbed equilibrium after a disturbance: Return (or recovery) time (also measured as return (or recovery) rate (1/τ return ) back to the desired attractor is defined as engineering resilience (Pimm 1984). Note that for large deviations from equilibrium, the eigenvalue approximation fails. In that case the response of the system back to equilibrium is determined by the rate of change along the potential (figure 3(B)). The above defined return time defines the long-term (or asymptotic) resilience (Arnoldi et al 2018). It is possible that after a disturbance the system might initially move away from the attractor before eventually returning back to it. This property that is also related to the curvature of the potential is called reactivity (Neubert and Caswell 1997) but it is defined only for higher than 1dimensional systems: with Csym = C + C T where C is the Jacobian matrix and λ xattr the dominant eigenvalue of C at equilibrium x attr . • Variance. Under small, random and continuous disturbances, the system moves in a way that is defined by the curvature of the potential around the attractor (figure 3(B)). One way to think about this is to imagine the ball to wiggle around the bottom of the valley while still staying within its basin of attraction. For small disturbances, the resulting stationary distribution will be described by its variance V ar(x) that depends on the strength of stochasticity σ and the dominant eigenvalue λ (Ives et al 2003): The flatter the bottom of the valley is, the smaller the absolute value of λ will be and consequently the higher the variance. Note that as this definition considers small levels of noise, it depends on the eigenvalue at the local equilibrium and it does not take into account the whole curvature of the potential. In that sense it is a local measure of resilience (figure 3(C)). • Autocorrelation. In addition to variance, the values of the resulting stationary distribution will be correlated in time Corr(x t , x t+1 ) defined by the eigenvalue λ: A flat valley bottom (i.e. small absolute value of λ) will lead to high values of autocorrelation. Note that as this definition (like the above for variance) considers small levels of noise, it does not take into account the whole curvature of the potential so that it is a local measure of resilience.

Metrics of ecological resilience (non-local stability)
• Potential depth. This metric reflects the amount of 'effort' needed to climb against the gradient of the potential in order to pass over the saddle point (hilltop). More precisely, it is the relative depth of the bottom of a valley to the hilltop (figure 3(A)).
where x attr is one of the alternative attractors (stable equilibrium) and x saddle is the saddle point (unstable equilibrium). This metric reflects how 'easy' it is for a system to be pushed away from its current state and shift to another basin of attraction. This metric has been also described as 'resistance' by , and it resembles the widely used notion of resistance defined as the ability of a system to persist despite a disturbance (Grimm et al 1997), although resistance is commonly quantified as the relative amount of system change given a specific size of disturbance (Ingrisch and Bahn 2018). • (Minimum) Distance to the basin threshold. This metric describes the distance between the valley bottom (stable equilibrium) and the limit of the basin of attraction (saddle point). Another way to think about this metric is the minimum amount of a single perturbation to the system state that one has to apply so that the systems shifts to the other state (figure 3(A)).
If the system is not at the bottom of the valley (at equilibrium) then the distance between the current state of the system and the saddle x un has been called by  precariousness. • Size of the basin of attraction. The size of the basin of attraction reflects the likelihood of the system to be in a given state after any random-non-smallperturbation (figure 3(A)). Mathematically, it can be defined as the set of initial conditions in state space that lead to a particular state. For example, for a one-dimensional system, the set W would be defined within a range of x values from x A to x B as: This set can be described by different geometries depending on the dimensionality of the system. For one-dimensional stability landscapes, it is simply the width of the basin of attraction (e.g. from equation (6) |x A − x saddle |). For two-dimensional stability landscapes, the size is described as an area, while for higher dimensions it is quantified as a (hyper-)volume (Menck et al 2013).
Although the size of the basin of attraction is a complete measure of the stability domain of a given attractor, sometimes it might be more relevant to use other metrics than size. For example, in a onedimensional stability landscape the width of the basin of attraction is not a useful metric when there are only two basins of attraction, as each basin has a threshold only on one side (figure 3(A)). In that case, it is more relevant to compare basins based on the minimum distance to the basin threshold. Similarly, in a two-dimensional stability landscape it might be that it is the minimum distance between the two basin thresholds (also referred as latitude (Walker and Meyers 2004)) (figure 2) that describes better the likelihood of the system to stay within a basin of attraction after a perturbation.
• Likelihood of bi-(multi)-modality. The existence of bi-(multi)-modality is a proxy of bistability and potential flipping between alternative states. Bi-(multi-)modality can be derived from the peaks of the probability distribution of the system state. The probability distribution describes the probability to be at any point along the stability landscape in the presence of (environmental or demographic) noise (figure 3(C)). The probability distribution depends on the potential function and the strength of noise.
The simplest way to imagine this is to consider the ball (system state) being 'pushed' around within the basin of attraction due to small and continuous perturbations. Where the system will spend more time will depend on the shape of potential landscape and how noise is acting on the system. The values of system state where the likelihood is the maximum are the modes of the distribution and correspond to the alternative attractors of the system. The probability P(x) of the system state being at a point x is given by: where N is a normalization factor, σ is the level of noise. This expression is the solution to the Fokker-Planck stochastic version of equation (1) (Gardiner 2003). If noise has an additive structure (e.g. Gaussian distribution added on the system state), the probability distribution is of the form (equation (12)). For more complex structures of noise (i.e. multiplicative), the probability distribution becomes a function of noise as well (see for instance Guttal and Jayaprakash 2008). • (Mean) Exit time to the alternative basin of attraction. The (mean) exit time τ is the average amount of time that it will take for small and continuous disturbances to 'push' the ball out of its current basin of attraction. This time depends on potential properties (potential depth and curvature potential) as well as on the strength of noise (Nolting and Abbott 2016): where λ x quantifies the curvature of the potential at point x (see also below under 'Curvature').
Obviously, if stochasticity is low (σ → 0) exit time becomes infinite, in other words it is impossible to escape from the current basin of attraction. • Skewness. In the presence of continuous, random disturbances that 'push' the system around the basin of attraction (but not outside of it), the resulting stationary distribution of the system state will reflect the shape of the basin of attraction. If the basin of attraction is asymmetric, the stationary distribution will also be asymmetric (figure 3(C)). Asymmetry suggests the existence of another attractor but it is not a proof of it. The simplest way to measure the basin's asymmetry is to quantify the skewness γ of the resulting stationary distribution (Guttal and Jayaprakash 2008): where µ is the mean value of the probability distribution P(x). Negative skewness would imply leftside asymmetry, whereas positive skewness would imply asymmetry to the right.

Resilience to perturbations on parameters
The above metrics reflect properties of a nonchanging (static) stability landscape. However, the stability landscape can also change due to changes in external (e.g. environmental) conditions (Suding et al 2004). Such changes are related to changes in system parameters. For instance, this would mean that the coefficients of the polynomial potential function (equation (3)) are changing. In such context of perturbations in parameter space, other metrics of resilience become relevant. One of them is the amount of change in an external parameter that is needed for the current attractor to disappear and be replaced by its alternative attractor (figure 3(D)). Mathematically, the appearance and disappearance of attractors as well as changes in their stability correspond to bifurcation points (Kuznetsov 1995). It is straightforward to imagine that such bifurcations can be caused by more than one external parameters. In that case, there is actually a parameter space (analogous to the statespace that defines the 'Size of Basin of Attraction') that defines the size of parameter changes that can be tolerated without leading to shifts in the number or in the stability of attractors. This parameter space defines the 'structural stability' of a system (Smale 1967, Thom 1989).
There are different ways for quantifying resilience to perturbations in parameter space.
• Distance to bifurcation. Analogous to the 'Distance to basin threshold' metric, we can find the minimum amount of change needed along one parameter p to incur the disappearance of the current valley bottom (figure 3(D)).
where p attr is the current value of the parameter and p bif is the value of the parameter at which one of the bifurcations occurs. • Size of hysteresis. Hysteresis defines the range of parameter space for which bistability exists along a single parameter p (figure 3(D)).
where p bifA and p bifB are the values of the parameter at which the bifurcations occurs.

Relationships between engineering (local) and ecological (non-local) resilience metrics derived from generic one-dimensional potentials
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.
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  (Carpenter et al 1999). Correlations are based on metrics estimated for both basins of attraction (18 074 cases for the polynomial function, 88 201 cases for the exponential function, 90 257 cases for the herbivory model, 76 515 cases for the desertification model and 16 145 cases for the eutrophication model). For the three last metrics a sign correction was done so that the metric increases as stability increases. Red colors represent positive correlations, and blue colors negative correlations. Values in the circles are the correlations measured. All correlations are significant (p < 0.05). potential described by a double exponential function, resilience metrics such as potential depth and curvature can be decoupled.

How to measure resilience in practice
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).
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 , whether the system recovers or not after a pulse perturbation), cv: coefficient of variation, var: variance, acorr: autocorrelation, sdev: standard deviation, react: reactivity, attract: type of attractor, dist bif: distance to bifurcation, skew: skewness, shift: detection of a regime shift, qpot: quasi-potential).
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.

Based on a model
Being able to derive (analytically or numerically) the potential function of the ecological model enables us to estimate all resilience metrics. Recall that this is usually possible for one-dimensional systems.
Similarly, numerical tools of bifurcation analysis enable us to explore the presence of thresholds, the likelihood of bistability and the range of hysteresis. For higher-dimensional systems, it can be possible to approximate the potential by the quasipotential (Nolting and Abbott 2016) (with varying degree of error depending on the model (Rodríguez-Sánchez et al 2020)). Recent works provide useful numerical techniques and brilliantly illustrate how the derived quasipotential can be used to quantify ecological resilience in ecological models (Nolting and Abbott 2016).
However, even after having defined a model, finding the potential function or performing a bifurcation analysis to explore the range of hysteresis might still be difficult. In that case, numerical methods can indirectly quantify some of the resilience metrics. A classical approach is to randomly sample initial conditions and record the number of alternative endpoint attractors (Beddington et al 1976). For instance, the  (2016) probability of reaching a given state relative to the total amount of trials is a proxy of the size of the basin of attraction of that state (Lundström 2017  2015)), and have not been widely adopted by ecologists especially in highdimensional systems, like multispecies communities.

Based on data
If we have no model at hand but lots of observations of the study system, we can reconstruct a 'probabilistic' stability landscape by assuming that the system has 'visited' most parts of the 'true' stability landscape. Using the probability distribution of system states could be achieved by simple histograms or by fitting uni-or multivariate density (e.g. Gaussian) functions to the data (Hirota et al 2011, Berdugo et al 2017. In this way, the dominant modes are seen as the potential attractors and resilience can be evaluated based on their number and distance. We could even go a step further and fit a generic function (like for instance a polynomial or exponential function appendix B) as approximation of the potential function. This approach is followed in potential analysis (Livina et al 2010) that can be used not only to describe bi-(multi)modality, but also describe potential properties. For example, based on standard model selection tools, the best-fitted polynomial function could be used to identify the number (polynomial degree), position and distance between alternative basins. Although these approaches appear promising, they still are limited to low dimensional application and require a lot of observations; that is why the best examples of the application of these approaches come from remote sensing data describing vegetation cover in forests, savannas and drylands (Hirota et al 2011, Berdugo et al 2017, van Belzen et al 2017. Yet, in multi-dimensional cases (e.g. multispecies communities), aggregate measures of community/ compositional state can be used through ordination techniques (such as PCA) for identifying alternative states and attraction basins through clustering based on densities (Génin et al 2020), or through building weighted networks of states calculated by pairwise maximum entropy models (Suzuki et al 2021).
Other methods attempt to reconstruct the parameter space where bistability exists by fitting the simplest ('normal') form function (similar to equation (3)) that produce cusp surfaces (Jones 1975, Grasman et al 2009 in three dimensional space made of the data itself sampled at different environmental conditions (Sguotti et al 2018).
Regime shift detection has been perhaps the mostused approach for inferring the likelihood of bimodality (i.e. the existence of alternative states). In that sense regime shift detection can been seen as a proxy of ecological resilience. However, the existence of a regime shift is not necessary proof of bi-(multi)stability (Petraitis 2013). Regime shifts can be caused by ways that are not related to loss of ecological resilience or the presence of alternative states (Dakos et al 2015, Ratajczak et al 2018. Thus, the detection of shifts in a dataset (in a timeseries or along a gradient) is only an indication that there could be underlying bistability. Regime shift detection is typically done using breakpoint analysis (Andersen et al 2009). There is a wealth of methods and statistical software for conducting breakpoint analysis based either on statistical tests (Beaulieu et al 2012), fitting threshold-general additive models (Ciannelli et al 2004) or threshold auto-regressive models to the data (Gröger et al 2011, Ives and Dakos 2012, Laitinen et al 2021. A recent approach is based on estimating changes in the pattern of system dynamics that can be used as proxies of distance to bifurcation (ecological resilience). In that sense, it is not the magnitude per se but the relative change in time, space, or along a gradient that is used as a proxy of loss (or gain) of resilience. These indicators (also referred to as 'early-warning signals' (Scheffer 2009) correspond to trends in variance, autocorrelation and skewness. Although they are often interpreted as measures of ecological resilience, they are primarily based on metrics of engineering resilience. As shown in figures 3(E) and (F), this interpretation is driven by the fact that changes in the overall (non-local) properties of the potential (due to the decreasing distance to bifurcation) are strongly correlated to changes in the local stability properties of the potential ( figure 4). Indeed, these changes are not specific to bistable potential landscapes, but signal proximity to any bifurcation at which the (local) resilience of the present state changes (Kéfi et al 2013).

Discussion
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.

Ecological and engineering resilience metrics can be strongly correlated, at least in one-dimensional systems
Using two generic functions to define the potential that mathematically describes a system's stability landscape, we estimated engineering and ecological resilience metrics and explored their correlations. We find that engineering and ecological resilience metrics can be positively correlated with each other (figures 4(A) and (B)). This result also seems to hold to some extent in classical one-dimensional ecological models (figures 4(C)-(E)). At first sight, this might be expected given that most of the resilience metrics we explored (but not all) are functions of each other. It also suggests that perhaps we can use at least some metrics of engineering resilience as proxies for some metrics of ecological resilience and the other way around in these simple models. However, this observation remains somewhat surprising because there is no theoretical reason that local stability should be related to non-local stability. Actually, we can easily imagine ways where the properties of the stability landscape can change independently, so that a steep basin of attraction (high curvature, high engineering resilience) should also be shallow (small potential depth, low ecological resilience) or narrow (small distance to threshold, low ecological resilience).
One question raised by these results is the extent to which our understanding of the relation between engineering and ecological resilience might be strongly biased by the type of models used in ecology to study ecological resilience. The low dimensional models (even of only a single dimension like the ones we presented here) are missing realism when approximating a much more complex reality, but they still provide some insight to otherwise intractable problems. For example, strong correlations between diverse stability metrics have been shown to exist in multispecies ecological communities as well (up to a 100 dimensions) (Dominguez-Garcia et al 2019), although the authors did not focus on all engineering and ecological resilience metrics studied here. When assuming that a linear and monotonic relationship exists between the amplitude of disturbances and their effects on the system state, engineering and ecological resilience become not only correlated but even equivalent (Zampieri 2021). Therefore, more generally, the question is whether and under which conditions we should expect strong correlations between ecological and engineering resilience metrics. For instance, we can construct peculiar one-dimensional models where the potential does not change smoothly but discretely (Menck et al 2013), or where the potential becomes steeper rather than shallower when approaching a bifurcation (Titus and Watson 2020). Even experimentally, it has been shown that in yeast populations ecological resilience (quantified as distance to the basin threshold) can be negatively related to engineering resilience (quantified as recovery rate), under multiple types of pressures that act simultaneously in different directions (Dai et al 2015). We can also construct two-dimensional potentials whose shape properties can change independently (Nolting and Abbott 2016). However we still do not know to what extent this independence affects the correlations between the resilience metrics that we find for onedimensional potentials. In higher dimensions (more than 2), possibly this effect becomes stronger but the problem is that we cannot easily explore this question since, in high-dimensional models, the potential often cannot be derived (or no potential exists) (Hastings and Wysham 2010). New approaches to construct a quasipotential (Nolting and Abbott 2016, Rodríguez-Sánchez et al 2020) in more than onedimensional systems might help us to systematically test to what extent the strong relationships between ecological and engineering resilience still hold.

Stability landscapes are not just metaphors
Interestingly, studies have suggested that empirical potentials can be reconstructed if the right type of data is available, such as replicates of ecosystem data in space (e.g. Hirota et al 2011, Berdugo et al 2017. In those cases, the reconstructed potentials resemble the ones derived from the simple one-dimensional models. The reconstructed potentials are described by high-order polynomial functions that, as shown in appendix B, translate into strongly correlated resilience metrics. More research is needed to address whether there is a priori a good reason to fit a polynomial potential to empirical data. Simulated data from high-dimensional models can serve as test ground of potential reconstruction before applying these approaches on the increasingly available empirical data.

Measuring ecological resilience will always be context dependent
The mathematical definitions of resilience we presented are theoretical approximations. In practice, we would need to make assumptions and simplifications to fit our empirical examples into these definitions. We need to put boundaries and assume a closed system in order to define the dimensionality of the stability landscape. We need to distinguish perturbations on state variables (resilience in state-space) from perturbation on state parameters (resilience of statespace). We need to identify and isolate the effect of multiple perturbations. These are all challenging to solve empirically and will affect our estimates of ecological resilience. Our knowledge of the structure of perturbations will affect the choice of resilience metrics (Grumm 1976). Which variable to use for characterizing a property of resilience will depend on data availability. A state variable can be considered as a state parameter and vice versa (for instance slow changes in soil water make it an external parameter, but if soil water dynamics are fast it can be seen as state variable). These are hard choices to make and there is no single-best answer. We could, however, use the mathematical definitions to adjust the resilience metrics specifically for different ecosystem types (Yi and Jackson 2021). Or we can be flexible and combine different metrics to produce approximations of ecological resilience.

Measuring ecological resilience in relative terms may be easier than measuring it in absolute terms
At the moment, widely-used proxies of ecological resilience, are directional changes in statistical patterns of the dynamics of a system around equilibrium (Scheffer 2009). They are proxies, because they are not absolute measures of resilience, but they reflect changes in resilience compared to a control or a baseline (Dakos et al 2015). They predominantly measure changes in variability, correlations, or general changes in spatio-temporal patterns (Fath et al 2003, Mayer et al 2006, Kéfi et al 2014, Sundstrom et al 2017. Some are direct metrics of engineering resilience (e.g. variance, recovery rate, autocorrelation), and thus are related to local rather than non-local stability properties. Although it is well-understood that these changes are not specific to multistable systems but characterize changes in local stability when a bifurcation is approached (Boettiger et al 2013, Kéfi et al 2013, they still capture slow responses caused by the more shallow potential surface in the vicinity of saddles and tipping points especially when we consider strong perturbations far from equilibrium ( figure 3(B)). Together with their advantage of being straightforward to quantify in empirical data, these proxies have been used to flag approaching shifts in natural systems between alternative states (e.g. Dai et al 2012, Rindi et al 2017. In that sense, they appear to stand somewhere between ecological and engineering resilience but do not fully bridge the gap.

Conclusion
No doubt, the fact that it is easy to conceptualize ecological resilience but difficult to operationalize it has contributed to its misuse (or even abuse) as a concept. Most would agree that resilience cannot be captured by a single number , but multiple metrics should be used to provide a composite assessment (Grafton et al 2019). However, the question remains of which metrics these are. We are not the first who ask these questions aiming to bridge the gap between theory and practice (Grumm 1976, Meyer 2016, nor do we claim that we succeed in answering them. However, we believe that it is fruitful to identify what we can do (even if approximately) contrary to what is impossible to do. For instance, we need to test to which extend the correlations between metrics of ecological and engineering resilience hold in the ecological models we use. At the same time, it is time to use big data to probabilistically identify alternative states and potential thresholds (e.g. Berdugo et al 2020). In this way, we may finally bridge the gap between the metaphors of potential landscapes and basins of attractions and the quantification of ecological resilience in practice.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors. 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: Note that the equation can be rewritten as follows (see also Saade et al in prep): The numerator can be rewritten as a third degree polynomial: with α 3 > 0. An illustration of the potential can be seen in figure A.1 panel A.

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: 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 dW dt = 0 and get an expression for W: W = A L+RN 2 . We thereby get a single equation: The potential function of this model is: Note that, in the same vein as for Noy-Meir's model, the equations of this model can be rewritten as: where the numerator is a third degree polynomial: with α i > 0. 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: 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): 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 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).
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)).

Appendix B. Polynomial and exponential functions to estimate resilience metrics correlations
We estimated correlations between ecological and engineering resilience metrics based on onedimensional potentials described by a fourth degree polynomial (equation (B.1)) and an an exponential function (equation (B.2)) (figure B.1).
Correlations between the six resilience metrics for the polynomial (A) and exponential (B) functions. Left panels correspond to correlations of the low equilibrium and right panels of the high equilibrium. For the three last metrics a sign correction was done so that the metric increases as stability increases.
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 (τ xun x ), 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).
All calculations were performed in R.