Ecological systems such as bird species competing for resources (Gibbs & Grant, 1987), plants and their mutualistic pollinators (Burkle et al., 2013), and microbes interacting in the human gut (Costello et al., 2012) face constant perturbations in changing environments. The resilience of such systems to different types of perturbations has been one of the most debated concepts in ecological research and has important implications for biodiversity and ecosystem functioning (Hodgson et al., 2015;Folke et al., 2016;Pimm et al., 2019). Broadly, resilience has been referred to as the ability of an ecological system to resist and recover from external perturbations (Hodgson et al., 2015;Capdevila et al., 2020). To understand these concepts and measures, epistemological work has established three necessary conditions to study ecological responses to perturbations (Justus, 2013): (i) a description of the ecological system, which mathematically takes the form of a population dynamics model, (ii) the definition of a reference state from which the system will be perturbed, and (iii) the definition of the type and magnitude of perturbations, which can act on state variables, model parameters (i.e., the system's structure), or both (see Glossary).
Following this rationale, recovery (although it has also been called resilience, see Justus (2013) and Pimm et al. (2019)) has been typically defined as how fast a system returns to a reference state after a given perturbation (Hodgson et al., 2015;Capdevila et al., 2020). In turn, resistance has been defined as how much change a system can exhibit after a given perturbation (Hodgson et al., 2015;Capdevila et al., 2020). Hence, a sufficient condition for a system to be resilient is to exhibit fast recovery and high resistance. However, in order to assess the resilience of a system, it is paramount to have informative indicators of recovery and resistance and to understand the interconnection between these two measures, which may not necessarily be complementary (i.e., provide different qualitative information about the resilience of a system) (Domínguez-García et al., 2019;Kéfi et al., 2019).
Recovery has had a strong tradition in the mathematical analysis of asymptotic dynamical stability (hereafter dynamical stability) (Strogatz, 2001). In ecology, recovery has been typically conceptualized as the return rate along the slowest direction to a reference equilibrium state (Pimm & Lawton, 1977;Dakos et al., 2015;Donohue et al., 2016). This return rate can be quantified by the real part of the largest eigenvalue (if negative) of the Jacobian matrix when evaluated at equilibrium (Novak et al., 2016;Logofet, 2018). The Jacobian matrix represents the linearized dynamical forces acting around a given state. This implies that recovery has been mathematically quantified by the long-term return rate (or the return time if inverted) of a system to a reference equilibrium state after small perturbations (per the linear validity of the Jacobian) acting on species abundances (i.e., the state variables of the system).
By contrast, the measurement of resistance has had a weaker tradition in the field of dynamical systems (Strogatz, 2001) and it is unclear whether it can be linked to dynamical stability in ecological systems (Justus, 2013;Domínguez-García et al., 2019). Resistance has been typically quantified as the magnitude of displacement of a system following a perturbation (Carpenter et al., 1992;Hillebrand et al., 2018). Differently from recovery, resistance is not necessarily limited to a particular quantitative equilibrium state (i.e., abundance distribution at equilibrium), and it is expected to reflect the response of a system to perturbations acting on both species abundances and on the structure of a system (Donohue et al., 2016;Pimm et al., 2019). Formally, this structure can be represented by a population dynamics model describing the governing laws of a system and its parameters (Justus, 2013). That is, while species abundances represent the state variables of a system, the structure represents the model parameters (and the model itself). Thus, it remains unclear whether recovery from abundance perturbations and resistance to structural perturbations can be quantified within the same formalism, and whether there are potential connections between these two aspects of resilience (Hodgson et al., 2015;Pimm et al., 2019).
Focusing on parameter perturbations, a dynamical system is said to be structurally stable if the topology of the phase portrait (i.e., the qualitative behavior of a dynamical system) is preserved under smooth parameter changes (Smale, 1967;Arnold, 1988). For example, a system with S species can be considered structurally stable if none of the species goes extinct after small changes in the parameters corresponding to species intrinsic growth rates (Saavedra et al., 2017). Importantly, it has already been shown that the sensitivity to stochastic perturbations of abundances in the vicinity of an equilibrium is equivalent to the sensitivity to fluctuations of the parameters within the same vicinity (Arnoldi & Haegeman, 2016). Specifically, the response to external infinitesimal perturbations acting on species abundances agrees with the minimal parameter perturbation able to render the system dynamically unstable. Furthermore, structurally unstable systems (i.e., systems that change their qualitative behavior after small parameter changes) are characterized by conditions under which the largest eigenvalue of the Jacobian matrix is equal to zero (Strogatz, 2001;Duan, 2015). These previous results have established a potential direction to link dynamical indicators (i.e., measures related to perturbations acting on species abundances) with structural indicators (i.e., measures related to perturbations acting on model parameters) to measure resilience in ecological systems (Dobrinevski & Frey, 2012;Constable & McKane, 2015;Cenci & Saavedra, 2018).
Here, we extend the connection between dynamical and structural indicators beyond the vicinity of an equilibrium state to study resilience in multispecies ecological systems. Specifically, we link recovery with dynamical stability and define it as the long-term return rate of a system to Accepted Article a quantitative reference state (i.e., abundance distribution at equilibrium) after small random perturbations on species abundances. We then separate recovery into full and partial recovery as to whether species abundances return fully or partially to such quantitative reference state, respectively. Next, we link resistance with structural stability and define it as the ability of a system to remain in a qualitative reference state (i.e., species composition at equilibrium) after random parameter perturbations of any magnitude. We then separate resistance into full and partial resistance as to whether all the species or at most one species remains in such qualitative reference state, respectively. Therefore, we define full resilience as the capacity of a system to maintain its full species composition through the recovery and resistance of all species. Then, we define partial resilience as the capacity of a system to maintain a partial species composition through the recovery and resistance of a subset of species (see Glossary).
We first illustrate our framework using theoretically parameterized ecological systems spanning multiple interaction types and number of species. In particular, we explore competition, mutualistic, and antagonistic systems with 3, 4, and 5 species. Overall, we find that when considering abundance and parameter perturbations together, recovery and resistance are interconnected measures of resilience. However, we show that these dynamical and structural indicators can provide complementary information about a system when separated into full and partial resilience. Then, we apply our framework to 17 experimentally parameterized microbial systems containing three interacting species each. We corroborate our theoretical results that full and partial resilience are complementary components of experimental systems. Finally, we discuss the implications of our findings and future avenues of research on resilience.
To study and measure resilience in ecological systems, it is necessary to define (i) a population dynamics model, (ii) a reference state, and (iii) the system's response as a function of the type and magnitude of perturbations (Justus, 2013). Focusing on the model, we assume that the dynamics of ecological systems are governed by any model topologically equivalent to the classic Lotka-Volterra (LV) dynamics. It has been shown that if the unstable and stable equilibria of the classic LV model can be mapped into the unstable and stable equilibria of a modified model, then this modified model is topologically equivalent to the classic LV (Cenci & Saavedra, 2018;
where N i is the abundance (or biomass) of species i, S is the number of species in the system, 132 and a ij is an element of the interaction matrix A = {a ij } representing the per capita effect of 133 species j on species i.
The phenomenological parameter r i is the intrinsic growth rate of species i, 134 representing how abiotic factors affect the balance between mortality and resource intake (Case, 135 2000; Coulson et al., 2017). 136 To investigate resilience in the classic LV model, we generate fully connected random interaction 137 matrices A by setting the diagonal elements to a ii = -1 ∀ i (introducing the biological principle of 138 self-regulation), whereas the off-diagonal elements a ij (i ̸ = j) are randomly drawn from a normal 139 distribution with mean zero and variance σ 2 (i.e., a ij ∼ N (0, σ 2 )). Note that the value of σ 2 sets 140 the relative strength of interspecific interactions (a ij ) given that intraspecific interactions are fixed 141 (i.e., a ii = -1). To explore different types of ecological systems, we introduce sign constraints: 142 (i) a ij < 0 for competition systems, (ii) a ij > 0 for mutualistic systems, and (iii) a ij a ji < 0 for 143 antagonistic systems (Murdoch et al., 2003; Allesina & Tang, 2012) (see Resilience in theoretical 144 systems). Applying these constraints is equivalent to sampling a ij from a half-normal distribution 145 and taking the appropriate sign reversal (Allesina & Tang, 2012; Song et al., 2020). For simplicity, 146 we only explore one type of antagonistic network structure given by a matrix A where a ij < 0 if 147 i > j and a ij > 0 if i < j (e.g., a trophic chain with omnivory). Note that, in these antagonistic 148 systems, the feasibility condition itself (i.e., N * i > 0 ∀ i; see Reference state) constraints the r i 149 values to be ecologically realistic as top predators and producers would necessarily have negative 150 and positive r i values, respectively. 151 To guarantee that reference equilibrium states are dynamically stable (see Reference state), we 152 follow a probabilistic criteria for random matrices and set the variance of the distribution of 153
interaction strengths proportional to system size (σ 2 = 1 S 2 ) (May, 1972). Note that having σ 2 154 to scale with the number of species S is an ecologically realistic assumption (Dougoud et al., 155 2018). We also tested a scenario in which interspecific interactions are stronger (i.e., σ 2 = 1 S ) and 156 obtained similar results (Figures S3 and S4). It is important to note that, under our framework, 157 we cannot increase the relative strength of interspecific interactions beyond a certain limit in 158 order to guarantee dynamical stability. Thus, we assume that every species in our theoretical 159 systems (see Resilience in theoretical systems) is self-regulated, including predators in antagonistic 160 systems (Barabás et al., 2017;Song & Saavedra, 2020).
161
We consider a feasible and dynamically stable equilibrium as our reference state (Song & Saavedra, 2018;Song et al., 2020). Feasibility corresponds to the capacity of a system to sustain all its constituent species over the long-run. Formally, feasibility implies the existence of a positive solution (i.e., N * i > 0 ∀ i under the equilibrium dN i /dt = 0 in Equation ( 1)). Note that feasibility is a necessary condition for persistence, permanence, and the existence of bounded orbits (Hofbauer & Sigmund, 1998). In turn, we define dynamical stability as the capacity of a system to return to its feasible equilibrium state after small random perturbations on species abundances. This is fulfilled when all the eigenvalues of the Jacobian matrix (J) when evaluated at the equilibrium
⊤ , as defined by dN * /dt = 0) have negative real parts (May, 1972;Case, 2000).
For the classic LV dynamics (Equation ( 1)), the Jacobian matrix at a feasible equilibrium is defined as:
where diag(N * ) is a diagonal matrix with N * 1 , . . . , N * S in its diagonal. The equilibrium of the classic LV model is given by the vector of species abundances N * = -A -1 r. This definition of reference state has the advantage of allowing us to represent it as a quantitative or qualitative reference state. Quantitatively, the focus is on the exact values of the feasible and dynamically stable equilibrium N * (i.e., the species abundance distribution at equilibrium). Qualitatively, the focus is on the existence of such a feasible equilibrium (N * i > 0 ∀ i) (i.e., the species composition at equilibrium).
Regarding the response of a system to abundance perturbations, we focus on the standard definition of recovery linked to dynamical stability (Strogatz, 2001). Specifically, we define recovery as the return rate of a system to a feasible and dynamically stable reference state after a small perturbation on abundances. Formally, we use as indicators the real part of the largest (λ 1 ) and second smallest (λ S-1 ) eigenvalues of the Jacobian matrix (J) evaluated at equilibrium (Equation ( 2)). The largest eigenvalue (if negative) measures the return rate along the slowest direction of recovery. Because species abundances will have recovered completely after going through the slowest direction, we use λ 1 as an indicator of full recovery. The second smallest eigenvalue (if negative) represents the return rate along the second fastest direction. Because species abundances will have only partially recovered after going through the second fastest direction, we
use λ S-1 as an indicator of partial recovery. Thus, for each possible feasible and dynamically stable state (N * ) of a system, we use λ 1 and λ S-1 as a measure of its full and partial recovery, respectively (Figures 1a and 1b). To be able to compare λ 1 and λ S-1 across different equilibrium states, we normalize N * to unit norm (i.e., N * ||N * || ) before computing these indicators.
Shifting our focus to the response of a system to random parameter perturbations of any magnitude, we base our analysis on the concept of structural stability and link it with resistance.
We define resistance as the smallest random parameter perturbation that a system can tolerate without affecting its qualitative reference state (defined as a feasible and dynamically stable equilibrium). For a given interaction matrix A, feasibility in the classic LV model will be satisfied as long as the r-vector falls inside the feasibility domain defined as (Song et al., 2018;Medeiros et al., 2020):
where -v i is the ith column vector of A and N * i is the feasible (i.e., positive) abundance of species i at equilibrium. Geometrically, D F (A) is a cone in S dimensions and any r-vector inside the cone gives rise to a feasible solution. Note that only the direction of the r-vector matters because if Equation ( 3) is satisfied for a given vector r, it is also satisfied for a scaled vector cr ∀ c > 0.
Therefore, the region inside the borders of D F (A) corresponds to the specific directions of r for which the system is feasible (Saavedra et al., 2017). Ecologically, this domain defines the range of conditions linked to abiotic factors, which are phenomenologically represented by the direction of the r-vector, compatible with the persistence of all species in the system characterized by A (Song et al., 2020;Medeiros et al., 2020). Note that we do not impose any restriction on the sign of the elements r i of the r-vector, they can be positive or negative depending on D F (A) (Song et al., 2018).
Importantly, starting from a feasible equilibrium state N * , a border of
with N * j = 0 and N * i > 0 ∀ i ̸ = j) represents a limit in the direction of r where at least one of the S species (i.e., species j) goes extinct (Rohr et al., 2016;Grilli et al., 2017;Tabi et al., 2020). Similarly, a vertex of D F (A) (i.e., r = N * i v i , with N * i > 0 and N * j = 0 ∀ j ̸ = i) represents a limit in the direction of r where a single species (i.e., species i) survives at most, depending on whether this species is self-sustained (Rohr et al., 2016;Grilli et al., 2017;Tabi et al., 2020).
For example, with three species, if one sets N * 1 = 0, the r-vector will lie on one of the borders
of the cone (D F (A)), which is given by r
N * 1 = 0 and N * 2 = 0, then the r-vector will lie on one of the vertices of D F (A), which is given by r = N * 3 v 3 , with N * 3 > 0. Thus, we focus on the shortest distances that a system can withstand before hitting a border or vertex under random perturbations to r. Note that we focus on random perturbations as we typically have no a priori information about the direction of environmental perturbations. Importantly, measuring such distances in the parameter space allows us to consider how resistant a system is to parameter perturbations of any magnitude.
Because it is only necessary to know the direction (not the magnitude) of r to know if a system
A is feasible, we normalize the intrinsic growth rates to unit norm (i.e., r ||r|| ). Then, the distance between an intrinsic growth rate vector inside the feasibility domain (r(N * )) and a border of vertices. We focus on the distance to the closest vertex (min{d v }). Therefore, the distance of a system to the closest vertex (min{d v }) represents the largest random parameter perturbation that a system can withstand before it is reduced to at most a single species-which we call partial resistance (Figures 1c and 1d).
First, we investigate the potential associations between recovery and resistance according to our measures described above across several theoretical systems (i.e., matrices A). To do so, we first generate three types of random matrices: competition systems, mutualistic systems, and antagonistic systems (see Population dynamics model). For each of these three types of system, we generate 100 random matrices A by sampling the interspecific interaction coefficients a ij from N (0, 1 S 2 ) for three different system sizes: S = 3, 4, and 5 species. Then, for each random matrix,
we sample 100 × 2 (S-2) feasible equilibria (i.e., N * > 0) uniformly on the unit S-dimensional hypersphere (i.e., ||N * || = 1) and solve for the intrinsic growth rates: r = -AN * . To guarantee dynamical stability in addition to feasibility, we only sample equilibrium states for which the real part of the largest eigenvalue of J = diag(N * )A is negative. We eliminate random matrices for which all sampled equilibrium states had a non-negative largest eigenvalue. Note that we increase the number of sampled equilibrium states exponentially with S to account for the fact that sparsity among points grows exponentially with the number of dimensions.
For each feasible and dynamically stable equilibrium state that we sample, we calculate its full recovery (largest eigenvalue, λ 1 ), partial recovery (second smallest eigenvalue, λ S-1 ), full resistance (distance to closest border, min{d b }), and partial resistance (distance to closest vertex, min{d v }).
To investigate whether recovery and resistance are linked, we compute the Pearson correlation coefficient (hereafter correlation) between full recovery and full resistance (ρ(λ 1 , min{d b })), as well as the correlation between partial recovery and partial resistance (ρ(λ S-1 , min{d v })). We compute such correlations separately for each random system A. A strong negative correlation indicates that the faster the recovery from abundance perturbations, the higher the resistance to parameter perturbations (Figure 1). Finally, we study the complementarity between our indicators by computing the partial correlation between λ 1 and λ S-1 controlling for the rank of λ 1 (ρ(λ 1 , λ S-1 | rank of λ 1 )). This control is necessary to account for the fact that λ 1 ≥ λ S-1 for each equilibrium state by definition. We also compute the partial correlation between min{d b } and min{d v } controlling for the rank of min{d v } (ρ(min{d b }, min{d v } | rank of min{d v })). Note that min{d v } ≥ min{d b } for each equilibrium state. A value close to zero of these partial correlations would imply the that full and partial components are complementary.
Next, we investigate the complementarity and potential application of our indicators of resilience using experimental ecological systems. For this purpose, we use 17 feasible and dynamically stable systems where each system consists of three interacting soil-dwelling bacteria species (Friedman et al., 2017). These publicly available data come from a very detailed and controlled study performing persistence experiments by co-inoculating different combinations of heterotrophic soildwelling bacteria species at varying initial fractions and propagating them through five growthdilution cycles (Friedman et al., 2017). Each 3-species experimental system is characterized by an interaction matrix A and an intrinsic growth rate vector r. The experimental values of r were inferred by fitting via least-squares the classic LV model (Equation (1)) to the observed abundance time series of species monocultures (Friedman et al., 2017) (Table S1). Each interaction matrix
A was inferred through pairwise tournaments by fitting via least-squares Equation (1) to the observed time series of species abundances (Friedman et al., 2017) (Table S2). Although more than 17 systems with three species are available from this data set, we only use the systems for which the equilibrium state obtained using the experimentally inferred matrix A and intrinsic growth rate vector r (N * = -A -1 r) is feasible and dynamically stable (Table S3).
Importantly, interaction matrices contained competition interactions (i.e., a ij , a ji < 0), antagonistic interactions (i.e., a ij < 0, a ji > 0), or a mix of both. It is also important to note that the strength of inferred interactions (a ij ) varies greatly across and within experimental systems (Table S2), allowing us to test our framework for different scenarios of interaction strengths. For each of the 17 systems (each combination of A and r), we calculate its full and partial resilience components following the same methodology as described above (see Resilience in theoretical systems). Additionally, for each experimental system, we randomly sampled 2,000 feasible and dynamically stable equilibria uniformly on the positive orthant of the unit sphere (i.e., N * > 0,
||N * || = 1, and λ 1 < 0) and measured our resilience indicators for each one of these sampled equilibria. The rationale behind sampling these 2,000 random equilibria was to compare the full/partial recovery and full/partial resistance observed for the experimental systems with other potential values of these indicators that could have been observed under different conditions.
Following our analysis with theoretical systems, for each of the 17 experimental systems, we study the complementarity of our indicators by computing the partial correlation between full recovery and partial recovery (λ 1 , λ 2 ) controlling for the rank of λ 1 (ρ(λ 1 , λ 2 | rank of λ 1 )) as well as the partial correlation between full resistance and partial resistance (min{d b }, min{d v }) controlling for the rank of min{d v } (ρ(min{d b }, min{d v } | rank of min{d v })). Then, in order to properly summarize the relationship between full and partial resilience across systems, we compute a relative value for each indicator ( λ1 , λ2 , min{ db }, and min{ dv }) by dividing each indicator by its minimum (recovery indicators) or maximum (resistance indicators) possible value that can be attained inside its corresponding feasibility domain. For example, the relative value of λ 1 ( λ1 ) for a given experimental equilibrium state N * is obtained by dividing λ 1 by the minimum value of λ 1 attained over all 2,000 randomly sampled equilibria N * , i.e., λ1 = λ 1 min{λ 1 } . Thus, relative values close to zero or one represent systems that have the minimum or maximum possible resilience, respectively. It is worth highlighting that this normalization does not allow us (and it is not intended) to compare the level of resilience between two different systems, only their resilience relative to the maximum possible value within the feasibility domain of the given system. Thus, a large variation in the experimental values of relative full and partial resilience can indicate a strong heterogeneity of resilience patterns across systems.
We found that fast recovery from abundance perturbations (i.e., state variables, N * i ) is associated with a high resistance to parameter perturbations (i.e., perturbations on intrinsic growth rates, r i ). Specifically, we found a strong negative correlation between full recovery (largest eigenvalue, λ 1 ) and full resistance (distance to closest border, min{d b }) across the three types of theoretical ecological systems (competition, mutualistic, and antagonistic systems; Figures 2a, S1, and S2).
Moreover, we found that this strong association between full recovery and full resistance holds for multiple random systems with different number of species (Figure 2b, Table S4). We obtained similar results for stronger interspecific interactions (i.e., σ 2 = 1 S ), but correlations can be weaker due to a smaller variation in λ 1 and highly asymmetric feasibility domains (i.e., high variation among the border lengths) (Figure S3). We can understand this result by noting that the Jacobian matrix associated with an r-vector on a border of the feasibility domain will have a row of zeros, which implies that the largest eigenvalue of this matrix will be zero (Figure 1a). In particular, all elements of the ith row of the Jacobian will be zero when the system is located at the border where N * i = 0. Recall that all eigenvalues of the Jacobian matrix need to be negative in order to guarantee dynamical stability and more negative values imply faster recovery.
Similarly, we found a strong negative correlation between partial recovery (second smallest eigenvalue, λ S-1 ) and partial resistance (distance to closest vertex, min{d v }) across the three types of theoretical systems (Figure 3a, S1, and S2). In addition, we confirmed that this strong relationship between partial recovery and partial resistance holds for multiple random systems with different number of species (Figure 3b, Table S4), even when interspecific interactions are stronger (Figure S4). Similarly to the association between λ 1 and min{d b }, the association between λ S-1 and min{d v } can be understood by noting that the Jacobian matrix associated with an r-vector on a vertex will have S -1 rows of zeros. This implies that the largest S -1 eigenvalues of the Jacobian will be zero at a vertex and, therefore, the second smallest eigenvalue will be associated with the distance to the closest vertex (Figures 1a and 1b).
We can additionally separate the relationship between recovery (full or partial) and resistance (full or partial) by the different vertices of the feasibility domain, noting that each vertex corresponds to the dominant species in the system (represented by different symbols in Figures 2a and 3a) (Rohr et al., 2016;Tabi et al., 2020). The extent to which the relationship between recovery and resistance may differ across vertices depends on the level of asymmetry in the feasibility domain (Rohr et al., 2016;Tabi et al., 2020). However, regardless of the asymmetry in the feasibility domain or the identity of the most abundant species, the relationship between recovery and result, we computed the partial correlation between recovery and resistance while controlling for 357 the identity of the most abundant species and, as expected, obtained slightly stronger correlation 358 values for λ S-1 and min{d v } but not for λ 1 and min{d b } (Figures S5 and S6; Table S5). In S5).
Moving to the experimental data, we confirmed the negative correlation between full recovery 363 (λ 1 ) and full resistance (min{d b }) (Figure 4a) as well as between partial recovery (λ 2 ) and partial 364 resistance (min{d v }) (Figure 4b) using 3-species microbial systems. In particular, we found 365 that ρ(λ 1 , min{d b }) and ρ(λ 2 , min{d v }) are strong and negative for all 17 experimental systems 366 (Figure 4c; Tables S6 and S7). Thus, the interconnections between recovery and resistance 367 remained strong for the experimental systems despite large differences in the type of interactions 368 (i.e., competition, antagonistic, or both interaction types) and interaction strengths (Figures S9
369 and S10). Furthermore, this complementarity was also observed in the diversity of combinations between full 375 and partial resilience across systems. That is, while some systems appear to exhibit high relative 376 full resilience and high relative partial resilience (i.e., relative values close to 1), others appear to 377 exhibit low relative values (i.e., relative values close to 0; Figures 5c and 5d). In addition, some 378 systems appear to exhibit an asymmetry between the indicators-i.e., high relative full resilience 379 and low relative partial resilience and vice versa (Figures 5c and 5d). Note that correlations 380 between full recovery and full resistance (ρ(λ 1 , min{d b })) as well as between partial recovery and 381 partial resistance (ρ(λ 2 , min{d v })) are strong and negative, but not perfectly correlated (Figure 382 4c). Therefore, we should expect a strong qualitative mapping, but not a perfect quantitative 383 match between the position of a system in Figure 5c and its position in Figure 5d. This is This article is protected by copyright. All rights reserved.
The resilience of ecological systems to perturbations is one of the most important yet broadly defined concepts in ecology and sustainability science (Hodgson et al., 2015;Folke et al., 2016;Pimm et al., 2019). Recently, the definition of resilience in ecology has been converging to the capacity of an ecological system to resist and recover from external perturbations (Hodgson et al., 2015;Capdevila et al., 2020). However, there are still many open questions about how to measure the response of multispecies systems to different types of perturbation, and the extent to which these responses provide complementary information (Donohue et al., 2016;Domínguez-García et al., 2019;Kéfi et al., 2019). Understanding the different components of resilience and their interconnections with respect to different types of perturbation is paramount to implement risk assessment and conservation strategies in ecological systems (Folke et al., 2004(Folke et al., , 2016)).
Here, we have introduced a new perspective in order to expand the way we measure resilience and understand the connections among its components in three major ways. First, by acknowledging that external perturbations can be of any type (i.e., abundance or parameter perturbations), we have extended and integrated dynamical (i.e., focused on abundance perturbations) and structural (i.e., focused on parameter perturbations) indicators of resilience (Arnoldi & Haegeman, 2016;Saavedra et al., 2017;Cenci & Saavedra, 2018). Specifically, we have found a clear link between a dynamical stability indicator (long-term return rate of a system to a dynamically stable equilibrium after small abundance perturbations) and a structural stability indicator (the largest random parameter perturbation that a system can withstand before losing species). Thus, differently from previous studies that have analyzed the associations between several resilience indicators (Hillebrand et al., 2018;Domínguez-García et al., 2019), here we have focused only on recovery and resistance and suggest a fundamental mathematical link between them when the dynamics follows the LV model. This suggests that other potential mathematical links may exist between frequently used dynamical and structural indicators of resilience (Domínguez-García et al., 2019;Arnoldi et al., 2016).
Second, we have found that recovery and resistance are interconnected when focusing on either full or partial components. We have defined full resilience as the capacity of a system to maintain its full species composition through the recovery and resistance of all species. In turn, we have defined partial resilience as the capacity of a system to maintain a partial species composition through the recovery and resistance of a subset of species. Specifically, we have shown that, under the LV model, fast (full or partial) recovery from abundance perturbations implies a high (full or partial) resistance to parameter perturbations. Therefore, we hypothesize that other hidden
connections between resilience indicators may be found if analyzed through the lens of full and partial components (Kéfi et al., 2019). From a practical point of view, this connection between recovery and resistance means that we can monitor both of these aspects of resilience using a small number of indicators (e.g., the eigenvalues of the Jacobian matrix of a system).
Third, we have found that full and partial resilience (either recovery or resistance) can be treated as complementary components. Thus, our study proposes a novel way to understand orthogonal dimensions of ecological resilience (Hillebrand et al., 2018;Domínguez-García et al., 2019).
Interestingly, full resilience is related to individual risk, whereas partial resilience is related to systemic risk, two important concepts in ecological, financial, and other complex systems (Levin, 1998;Beale et al., 2011). In fact, because full and partial resilience are complementary, focusing on only one of these components may provide a misleading risk assessment of ecological systems (Beale et al., 2011). For example, our results show that a system that minimizes systemic risk (i.e., is located as far as possible from the vertices of the feasibility domain) may still have a high risk of losing individual species (i.e., is located close to a border of the feasibility domain;
Figures 1c and 1d). In particular, the existence of a wide diversity of relationships between full and partial resilience in the studied experimental microbial systems suggests that future work may investigate the extent to which these systems may be more exposed to individual or systemic risk in the face of perturbations.
Our theoretical results can in principle be validated by performing perturbation experiments with simple experimental systems (e.g., microbial systems). For instance, one can measure the recovery rate of a system after small abundance perturbations (e.g., a single removal/addition of individuals/biomass (Steiner et al., 2006)) and the resistance to extinctions after structural perturbations (e.g., a fixed change in temperature (Tabi et al., 2020) or resource availability (Hoek et al., 2016)) over a set of baseline environmental conditions. Then, one can measure the association between dynamical and structural indicators, that is, recovery rate and resistance to extinctions. Furthermore, both indicators can be divided into their full (e.g., abundances fully recover, no single-species extinctions) and partial (e.g., abundances partially recover, no systemic collapse) components. For a given system, our theory purports an association between full (partial) recovery and full (partial) resistance (Figures 4a and 4b), and no particular association between full and partial components (recovery or resistance; Figures 5a and 5b).
Finally, it is worth mentioning that our indicators of resilience are formally derived under the assumptions of feasible and dynamically stable equilibria with the classic LV model (Case, 2000).
Thus, these assumptions have to be fulfilled in order to apply our framework using empirical data.
Nevertheless, multiple methods have been developed to estimate the structure of nonequilibrium 455 systems (e.g., Jacobian matrix) without assuming a parameterized population dynamics model 456 (i.e., a nonparametric approach) (Ives et al., 2003;Deyle et al., 2016;Ushio et al., 2018). In partic-457 ular, previous work has shown that the divergence of a nonequilibrium vector field (characterized 458 by the trace of the Jacobian matrix) is associated with the extent to which the trajectory of the 459 system changes after parameter perturbations (Cenci & Saavedra, 2019;Cenci et al., 2020)-a 460 measure of resistance. Thus, these previous results indicate that it may be possible to merge 461 concepts from dynamical and structural stability also under a nonparametric approach. In this 462 sense, we believe that our study may also add to the unification of parametric and nonparametric 463 approaches (Song & Saavedra, 2020) and could be expanded to be used with different types of 464 empirical data.
Glossary Asymptotic dynamical stability: The capacity of a dynamical system to return to a quantitative equilibrium reference state after a given perturbation acting on its state variables (i.e., species abundances) (Strogatz, 2001).
A measure of resilience related to perturbations acting on species abundances.
The set of directions of parameter values (here intrinsic growth rates) compatible with positive solutions for all species (i.e., feasible system).
The return rate along the slowest direction following a perturbation acting on species abundances. We measure full recovery as the largest eigenvalue (λ 1 ) of the Jacobian matrix evaluated at an equilibrium state.
The capacity of a system to maintain its full species composition through the recovery and resistance of all species. Full resilience is partitioned into full recovery, which is related to abundance perturbations (i.e., asymptotic dynamical stability), and full resistance, which is related to parameter perturbations (i.e., structural stability).
The largest random parameter perturbation that the system can withstand before losing any species. We measure full resistance as the distance of an equilibrium state to the closest border (min{d b }) of the feasibility domain (D F (A)).
The return rate along the second fastest direction following a perturbation acting on species abundances. We measure partial recovery as the second smallest eigenvalue (λ S-1 ) of the Jacobian matrix evaluated at an equilibrium state.
The capacity of a system to maintain a partial species composition through the recovery and resistance of a subset of species. Partial resilience is partitioned into partial recovery, which is related to abundance perturbations (i.e., asymptotic dynamical stability), and partial resistance, which is related to parameter perturbations (i.e., structural stability).
The largest random parameter perturbation that a system can withstand before at most a single species survives. We measure partial resistance as the distance of an equilibrium state to the closest vertex (min{d v }) of the feasibility domain (D F (A)).
Perturbation: Any event that impacts the species abundances directly or the rules that govern population dynamics in an ecological system (i.e., the structure of the ecological system) (Justus, 2013).
A mathematical idealized description of the causal relationships (mechanistic or phenomenological) connecting the change of a population through time as a function of abiotic and biotic factors (Case, 2000).
The rate of return (or time, if inverted) of an ecological system to a reference state after a given perturbation (Hodgson et al., 2015).
Reference state: The state of an ecological system against which the perturbed system will be compared to (Justus, 2013). The reference state can be quantitative (e.g., the species abundances at equilibrium) or qualitative (e.g., the species composition of the system).
Resilience: Ability of an ecological system (i.e., set of interacting species) to resist and recover from external perturbations (Hodgson et al., 2015).
The capacity of an ecological system to resist changes in relation to its reference state after a given perturbation (Hodgson et al., 2015).
A measure of resilience related to perturbations acting on model parameters.
The capacity of a dynamical system to retain the topology of the phase portrait (i.e., the qualitative behavior of the trajectories) after a given perturbation to its structure (i.e., governing laws or model parameters) (Strogatz, 2001).
Competition system Mutualistic system Antagonistic system
This article is protected by copyright. All rights reserved.