Ecological resilience to pulse disturbances

1.1. Ecological resilience and the EDR framework

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:

• Sánchez-Pinillos M., Dakos, V., Kéfi, S. 2024. Ecological dynamic regimes: A key concept for assessing ecological resilience. Biological Conservation. https://doi.org/10.1016/j.biocon.2023.110409

Additional information about the EDR framework can be found in this publication:

• Sánchez-Pinillos M., Kéfi, S., De Cáceres, M., Dakos, V. 2023. Ecological Dynamic Regimes: Identification, characterization, and comparison. Ecological Monographs. https://doi.org/10.1002/ecm.1589

1.2. About this vignette

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:

library(ecoregime)

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)'.

2.1. Identifying the reference ecological dynamic regime

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.

2.2. Define the main dynamic trends through representative trajectories

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")

3.1. Resistance

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:

3.2. Amplitude

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.

3.3. Recovery

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")

3.4. Net change

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")