System analysis: Predation

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

Intro

  • Predation is a major cause of mortality.

  • An understanding of predation is necessary to interpret simultaneous data of, e.g., phyto- and zooplankton and trophic pyramids in general.

  • Predation may be linked to oscillations of species abundances (see classical examples).

About this example

Load the model

  • The model is is provided as a workbook.
rm(list=ls())                              # initial clean-up
library("rodeo")                           # load package

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

State variables

name unit description default
R µg / mL resource for bacterial growth 1e+00
B cells / mL abundance of prey bacteria 1e+06
F cells / mL abundance of flagellates 1e+03

Model equations

description rate R B F
growth of bacteria muB * (R / (R+hR)) * B -R_per_B 1 0
growth of flagellates muF * max(0, (B-minB) / (B-minB+hB)) * F 0 -B_per_F 1
export of resources 1/tau * R -1 0 0
export of bacteria 1/tau * B 0 -1 0
export of flagellates 1/tau * F 0 0 -1
import of resources 1/tau * Rin 1 0 0

Parameters and functions

name unit description default
muB 1 / h max. growth rate constant of bacteria 1.0e+00
muF 1 / h max. growth rate constant of flagellates 1.0e-01
hR µg / mL half saturation const. in bacterial growth 5.0e-01
hB cells / mL half saturation const. in flagellate growth 1.0e+05
R_per_B µg / cell R consumed in bacterial reproduction 1.0e-06
B_per_F cells / cell B consumed per flagellate reproduction 3.0e+02
Rin µg / mL resource concentration in inflow 1.0e+00
tau h residence time 4.8e+01
minB cells / mL min. bacterial density for flagellate growth 1.0e+03
max - intrinsic function NA

Code for convenient plotting (1)

Comparison of scenarios for a particular variable

color.scen <- function(scen) { c(`very low`="steelblue4",
  `low`="steelblue1", `moderate`="olivedrab", `high`="orange",
  `very high`="firebrick")[scen] }

plot.scenarios <- function(x, varname) {
  plot(range(x[,"time"]), range(x[,varname]), type="n",
    ann=F, yaxt="n", log="y")
  mtext(side=3, varname, cex=par("cex"))
  axis(2, las=2)
  fn <- function(x) {
    lines(x[,"time"], x[,varname], col=color.scen(x[,"scenario"]))
  }
  by(x, x[,"scenario"], fn)
}

Code for convenient plotting (2)

State variables and legend arranged in 2 x 2 layout

plot.all <- function(x) {
  oma <- par("mar")
  par(mar=c(4,4,1,1))
  par(mfrow=c(2,2), cex=1.3)
  plot.scenarios(x, "B")
  plot.scenarios(x, "F")
  plot.scenarios(x, "R")
  plot(0, 0, type="n", ann=F, axes=F)
  legend("center", bty="n", lty=1,
    col=color.scen(unique(x[,"scenario"])),
    legend=unique(x[,"scenario"]), title="Scenario")
  par(mfrow=c(1,1), cex=1)
  par(mar=oma)
}

Scenarios 1: Impact of resources

scen <- list(
  `very low` =c(Rin=0.001),
  `low`      =c(Rin=0.01),
  `moderate` =c(Rin=.1),
  `high`     =c(Rin=1),
  `very high`=c(Rin=10)
)
x <- m$scenarios(times=0:(20*24), scenarios=scen, plot.vars=F)
plot.all(x)

Scenarios 1: Impact of resources

Interpretation

Very low R supply: Even bacteria go extinct

Interpretation

Low R supply: Predator still does not survive

Interpretation

Moderate R supply: Classic steady-state

Interpretation

High & beyond: Stable oscillations

Scenarios 2: Density dependence

scen <- list(
  `very low` =c(hB=1e4),  # hB: Density of bacteria supporting
  `low`      =c(hB=1e5),  #     growth of the predator at 50%
  `moderate` =c(hB=1e6),  #     of the possible maximum rate
  `high`     =c(hB=3e6),
  `very high`=c(hB=4e6)
)
x <- m$scenarios(times=0:(60*24), scenarios=scen, plot.vars=F)
plot.all(x)

Scenarios 2: Density dependence

Attactor plot for default settings

# long-term simulation
x <- m$scenarios(times=c(0, 10^seq(0, 4, 0.01)), plot.vars=FALSE)

# state-space plot with color coding for time
col <- colorRampPalette(c("cyan","steelblue"))(nrow(x))

omar <- par("mar"); par(mar=c(4,4,1,1))
plot(x$B, x$F, ylim=range(x$F)*c(0.7,1),
  type="p", pch=20, col=col, log="xy",
  xlab="Bacteria (cells/mL)", ylab="Flagellates (cells/mL)")
legend("bottomleft", bty="n", pch=20, col=c("cyan","steelblue"),
  legend=c("Initial state", "Steady state"), horiz=TRUE)
par(mar=omar)

  • The attractor is no longer a fixed point.

  • The closed trajectory corresponding to the steady state is known as a limit cycle.

Attactor plot for default settings