System analysis: Alternative steady-states and hysteresis

Intro

• For given parameters and forcings, many system evolve towards a unique steady state, independent of the system’s initial state.

• This is not necessarily so. If the (steady) state exhibits a particular dependence on history, this is called hysteresis.

• Marten Scheffer’s book Ecology of Shallow Lakes (p. 107/108) provides a nice example.

• The example considers green algae and cyanobacteria competing for light and phosphorus.

Load the model

• The model is is provided as a workbook.

rm(list=ls())                              # initial clean-up
library("rodeo")                           # load package

workbook <- "models/hysteresis.xlsx"       # adjust path as necessary
m <- buildFromWorkbook(workbook)           # build model

Resource dependence (1)

• Phosphorus limitation is described by a classical Monod term

• More P promotes growth.

\(\mu = \mu_{max} \cdot \frac{DIP}{DIP + h}\)

Resource dependence (2)

• Growth is modeled as a function of shadow rather than light.

• Tolerance to shadow (\(S\)) is expressed by an inverse Monod term, since more shadow inhibits growth.

\(\mu = \mu_{max} \cdot \left( 1 - \frac{S}{S + h} \right) = \mu_{max} \cdot \frac{h}{h + S}\)

• Shadow is approximated by the product of turbidity and depth. Turbidity is controlled by algae density.

Parameters

Differences between algae groups are reflected by contrasting parameter values.

1. Green algae have a higher intrinsic growth rate than cyanobacteria.

2. Green algae suffer from higher mortality (e.g. due to grazing and settling).

3. Cyanobacteria have a higher shade tolerance than green algae but also cause more shade (represented by a higher specific extinction coefficient).

Considered scenarios

• We study two distinct initial value scenarios representing different system history:

  • oligotrophic steady state
  • hypertrophic steady state

• The two initial states are confronted with a gradient of phosphorus inputs to model eutrophication or re-oligotrophication, respectively.

Creation of initial states

scen <- list(
  oligo = c(DIPin = 0.005),
  hyper = c(DIPin = 0.22)
)
x <- m$scenarios(times=c(0, 1e7),        # long-term sim.
  scenarios=scen, plot.vars=F)

x <- x[x[,"time"] == max(x[,"time"]),]   # final values as matrix
rownames(x) <- x[,"scenario"]
x <- as.matrix(x[,c("G","B")])

x[x <= 1e-6] <- 1e-6                     # ensure a small inoculum of
                                         # inferior competitor

print(x)

               G        B
oligo 0.07827241 0.000001
hyper 0.00000100 1.798862

Scenarios with different forcing

# scenarios as a data frame
scen <- rbind(
  data.frame(as.list(x["oligo",c("G","B")]),  # clear initial state
    DIPin=seq(0.005, 0.22, 0.005)),           # gradient of P supply
  data.frame(as.list(x["hyper",c("G","B")]),  # turbid initial state
    DIPin=seq(0.005, 0.22, 0.005))            # identical gradient
)

# turn scenario data frame into suitable list
scen <- apply(scen, 1, unlist, simplify=FALSE)
names(scen) <- paste0("s",1:length(scen))

# simulation of steady states
x <- m$scenarios(times=c(0, 1e7),
  scenarios=scen, plot.vars=FALSE)

Compile results for plotting

• We put initial and steady-state values side by side

x <- data.frame(
  G.std=x[x[,"time"] > 0, "G"],
  B.std=x[x[,"time"] > 0, "B"],
  G.ini=x[x[,"time"] == 0, "G"],
  B.ini=x[x[,"time"] == 0, "B"],
  DIPin=sapply(scen, function(z) {z[["DIPin"]]})
)

print(head(x))

        G.std         B.std      G.ini B.ini DIPin
s1 0.07827241  8.330422e-44 0.07827241 1e-06 0.005
s2 0.25478077  4.774641e-43 0.07827241 1e-06 0.010
s3 0.42857921  7.740739e-43 0.07827241 1e-06 0.015
s4 0.59934579 -2.019687e-40 0.07827241 1e-06 0.020
s5 0.76672276  1.635046e-41 0.07827241 1e-06 0.025
s6 0.93031616  9.373966e-40 0.07827241 1e-06 0.030

Visualize the hysteresis

omar <- par("mar"); par(mar=c(4,4,1,1))
with(x, {
  plot(range(DIPin), c(-0.1, max(G.std,B.std)), type="n",
    xlab="Phosphorus input (DIPin) in mg / L",
    ylab="Algae in steady state (mg C / L)")
  points(DIPin, G.std, col="olivedrab",
    pch=ifelse(G.ini > B.ini, 3, 0))
  points(DIPin, B.std, col="steelblue4",
    pch=ifelse(G.ini > B.ini, 4, 5))
})

legend(x=0.11, y=1.5, bty="n", title="Initially clear state",
  col=c("olivedrab", "steelblue4"), pch=c(3, 4),
  legend=c("Green algae", "Blue-greens"))
legend(x=0.11, y=0.8, bty="n", title="Initially turbid state",
  col=c("steelblue4", "olivedrab"), pch=c(0, 5),
  legend=c("Blue-greens", "Green algae"))
par(mar=omar)

Lessons from the hysteresis plot

• Steady states are exclusive: Green algae dominate at low P input, otherwise cyanobacteria take over

• There is a range of P inputs where the system can take one of two alternative steady states. This is known as bi-stability.

• Within the range of bi-stability, real-world ecosystems tend to be sensitive to perturbations (switch between alternatives).

• The predicted steady state obviously depends on the system’s history (represented by the initial state).

State space plot

# initial state scenarios as data frame
scen <- expand.grid(
  DIPin=0.065,            # fixed within the critical range
  G=c(0.1, 0.5, 1, 2, 3),
  B=c(0.1, 0.5, 1, 2, 3)
)

# cast scenarios as suitable list
scen <- apply(scen, 1, unlist, simplify=FALSE)
names(scen) <- paste0("s",1:length(scen))

# long-term simulation with sparse output
x <- m$scenarios(times=c(0, 10^seq(0,7,0.1)),
  scenarios=scen, plot.vars=FALSE)

State space plot, ctd.

# start with empty plot
omar <- par("mar"); par(mar=c(4,4,2,1))
with(x, plot(range(G), range(B), type="n",
    xlab="Green algae", ylab="Blue-greens"))

# draw trajectories, end points highlighted
fn <- function(x) {
  last <- nrow(x)
  lines(x$G, x$B, col="grey")
  points(x$G[last], x$B[last], pch=20, cex=2,
    col=if(x$G[last] > x$B[last]) "olivedrab" else "steelblue4")
  points(x$G[1], x$B[1], col="grey")
}
by(x, x[,"scenario"], fn)

legend(0, 3.5, xpd=TRUE, bty="n", horiz=TRUE, pch=c(1,20),
  pt.cex=1:2, legend=c("Initial states", "Steady-states"))
par(mar=omar)

Summary 1

• Eutrophication (increased P input) will, at some point, leads to a re-structuring of the phytoplankton community from green algae to cyanobacteria.

• More on critical thresholds can be found in this paper by M. Scheffer et al..

• To restore a cyanobacteria-dominated lake and make it switch back to a green algae state, we must reduce P input beyond the intuitive threshold (hysteresis).

• In reality, further mechanisms impede lake restoration (e.g. internal phosphorus loading).

Summary 2

• Hysteresis also occurs in other contexts. The phenomenon is not necessarily limited to steady states.

• For example, there is no unique relationship between flow rate (Q) and water level (H).

• For a given value of Q, the corresponding H is lower during the onset of a flood wave while it is higher when water levels are falling.

• This results in loop-shaped rating curves.