Introduction

The original paper by Scheffer et al (a review in Nature in 2001) describes how smooth changes in environmental conditions can create drastic / catastrophic changes in ecosystem state. These changes can display hysteresis, whereby responses to changing environmental conditions strongly depend on historical events. The paper reviews evidence of such shifts in lakes, coral reefs, woodlands, deserts, and oceans. It goes on to discuss emerging patterns and implications for management.

The initial purpose of this reproduction is to use the model in Box 1, described as a minimal model of an ecosystem showing hysteresis, to create versions of figures 1, 2, and 3. Some time series of dynamics might also be nice.

There is much more that could be added:

A minimal model of an ecosystem showing hysteresis

\(\frac{dx}{dt} = a - bx + rf(x)\)

\(x\) is an “unwanted” ecosystem property. \(a\) is an environmental factor that “promotes” \(x\). \(b\) is the rate of decay of \(x\). \(r\) is the rate at which \(x\) recovers, as a function of \(x\).

“For a lake, one can think of \(x\) as nutrients suspended in phytoplankton causing turbidity, of \(a\) as nutrient loading, and \(b\) as nutrient removal rate, and of \(r\) as internal nutrient recycling. For desertification, one could interpret \(x\) as barren soil, \(a\) as vegetation destruction, \(b\) as recolonization of barren soil by plants and \(r\) as erosion by wind and runoff.”

Hysteresis and alternative stable states can occur if \(f(x)\) is a function with a threshold, e.g. the hill function:

\(f(x) = \frac{x^p}{x^p + h^p}\)

Explore f(x), the function that determines the internal feedback

The non-linearity of this function is all important for the system dynamics. If the function is smooth (little non-linearity when \(p\) is small) then the rate \(r * f(x)\) varies smoothly with \(x\) and dynamics are quite linear with the environmental condition \(a\). If the function is very nonlinear (large \(p\)), then the system responds quite nonlinearly to changes in the environmental condition. Here is the function for some different values of \(p\).

Figure 1: environment - equilibrium state relationships

Below are some relationships between environmental conditions (\(a\)) and equilibrium ecosystem state (\(x\)), going from smooth relationship (low values of \(p\)), to nonlinear relationship (medium values of \(p\)), to folded relationship (high values of \(p\)). Solid lines indicate a stable equilibrium, dashed lines an unstable one.

## Warning: `cols` is now required.
## Please use `cols = c(roots)`

Figure 2: dynamically shifting between states

Here we simulate the model through time, with increasing and then decreasing environmental condition \(a\) (green line in the first graph) driving the dynamics. A relatively large value of \(p\) is used, so we get a folded relationship between environmental conditions (\(a\)) and ecosystem state (\(x\)) (light blue line in second graph). The system dynamics are the thin black line; there are two ways of shifting between the alternate states. One way as the environmental condition increases from low to high, and the other as the environmental conditions decreases from high to low.

Figure 3: basins of attraction

A locally stable equilibrium is sometimes viewed as a depression, or basin, in a surface, and a ball on the surface is attracted to and sits at the bottom of the basin. A locally unstable equilibrium is a peak, where a ball can sit, though will fall off given the slightest push. Below we can see the surface, basins of attraction, and the peaks, for some different values of \(a\). Imagine the position of a ball when \(a\) is low: it will sit in the single basin around \(x=0.2\). An increase in \(a\) moves the bottom of the basin slowly to the right (i.e. higher equilibrium value of \(x\)). At intermediate values of \(a\) (e.g. 0.45, 0.5) there are two basins with a peak in between. Which basin the ball is in depends on where it was before. At high values of \(a\) there is again only one basin.

Effects of \(p\) and rate of change of \(a\)

Show relationship between feedback strength and a measure of the amount of nonlinearity in the state-environmental relationship, and the similarity in state dynamics between forward and reverse direction of environmental change.

## get non-linearity and hysteresis for a value of p
get_nl_hyst_by_p <- function(p, inv_rate=20000, a_forcing) {
  
  #print(p)
  
  this.p <- p
  parameters <- c(a = a, b = b, r = r, p = this.p, h = h)
  state <- c(x = 0.1)
  
  model <- function(t, state, parameters) {
    with(as.list(c(state, parameters)), {
      dx <- dx_dt(x, a_forcing2(t), b, r, p, h)
      list(c(dx), a=a_forcing2(t))
    })
  }
  
  times <- seq(0, inv_rate, length = 2001)
  
  a_forcing1 <- matrix(ncol=2, byrow=T, data=c(0, a_forcing[1],
                                               mean(times), a_forcing[2],
                                               max(times), a_forcing[3]))
  a_forcing2 <- approxfun(x = a_forcing1[,1], y = a_forcing1[,2], method = "linear", rule = 2)
  out <- as.data.frame(ode(y = state, times = times, func = model, parms = parameters))
  
  ## nonlinearity
  get_nl <- function(x, y)
  {
    lin_pred <- predict(lm(y ~ x))
    gam_pred <- predict(mgcv::gam(y ~ s(x)))
    L <- sqrt(sum((lin_pred-gam_pred)^2)/length(x)) ## RMSE calculation
    L
  }
  nl_up <- get_nl(out$a[1:1000], out$x[1:1000])
  nl_down <- get_nl(out$a[1001:2000], out$x[1001:2000])
  
  ## hysteresis
  out2 <- data.frame(up_x=out$x[1:1000],
                   down_x=out$x[2000:1001])
  hyst = mean(abs(out2$up_x - out2$down_x))
  
  return(list(nl_up=nl_up, nl_down=nl_down, hyst=hyst))
}
#ps <- c(seq(1, 5, by=1), 10, 20, 50)
ps <- 10^seq(0, log10(10), length=20)
inv_rates <- c(200, 20000)
#rezs <- lapply(ps, function(x) get_nl_hyst_by_p(x))

expt3 <- expand.grid(ps=ps,
                     inv_rates=inv_rates)

## Here a starts low, increases, then decreases
a_forcing <- c(0.001,1,0.001)

temp_res3 <- expt3 %>%
  group_by(ps, inv_rates) %>%
  do(rez = as.data.frame(get_nl_hyst_by_p(p=.$ps, inv_rate = .$inv_rates, a_forcing)))

results3 <- temp_res3 %>%
  tidyr::unnest() %>%
  gather(key=Variable, value=Value, 3:5) %>%
  mutate(Variable=case_when(Variable=="nl_up" ~ "Non-linearity (up phase)",
                            Variable=="nl_down" ~ "Non-linearity (down phase)",
                            Variable=="hyst" ~ "Hysteresis"),
         `Simulation duration`=as.character(1/inv_rates))
## Warning: `cols` is now required.
## Please use `cols = c(rez)`
ggplot(results3, aes(x=ps, y=Value, col=`Simulation duration`)) +
  facet_wrap(~Variable, scales = "free") +
  geom_line() + geom_point() +
  scale_x_log10() +
  xlab("Strength of positive feedback") +
  guides(colour=guide_legend(title="Rate of\nenvironmental\nchange"))

#ps <- c(seq(1, 5, by=1), 10, 20, 50)
ps <- 10^seq(0, log10(10), length=20)
inv_rates <- c(200, 20000)
#rezs <- lapply(ps, function(x) get_nl_hyst_by_p(x))

expt3 <- expand.grid(ps=ps,
                     inv_rates=inv_rates)

## Here a starts high, decreases, then increases
a_forcing <- c(1,0.001,1)

temp_res3 <- expt3 %>%
  group_by(ps, inv_rates) %>%
  do(rez = as.data.frame(get_nl_hyst_by_p(p=.$ps, inv_rate = .$inv_rates, a_forcing)))

results3 <- temp_res3 %>%
  tidyr::unnest() %>%
  gather(key=Variable, value=Value, 3:5) %>%
  mutate(Variable=case_when(Variable=="nl_up" ~ "Non-linearity (down phase)",
                            Variable=="nl_down" ~ "Non-linearity (up phase)",
                            Variable=="hyst" ~ "Hysteresis"),
         `Simulation duration`=as.character(1/inv_rates))
## Warning: `cols` is now required.
## Please use `cols = c(rez)`
ggplot(results3, aes(x=ps, y=Value, col=`Simulation duration`)) +
  facet_wrap(~Variable, scales = "free") +
  geom_line() + geom_point() +
  scale_x_log10() +
  xlab("Strength of positive feedback") +
  guides(colour=guide_legend(title="Rate of\nenvironmental\nchange"))