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 measurements of, e.g., phyto- and zooplankton densities and trophic pyramids in general.

  • Predation is typically but not necessarily linked to cyclic dynamics, i.e. oscillations of species abundances.

About the example

  • We consider predation of the heterotrophic nano-flagellate Poteriospumella lacustris on the bacterium Pseudomonas putida.

  • Several parameters were adopted from a case study.

  • The example does not consider predation defense or predator-prey co-evolution for the sake of simplicity.

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/predation.xlsx",  # contains the model
  package="rodeoEasy")

m <- build(workbook)                              # build model object

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

name unit description rate R B F
grB cells/mL/h growth of bacteria muB * (R / (R+hR)) * B -R_per_B 1 0
grF cells/mL/h growth of flagellates muF * max(0, (B-minB) / (B-minB+hB)) * F 0 -B_per_F 1
exR µg/mL/h export of resources 1/tau * R -1 0 0
exB cells/mL/h export of bacteria 1/tau * B 0 -1 0
exF cells/mL/h export of flagellates 1/tau * F 0 0 -1
imR µg/mL/h import of resources 1/tau * Rin 1 0 0

Parameters

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 resource consumed per bacterial reproduction 1.0e-06
B_per_F cells / cell bacteria 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 abundance for flagellate growth 1.0e+03

Functions

name code
max # intrinsic function

Code for convenient plotting

  • State variables and legend arranged in 2 x 2 layout.

  • Axes with log-scaling and pretty labels.

custom.plot <- function(x) {
  oma <- par("mar")
  par(mar=c(4,4,1,1))
  par(mfrow=c(2,2), cex=1.3)
  cc <- show.scenarios(x, "B", ylog=TRUE)
  cc <- show.scenarios(x, "F", ylog=TRUE)
  cc <- show.scenarios(x, "R", ylog=TRUE)
  plot(0, 0, type="n", ann=F, axes=F)
  legend("center", bty="n", lty=1, col=cc,
    legend=names(cc), 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 <- run.scenarios(m, times=0:(20*24),
  scenarios=scen, plot.vars=F)
custom.plot(x)

Scenarios 1: Impact of resources

Interpretation

Very low: Even bacteria go extinct

Low: Predator does not survive

Moderate: Classic steady-state

High & beyond: Stable oscillations

Scenarios 2: Capabilities of predator

scen <- list(
  `very low` =c(hB=1e4),
  `low`      =c(hB=1e5),
  `moderate` =c(hB=1e6),
  `high`     =c(hB=3e6),
  `very high`=c(hB=4e6)
)
x <- run.scenarios(m, times=0:(60*24),
  scenarios=scen, plot.vars=F)
custom.plot(x)

Scenarios 2: Capabilities of predator

Attactor plot for default settings

# long-term simulation
x <- run.scenarios(m, 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