Ecological resilience is defined as the ability of ecological systems to tolerate disturbances and still maintain the same relationships between state variables (Holling, 1973; Annu. Rev. Ecol. Evol. Syst., https://doi.org/10.1146/annurev.es.04.110173.000245). Measuring ecological resilience requires considering two important components: the dynamic trends of the system, including its cyclic behaviors, and the random forces representing positive and negative feedback relationships between the components of the system (Holling, 1973). As such, assessing ecological resilience must account for the system dynamic regimes.
An ecological dynamic regime (EDR) is defined as the “fluctuations of an ecological system around some trend or average resulting from the interaction between internal processes and external forces that, in the absence of perturbations, keep the system within the basin of attraction” (Sánchez-Pinillos et al., 2023). The main dynamic trends characterizing an EDR are useful to identify cyclic behaviors and other more complex dynamics, such as transient dynamics. Additionally, the fluctuations and variability of ecological dynamics resulting from the interaction of multiple factors define the shape, size, and characteristics of dynamic regimes and potential domains of attraction (Sánchez-Pinillos et al., 2024).
The EDR framework is a set of algorithms and metrics to characterize and compare ecological dynamic regimes from empirical data so they can be used as the reference to assess ecological resilience accounting for both the dynamic trends of the system and the feedback relationships between its components.
To know more about the relevance of dynamic regimes for assessing ecological resilience, see this publication:
Additional information about the EDR framework can be found in this publication:
This vignette aims to illustrate how ecological resilience can be
assessed using the EDR framework implemented in ecoregime
.
In particular, this vignette focuses on the quantitative indices
proposed in Sánchez-Pinillos et al. (2024) and their geometric
rationale, taking ecological dynamic regimes as the reference to
evaluate ecological resilience.
You can install ecoregime
directly from CRAN or from my
GitHub account (development version):
# install.packages("ecoregime")
# devtools::install_github(repo = "MSPinillos/ecoregime", dependencies = T, build_vignettes = T)
Once you have installed ecoregime
you will have to load
it:
citation("ecoregime")
#> To cite 'ecoregime' in publications use:
#>
#> Sánchez-Pinillos M, Kéfi S, De Cáceres M, Dakos V (2023). "Ecological
#> dynamic regimes: Identification, characterization, and comparison."
#> _Ecological Monographs_, e1589. <https://doi.org/10.1002/ecm.1589>.
#>
#> Sánchez-Pinillos M, Dakos V, Kéfi S (2024). "Ecological dynamic
#> regimes: A key concept for assessing ecological resilience."
#> _Biological Conservation_, 110409.
#> <https://doi.org/10.1016/j.biocon.2023.110409>.
#>
#> Sánchez-Pinillos M (2023). _ecoregime: Analysis of Ecological Dynamic
#> Regimes_. <https://doi.org/10.5281/zenodo.7584943>.
#>
#> To see these entries in BibTeX format, use 'print(<citation>,
#> bibtex=TRUE)', 'toBibtex(.)', or set
#> 'options(citation.bibtex.max=999)'.
To assess the ecological resilience of a system, we need to identify its ecological dynamic regime (EDR) so it can be used as the reference to compare the disturbed dynamics. That is, the EDR to which the disturbed trajectories belonged before being disturbed.
Given a stability landscape composed of alternative dynamic regimes,
we can compute the dynamic dispersion (dDis) of the pre-disturbance
portion of the disturbed trajectories in relation to each alternative
EDR. The function dDis()
allows us to quantify the degree
of membership of any trajectory within an EDR. You can see more
information about dDis()
in
vignette("EDR_framework")
.
In the data included in ecoregime
(i.e.,
EDR_data
), you will find three abundance matrices
corresponding to a set of community trajectories forming three
ecological dynamic regimes (EDR1, EDR2, EDR3). Additionally, you will
find an abundance matrix associated with three disturbed communities
(EDR3_disturbed
). Assuming that EDR1, EDR2, and EDR3 are
three alternative dynamic regimes in a stability landscape, we will see
how dDis()
can be used to identify the EDR with which the
disturbed trajectories are associated.
# Species abundances in undisturbed states of the disturbed trajectories
# (Undisturbed states are identified by disturbed_states = 0)
abun_undist <- EDR_data$EDR3_disturbed$abundance[disturbed_states == 0]
selcols <- names(EDR_data$EDR1$abundance)
## EDR1 ------------------------------------------------------------------------
# Species abundances in EDR1 and the undisturbed states of disturbed trajectories
abun1_undist <- rbind(EDR_data$EDR1$abundance, abun_undist[, ..selcols])
# State dissimilarities in EDR1 and the undisturbed states of disturbed trajectories
d1_undist <- vegan::vegdist(x = abun1_undist[, paste0("sp", 1:12)], method = "bray")
# dDis of the disturbed trajectories in relation to EDR1
dDis1 <- sapply(unique(abun_undist$traj), function(iundist){
dDis(d = d1_undist, d.type = "dStates",
trajectories = abun1_undist$traj, states = abun1_undist$state,
reference = as.character(iundist))
})
## EDR2 ------------------------------------------------------------------------
# Species abundances in EDR2 and the undisturbed states of disturbed trajectories
abun2_undist <- rbind(EDR_data$EDR2$abundance, abun_undist[, ..selcols])
# State dissimilarities in EDR2 and the undisturbed states of disturbed trajectories
d2_undist <- vegan::vegdist(x = abun2_undist[, paste0("sp", 1:12)], method = "bray")
# dDis of the disturbed trajectories in relation to EDR2
dDis2 <- sapply(unique(abun_undist$traj), function(iundist){
dDis(d = d2_undist, d.type = "dStates",
trajectories = abun2_undist$traj, states = abun2_undist$state,
reference = as.character(iundist))
})
## EDR3 ------------------------------------------------------------------------
# Species abundances in EDR3 and the undisturbed states of disturbed trajectories
abun3_undist <- rbind(EDR_data$EDR3$abundance, abun_undist[, ..selcols])
# State dissimilarities in EDR3 and the undisturbed states of disturbed trajectories
d3_undist <- vegan::vegdist(x = abun3_undist[, paste0("sp", 1:12)], method = "bray")
# dDis of the disturbed trajectories in relation to EDR3
dDis3 <- sapply(unique(abun_undist$traj), function(iundist){
dDis(d = d3_undist, d.type = "dStates",
trajectories = abun3_undist$traj, states = abun3_undist$state,
reference = as.character(iundist))
})
## Compare dynamic dispersion --------------------------------------------------
# Compare dDis values for the three EDRs
dDis_df <- data.frame(EDR1 = dDis1, EDR2 = dDis2, EDR3 = dDis3)
If we compare the values of dDis for each disturbed trajectory (ref. 31, ref. 32, ref. 33) and EDR (EDR1, EDR2, EDR3), we see that the lowest dDis values are associated with EDR3. We can take EDR3 as the reference to assess the ecological resilience of the three disturbed communities.
EDR1 | EDR2 | EDR3 | |
---|---|---|---|
dDis (ref. 31) | 0.649 | 0.748 | 0.356 |
dDis (ref. 32) | 0.734 | 0.556 | 0.268 |
dDis (ref. 33) | 0.711 | 0.436 | 0.412 |
Once we have identified the EDR of reference, we can compute the
representative trajectories using the function retra_edr()
.
See vignette("EDR_framework")
for more information.
# State dissimilarities for EDR3 (considering only the undisturbed trajectories)
d_EDR3 <- vegan::vegdist(EDR_data$EDR3$abundance[, paste0("sp", 1:12)])
# Representative trajectories
retra <- retra_edr(d = d_EDR3,
trajectories = EDR_data$EDR3$abundance$traj,
states = EDR_data$EDR3$abundance$state, minSegs = 5)
Although there are five representative trajectories, we will select “T4” as the reference for being the longest and covering all regions of the EDR relatively well.
# Summarize retra
summary(retra)
#> ID Size Length Avg_link Sum_link Avg_density Max_density Avg_depth
#> T1 T1 8 0.56912543 0.12435644 0.24871287 7.250000 8 4.750000
#> T2 T2 6 0.39900498 0.10718905 0.21437811 7.666667 8 4.666667
#> T3 T3 4 0.08970149 0.03482587 0.03482587 6.500000 7 6.500000
#> T4 T4 16 1.11654181 0.08961127 0.53766760 6.625000 8 5.125000
#> T5 T5 6 0.35468137 0.08723881 0.17447761 7.333333 10 5.333333
#> Max_depth
#> T1 6
#> T2 6
#> T3 7
#> T4 7
#> T5 7
# Define T4 as the unique representative trajectory and generate an object of class 'RETRA'
retra_ref <- define_retra(data = retra$T4$Segments, d = d_EDR3,
trajectories = EDR_data$EDR3$abundance$traj,
states = EDR_data$EDR3$abundance$state,
retra = retra)
# Plot EDR3 and its representative trajectories
plot(retra, d = d_EDR3,
trajectories = EDR_data$EDR3$abundance$traj,
states = EDR_data$EDR3$abundance$state, select_RT = "T4",
main = "Representative trajectories in EDR3")
legend("topleft", c("Representative trajectory 'T4'",
"Other representative trajectories",
"Individual trajectories in EDR3"),
lty = 1, col = c("red", "black", "grey"), bty = "n")
When an ecological system is disturbed, its dynamics can deviate from the reference EDR and its representative trajectories. Depending on the immediate impact of the disturbance and the changes in the state variables following the release of the disturbance, post-disturbance trajectories can be characterized through four complementary indices: resistance (Rt), amplitude (A), recovery (Rc), and net change (NC).
The resistance index (Rt) quantifies the impact of the disturbance on the system based on the immediate changes in the state variables. That is, it is a measure of how similar the disturbed (1) and the pre-disturbance (0) states are (Sánchez-Pinillos et al., 2019). As such, the resistance index does not depend on the position of the system within the EDR.
\[Rt = 1 - d_{pre, dist}\]
where \(d_{pre,dist}\) is the dissimilarity between the pre-disturbance (0) and the disturbed (1) states.
Given the three disturbed trajectories in EDR3 (i.e.,
EDR3_disturbed
), we can compute the resistance index using
the function resistance()
:
# To calculate resistance, we need a state dissimilarity matrix for the disturbed trajectories
d_disturbed <- vegan::vegdist(EDR_data$EDR3_disturbed$abundance[, paste0("sp", 1:12)],
method = "bray")
# Compute resistance
# Note that the disturbed states are identified by disturbed_states = 1
Rt <- resistance(d = d_disturbed,
trajectories = EDR_data$EDR3_disturbed$abundance$traj,
states = EDR_data$EDR3_disturbed$abundance$state,
disturbed_trajectories = unique(EDR_data$EDR3_disturbed$abundance$traj),
disturbed_states = EDR_data$EDR3_disturbed$abundance[disturbed_states == 1]$state)
The three hypothetical systems show a relatively high resistance (close to 1) to the disturbance:
Disturbed trajectories | Rt |
---|---|
31 | 0.958 |
32 | 0.812 |
33 | 0.693 |
The amplitude (A) quantifies how much the system is deviated from its trajectory during the disturbance assuming that, in the absence of disturbances, the system would keep a constant distance to a representative trajectory taken as the reference.
The amplitude can be calculated in absolute terms as the difference of the dissimilarity between the disturbed state (1) and the representative trajectory (\(d_{dist,RT}\)) and the dissimilarity between the pre-disturbance state (0) and the representative trajectory (\(d_{pre,RT}\)):
\[A_{abs} = d_{dist,RT} - d_{pre,RT}\]
Alternatively, the amplitude can be calculated in relation to the impact of the disturbance (\(d_{pre,dist}\)). In this case, the amplitude quantifies the ability of the system to remain close to the representative trajectory in relation to the changes in the state variables provoked by the disturbance.
\[A_{rel} = \frac{d_{dist,RT} - d_{pre,RT}}{d_{pre,dist}}\]
In any case, positive amplitude values indicate that the system is deviated towards the boundaries of the EDR, whereas negative values indicate the deviation towards the representative trajectory taken as the reference.
Amplitude can be calculated using the function
amplitude()
:
# We need a state dissimilarity matrix containing the states of the disturbed
# trajectories and the representative trajectory taken as the reference
abun <- rbind(EDR_data$EDR3$abundance, EDR_data$EDR3_disturbed$abundance, fill = T)
d <- vegan::vegdist(abun[, paste0("sp", 1:12)], method = "bray")
# Compute amplitude
A <- amplitude(d = d,
trajectories = abun$traj,
states = abun$state,
disturbed_trajectories = abun[disturbed_states == 1]$traj,
disturbed_states = abun[disturbed_states == 1]$state,
reference = retra_ref, method = "nearest_state")
The three considered systems show positive amplitude values, indicating that the disturbance lead them towards the boundaries of the EDR. Whereas the absolute amplitude of trajectory 31 (\(A_{abs}(31) = 0.027\)) is smaller than the absolute amplitude of trajectory 32 (\(A_{rel}(32) = 0.119\)), both have similar relative values (\(A_{rel}(31) = 0.653\); \(A_{rel}(32) = 0.631\)). This result indicates that despite the small deviation of trajectory 31 from the representative trajectory during the disturbance, such deviation is disproportionately high in relation to the impact of the disturbance on the state variables.
Disturbed trajectories | Reference | Aabs | Arel |
---|---|---|---|
31 | T4.1 | 0.027 | 0.653 |
32 | T4.1 | 0.119 | 0.631 |
33 | T4.1 | 0.268 | 0.871 |
The recovery index (Rc) quantifies the ability of the system to reorganize itself after a disturbance and evolve towards the representative trajectory representing its dominant dynamic trends.
The recovery can be calculated in absolute terms as the difference of the dissimilarity between the disturbed state (1) and the representative trajectory (\(d_{dist,RT}\)) and the dissimilarity between one of the post-disturbance states (> 1) and the representative trajectory (\(d_{post,RT}\)):
\[Rc_{abs} = d_{dist,RT} - d_{post,RT}\]
Alternatively, the recovery can be calculated in relation to the changes in the state variables that the system must perform to return towards the representative trajectory from the disturbed state (\(d_{dist, post}\)). In this case, the index penalizes the systems that require major restructuring of the state variables to return towards the representative trajectory (positive recovery) and gives less negative values to the systems minimizing the escape from the dynamic regime through major changes in its state variables.
\[Rc_{rel} = \frac{d_{dist,RT} - d_{post,RT}}{d_{dist, post}}\]
Both absolute and relative recovery indices are positive when the system evolves towards the representative trajectory. Otherwise, negative values indicate that the system evolves towards the boundaries of the EDR.
Recovery can be calculated using the function recovery()
considering all states after the disturbed state (> 1):
# Compute recovery using the same data used for amplitude
Rc <- recovery(d = d,
trajectories = abun$traj,
states = abun$state,
disturbed_trajectories = abun[disturbed_states == 1]$traj,
disturbed_states = abun[disturbed_states == 1]$state,
reference = retra_ref, method = "nearest_state")
We can plot the variation of recovery over time considering all post-disturbance states:
# Number of states after the disturbed state
Rc <- data.table::data.table(Rc)
Rc[, ID_post := 1:(.N), by = disturbed_trajectories]
# Plot absolute recovery over time
plot(x = range(Rc$ID_post), y = range(Rc$Rc_abs), type = "n",
xlab = "Nb. states after disturbance", ylab = "Absolute recovery",
main = "Variation of absolute recovery")
for (i in unique(Rc$disturbed_trajectories)) {
lines(Rc[disturbed_trajectories == i, c("ID_post", "Rc_abs")],
col = which(unique(Rc$disturbed_trajectories) %in% i) + 1)
}
legend("bottomleft", legend = unique(Rc$disturbed_trajectories), lty = 1,
col = seq_along(unique(Rc$disturbed_trajectories)) + 1, bty = "n")
# Plot relative recovery over time
plot(x = range(Rc$ID_post), y = range(Rc$Rc_rel), type = "n",
xlab = "Nb. states after disturbance", ylab = "Relative recovery",
main = "Variation of relative recovery")
for (i in unique(Rc$disturbed_trajectories)) {
lines(Rc[disturbed_trajectories == i, c("ID_post", "Rc_rel")],
col = which(unique(Rc$disturbed_trajectories) %in% i) + 1)
}
legend("topright", legend = unique(Rc$disturbed_trajectories), lty = 1,
col = seq_along(unique(Rc$disturbed_trajectories)) + 1, bty = "n")
Net change (NC) quantifies how much the system is deviated from its trajectory after the release of the disturbance assuming that, in the absence of disturbances, the system would keep a constant distance to a representative trajectory taken as the reference.
Like amplitude and recovery, net change can be calculated in absolute or relative terms. The absolute net change is expressed as the difference of the dissimilarity between one of the post-disturbance states (> 1) and the representative trajectory (\(d_{post,RT}\)) and the dissimilarity between the pre-disturbance state (0) and the representative trajectory (\(d_{pre,RT}\)):
\[NC_{abs} = d_{post,RT} - d_{pre,RT}\]
The relative net change is calculated in relation to the changes produced in the state variables between the pre-disturbance and the post-disturbance states (\(d_{pre,post}\)). In this way, the index penalizes the systems that deviated from the expected trajectory despite being very similar to the pre-disturbance state.
\[NC_{rel} = \frac{d_{post,RT} - d_{pre,RT}}{d_{pre,post}}\]
As in the amplitude index, positive values indicate that the system is deviated towards the boundaries of the EDR and negative values indicate the deviation towards the representative trajectory taken as the reference.
The function net_change()
included in
ecoregime
can be used to calculate the net change of all
post-disturbance states in relation to the pre-disturbance state:
# Compute net change using the same data used for amplitude
NC <- net_change(d = d,
trajectories = abun$traj,
states = abun$state,
disturbed_trajectories = abun[disturbed_states == 1]$traj,
disturbed_states = abun[disturbed_states == 1]$state,
reference = retra_ref, method = "nearest_state")
We can plot the variation of net change over time considering all post-disturbance states:
# ID post-disturbance states
NC <- data.table::data.table(NC)
NC[, ID_post := 1:(.N), by = disturbed_trajectories]
# Plot absolute net change over time
plot(x = range(NC$ID_post), y = range(NC$NC_abs), type = "n",
xlab = "Nb. states after disturbance", ylab = "Absolute net change",
main = "Variation of absolute net change")
for (i in unique(NC$disturbed_trajectories)) {
lines(NC[disturbed_trajectories == i, c("ID_post", "NC_abs")],
col = which(unique(NC$disturbed_trajectories) %in% i) + 1)
}
legend("topleft", legend = unique(NC$disturbed_trajectories), lty = 1,
col = seq_along(unique(NC$disturbed_trajectories)) + 1, bty = "n")
# Plot relative net change over time
plot(x = range(NC$ID_post), y = range(NC$NC_rel), type = "n",
xlab = "Nb. states after disturbance", ylab = "Relative net change",
main = "Variation of relative net change")
for (i in unique(NC$disturbed_trajectories)) {
lines(NC[disturbed_trajectories == i, c("ID_post", "NC_rel")],
col = which(unique(NC$disturbed_trajectories) %in% i) + 1)
}
legend("bottomleft", legend = unique(NC$disturbed_trajectories), lty = 1,
col = seq_along(unique(NC$disturbed_trajectories)) + 1, bty = "n")
Altogether, the metrics of amplitude, recovery, and net change define the geometry of disturbed trajectories in relation to the reference EDR. Additionally, the resistance index helps interpret the response of the system considering the immediate impact of the disturbance.
As an example, we will take the three disturbed trajectories used in the previous section. For simplicity, we will only focus on resistance and the absolute values of amplitude, recovery, and net change and we will consider the 14th post-disturbance state (last post-disturbance state of Trajectory 33) of all trajectories to calculate recovery and net change.
# Merge the results for resistance, amplitude, recovery, and net change
results <- Reduce(function(x, y) merge(x, y, all = T),
list(Rt, A, Rc, NC))
results <- results[which(results$ID_post == 14),
c("disturbed_trajectories", "Rt", "A_abs", "Rc_abs", "NC_abs")]
Disturbed trajectories | Rt | Aabs | Rcabs | NCabs |
---|---|---|---|---|
31 | 0.958 | 0.027 | -0.474 | 0.502 |
32 | 0.812 | 0.119 | 0.099 | 0.020 |
33 | 0.693 | 0.268 | 0.283 | -0.016 |
Based on the values of the indices, we can figure out the shape of the disturbed trajectories. For example, Trajectory 31 was very resistant to the immediate impact of the disturbance (in terms of changes in the state variables) (\(Rt = 0.958\)) and the amplitude provoked by the disturbance was very small (\(A = 0.027\)). However, unlike Trajectory 32 and Trajectory 33, the recovery was negative (\(Rc = -0.474\)) and the net change positive and relatively high (\(NC = 0.502\)). Thus, despite the relatively low impact of the disturbance, the values of recovery and net change indicate that Trajectory 31 goes away from the reference EDR and the system potentially changes to an alternative dynamic regime.
Both Trajectory 32 and Trajectory 33 showed positive values of recovery and net change values close to zero. In both cases, amplitude and recovery have similar values, respectively (Trajectory 32: \(A = 0.119\), \(Rc = 0.099\); Trajectory 33: \(A = 0.268\), \(Rc = 0.283\)), indicating that these systems are able to reorganize and remain within the EDR.
Trajectory 32 was relatively resistant to the impact of the disturbance (\(Rt = 0.812\)) and showed a relatively low deviation from the representative trajectory (\(A = 0.119\)). In contrast, Trajectory 33 was more severely impacted by the disturbance (\(Rt = 0.693\); \(A = 0.268\)) but showed a high recovery capacity (\(Rc = 0.283\)). Thus, while both communities could be considered resilient, they showed different strategies. Whereas Trajectory 32 represents a resistant dynamic within the EDR, Trajectory 33 represents a more variable system whose dynamics visit the borders of the EDR to eventually return towards the representative trajectory.
To support our evaluation of the ecological resilience of disturbed trajectories, we can calculate the dynamic dispersion (dDis) of the post-disturbance portion of the disturbed trajectories in relation to the trajectories in the EDR.
# Species abundances in post-disturbance states of the disturbed trajectories
# The states after the release of the disturbance are identified by disturbed_states > 1
abun_post <- EDR_data$EDR3_disturbed$abundance[disturbed_states > 1]
selcols <- names(EDR_data$EDR3$abundance)
# Species abundances in EDR3 and the post-disturbance states of disturbed trajectories
abun3_post <- rbind(EDR_data$EDR3$abundance, abun_post[, ..selcols])
# State dissimilarities in EDR3 and the post-disturbance states of disturbed trajectories
d3_post <- as.matrix(vegan::vegdist(x = abun3_post[, paste0("sp", 1:12)], method = "bray"))
# dDis of each disturbed trajectory
dDis_dist <- sapply(unique(abun_post$traj), function(idist){
rm_disturbed <- unique(abun_post$traj[which(abun_post$traj != idist)])
irm <- which(abun3_post$traj %in% rm_disturbed)
dDis(d = d3_post[-irm, -irm], d.type = "dStates",
trajectories = abun3_post$traj[-irm],
states = abun3_post$state[-irm],
reference = as.character(idist))
})
The values of dDis confirm the resilience of Trajectory 32 and show that, despite being able to recover after the disturbance, the dynamic represented by Trajectory 33 is anomalous in relation to the other dynamics in the EDR.
Disturbed trajectories | dDis |
---|---|
31 | 0.463 |
32 | 0.257 |
33 | 0.474 |
Finally, we can illustrate the responses of the three disturbed communities by representing their trajectories in relation to the reference EDR and its representative trajectory:
# Plot EDR3 and its representative trajectories
plot(retra_ref, d = d,
trajectories = abun$traj,
states = abun$state,
traj.colors = c(rep("grey", 30), 2:4),
main = "Comparison of disturbed trajectories and EDR3")
legend("topleft", c("Representative trajectory",
"Individual trajectories in EDR3",
"Trajectory 31",
"Trajectory 32",
"Trajectory 33"),
lty = 1, col = c("black", "grey", 2:4), bty = "n")
This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No 891477 (RESET project).