Dynamical Systems in Biology

Exercise 1: Enzyme dynamics

Implement the enzyme reaction

\[ \begin{align*} S + E {{k_{ {\scriptscriptstyle}1} \atop \longrightarrow} \atop {\longleftarrow \atop k_{{\scriptscriptstyle}-1}}} C \stackrel{k_{2}}{\longrightarrow} P + E. \end{align*} \]

using mass action kinetics:

\[ \begin{align} \dot{s} & = k_{{\scriptscriptstyle}-1} \cdot c - k_{{\scriptscriptstyle}1} s \cdot e \\ \dot{e} & = (k_{{\scriptscriptstyle}-1} + k_{{\scriptscriptstyle}2}) \cdot c - k_{{\scriptscriptstyle}1} \cdot s \cdot e \\ \dot{c} & = k_{{\scriptscriptstyle}1} \cdot s \cdot e - (k_{{\scriptscriptstyle}-1} + k_{{\scriptscriptstyle}2}) \cdot c \\ \dot{p} & = k_{{\scriptscriptstyle}2} \cdot c, \end{align} \]

fun1 <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    
    dS <- koff*C-kon*S*E
    dP <- k2*C
    dE <- (koff+k2)*C - kon*S*E
    dC <- kon*S*E-(koff+k2)*C
    
    
    return(list(c(dS, dP, dE, dC)))
  })
}

using the steady-state approximation:

\[ \begin{align} \dot{s} & = - \dot{p} = - \frac{k_{{\scriptscriptstyle}2} \mathrm{e}_{{\scriptscriptstyle}T} \cdot s}{\frac{k_{{\scriptscriptstyle}-1}}{k_{{\scriptscriptstyle}1}} + s} \\ c & = \frac{\mathrm{e}_{{\scriptscriptstyle}T} \cdot s}{\frac{k_{{\scriptscriptstyle}-1}}{k_{{\scriptscriptstyle}1}} + s}. \end{align} \]

fun2 <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    
    dS <- -(k2*E*S)/(S+koff/kon)
    dP <- -dS
    dE <- 0
    
    C <- E*S/(S + koff/kon)
    
    return(list(c(dS, dP, dE), C = C))
  })
}

using the quasi-steady state approximation:

\[ \begin{align} \dot{s} & = - \dot{p} = - \frac{k_{{\scriptscriptstyle}2} \mathrm{e}_{{\scriptscriptstyle}T} \cdot s}{\frac{k_{{\scriptscriptstyle}-1}}{k_{{\scriptscriptstyle}1}} + s} \\ c & = \frac{\mathrm{e}_{{\scriptscriptstyle}T} \cdot s}{K_{{\scriptscriptstyle}m} + s} \quad \text{mit} \quad K_{{\scriptscriptstyle}m} = \frac{k_{{\scriptscriptstyle}-1} + k_{{\scriptscriptstyle}2}}{k_{{\scriptscriptstyle}1}} \end{align} \]

fun3 <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    
    dS <- -(k2*E*S)/(S+koff/kon)
    dP <- -dS
    dE <- 0
    
    Km <- (koff + k2)/kon
    C <- E*S/(S + Km)
    
    return(list(c(dS, dP, dE), C = C))
  })
}

Simulate the three dynamic systems with parameters \(k_1 = 0.1\), \(k_{-1} = 0.1\), and \(k_2 = 0.1\). Let the initial concentrations be \(e(0) = 10\), \(c(0) = 0\), and \(s(0) = 10, 20, ..., 100\). Hint: The return value of the ODE function is a list, the first entry being the derivatives. The list may contain more named elements, e.g., to return values of algebraic expressions of the ODE states and parameters.

library(deSolve)
library(tidyr)
library(ggplot2)

pars <- c(kon = .1,
          koff = .1,
          k2 = .1)

   # here for S=50:
yini  <- c(S = 50,
           P = 0,
           E = 10,
           C = 0)

times <- seq(0, 100, by = 0.1)

## Different dynamics
out1 <- ode(yini, times, fun1, pars)
out2 <- ode(yini[c("S", "P", "E")], times, fun2, pars)
out3 <- ode(yini[c("S", "P", "E")], times, fun3, pars)

# Put everything in one data.frame
# After 4 tutorials of writing your own gather functions, have a look at
# gather from the tidyr package.
out <- rbind(
  cbind(gather(as.data.frame(out1), key = "species", value = "value", -time), method = "MA"),
  cbind(gather(as.data.frame(out2), key = "species", value = "value", -time), method = "SS"),
  cbind(gather(as.data.frame(out3), key = "species", value = "value", -time), method = "QSS")
)


ggplot(out, aes(x = time, y = value, color = method)) + facet_wrap(~species, scales = "free") +
  geom_line()

# now several s(0) values:
svalues <- c(10,25,50,75,100,500)

out <- do.call(rbind, lapply(svalues, function(sval){
  yini  <- c(S = sval,
           P = 0,
           E = 10,
           C = 0)
  out1 <- ode(yini, times, fun1, pars)
  out2 <- ode(yini[c("S", "P", "E")], times, fun2, pars)
  out3 <- ode(yini[c("S", "P", "E")], times, fun3, pars)
  out <- rbind(
    cbind(gather(as.data.frame(out1), key = "species", value = "value", -time), method = "MA", Sini = as.character(sval)),
    cbind(gather(as.data.frame(out2), key = "species", value = "value", -time), method = "SS", Sini = as.character(sval)),
    cbind(gather(as.data.frame(out3), key = "species", value = "value", -time), method = "QSS", Sini = as.character(sval))
  )
  return(out)
  })
)

ggplot(out, aes(x = time, y = value, color = method)) + facet_grid(Sini~species, scales = "free") + geom_line()

Take a closer look at the initial time frame of the full system implemented by mass-action kinetics. What do you observe?
Have a look at the dynamics of all implementations in phase space \(s\) vs. \(c\).

outPS <- rbind(
  cbind(as.data.frame(out1)[ , c("time", "S", "C")], method = "MA"),
  cbind(as.data.frame(out2)[ , c("time", "S", "C")], method = "SS"),
  cbind(as.data.frame(out3)[ , c("time", "S", "C")], method = "QSS")
)


ggplot(outPS, aes(x = S, y = C, color = method)) +
  geom_path()

Dynamical Systems in Biology

Tutorial 6

Daniel Kaschek, Mirjam Fehling-Kaschek

Summer term 2017

Exercise 1: Enzyme dynamics

Cathedral exercise: