System analysis: Alternative steady-states and hysteresis

David Kneis (firstname.lastname @ tu-dresden.de)

Intro

  • For given parameters and forcings, we usually expect a system to converge to a unique steady state, independent of the system’s initial state.

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

  • We use an example from Marten Scheffer’s book Ecology of Shallow Lakes, p. 107/108 to demonstrate hysteresis.

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

Load the model

  • The model is part of the rodeoEasy package. You may open the respective workbook on your local computer.
rm(list=ls())                                     # initial clean-up

#remotes::install_github("dkneis/rodeoEasy")      # if not installed yet
library("rodeoEasy")

workbook <- system.file("models/hysteresis.xlsx", # contains the model
  package="rodeoEasy")

m <- build(workbook)                              # build model object

State variables

name unit description default
G mg C / L green algae 10.0
B mg C / L bluegreen algae 10.0
DIP mg P / L dissolved inorganic phosphorus 0.1

Model equations

name unit description rate G B DIP
grw_G mg C / L / h growth of green algae rg * G * hsg / (hsg + D * (eg * G + eb * B)) * DIP / (DIP + hpg) 1 0 -1/cpg
grw_B mg C / L / h growth of blue-green algae rb * B * hsb / (hsb + D * (eg * G + eb * B)) * DIP / (DIP + hpb) 0 1 -1/cpb
loss_G mg C / L / h loss of green algae xg * G -1 0 1/cpg
loss_B mg C / L / h loss of blue-green algae xb * B 0 -1 1/cpb
flush mg / L / h flushing 1/tau -G -B DIPin - DIP

Monod vs. inverse monod

  • Phosphorus limitation is described by a classical Monod term: More P promotes growth.

\(\frac{DIP}{DIP + h}\)

  • Shadow (\(S\)) is expressed as depth times turbidity. Tolerance to shadow (\(S\)) is modeled using an inverse Monod term since more shadow inhibits growth.

\(1 - \frac{S}{S + h} = \frac{S + h}{S + h} - \frac{S}{S + h} = \frac{h}{h + S}\)

Parameters

Differences between the two groups of algae are reflected by different 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).

Parameters, part 1

name unit description default
rg 1 / h growth rate constant of green algae 1.2
rb 1 / h growth rate constant of blue-green algae 0.6
xg 1 / h loss rate constant of green algae 0.12
xb 1 / h loss rate constant of blue-green algae 0.06
hsg dimenionless controls shade tolerance of green algae 1.5
hsb dimenionless controls shade tolerance of blue-green algae 3
hpg mg P / L controls P limitation of green algae 0.01
hpb mg P / L controls P limitation of blue-green algae 0.01
eg 1 / m / (mg C / L) specific light attenuation const. for green algae 0.5
eb 1 / m / (mg C / L) specific light attenuation const. for blue-green algae 1

Parameters, part 2

name unit description default
cpg mass ratio C:P ratio in green algae 41.03
cpb mass ratio C:P ratio in blue-green algae 41.03
tau days residence time 7
DIPin mg P / L phosphorus in inflow 0.02
D m water depth 3

Considered scenarios

  • We try two distinct initial value scenarios: oligotrophic vs. hypertrophic. This represents the system’s history.

  • For each initial value scenario, we consider a gradient of phosphorus inputs.

Creation of contrasting initial states

scen <- list(
  oligo = c(DIPin = 0.005),
  hyper = c(DIPin = 0.22)
)
x <- run.scenarios(m, 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 small inoculum

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))            # same 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 <- run.scenarios(m, times=c(0, 1e7),
  scenarios=scen, plot.vars=FALSE)

Compile results for plotting

  • We put the 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  3.157966e-39 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)

Visualize the hysteresis

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 <- run.scenarios(m, 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)

State space plot, ctd.

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 system and make it switch back to a green algae state, we need to decrease P input beyond the intuitive threshold (hysteresis).

  • In reality, further mechanisms impede lake restoration, most notably internal phosphorus loading.

Summary 2

  • Hysteresis effects also occur 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.