Dynamical Systems in Biology

Exercise 1: Enzyme reaction chain

Consider a chain of Michaelis-Menten enzyme reactions:

\[ S \stackrel{E_{1}}{\longrightarrow} S_{1} \stackrel{E_{2}}{\longrightarrow} S_{2} \stackrel{E_{3}}{\longrightarrow} S_{3} \stackrel{E_{4}}{\longrightarrow} P \]

for a constant concentration \(S=1\) and given \(V_\mathrm{max}\) and \(K_{\mathrm{M}}\) values:

\[ \begin{align*} & &V_{{\scriptscriptstyle}\mathrm{max}} &\quad K_{{\scriptscriptstyle}\mathrm{M}} \\ & E_{{\scriptscriptstyle}1} &0.1 &\quad 0.1 \\ & E_{{\scriptscriptstyle}2} &1.0 &\quad 1.0 \\ & E_{{\scriptscriptstyle}3} &1.0 &\quad 0.1 \\ & E_{{\scriptscriptstyle}4} &5.0 &\quad 5.0 \\ \end{align*} \]

Implement the system and plot the solution in config-space.

## enzyme reaction chain
fn <- function(time, state, pars) {
  with(as.list(c(state, pars)), {
  
    dS1 <- v0*S /(K0 + S ) - v1*S1/(K1 + S1)
    dS2 <- v1*S1/(K1 + S1) - v2*S2/(K2 + S2)
    dS3 <- v2*S2/(K2 + S2) - v3*S3/(K3 + S3)
    dP <- v3*S3/(K3 + S3)
    
    return(list(c(dS1, dS2, dS3, dP)))
  })  
  
}

yini <- c(S1 = 0, S2 = 0, S3 = 0, P = 0)
pars <- c(S=1, K0 = 0.1, K1 = 1, K2 = 0.1, K3 = 5, 
          v0 = 0.1, v1 = 1.0, v2 = 1.0, v3 = 5.0)

times <- seq(0, 10, len=200)

result <- ode(yini, times, fn, pars)

out <- gather(as.data.frame(result), key = "species", value = "value", -time)

ggplot(out, aes(x = time, y = value, color = species)) + geom_line()

Determine the steady state concentrations \(S_{{\scriptscriptstyle}1}, S_{{\scriptscriptstyle}2}\) und \(S_{{\scriptscriptstyle}3}\) and the steady state flux \(J\) by simulating the system for long time periods.

# system in steady state for time < 10min
times <- seq(0, 20, len=200)
result <- ode(yini, times, fn, pars)[length(times),]
print(round(result[c("S1","S2","S3")], digits = 3))

##    S1    S2    S3 
## 0.100 0.010 0.093

# get flux by evaluating fn (get dP value)
fluxes <- unlist(fn(200, result[-1], pars))
print(fluxes[4])

## [1] 0.09090905

Have a look at the function rootSolve::steady() and use it to calculate the steady state concentrations and flux.

library(rootSolve)

## enzyme reaction chain without dP
fnSS <- function(time, state, pars) {
  with(as.list(c(state, pars)), {
    dS1 <- v0*S /(K0 + S ) - v1*S1/(K1 + S1)
    dS2 <- v1*S1/(K1 + S1) - v2*S2/(K2 + S2)
    dS3 <- v2*S2/(K2 + S2) - v3*S3/(K3 + S3)
    return(list(c(dS1, dS2, dS3)))
  })  
}

yini <- c(S1 = 0, S2 = 0, S3 = 0)

# steady state flux
J <- function(state, pars) {
  with(as.list(c(state, pars)), v3*S3/(K3 + S3))
}

# steady state concentrations and flux
SJ <- function(pars) {
  S <- steady(yini, func=fnSS, parms=c(pars))$y
  J <- J(S, c(pars))
  return(c(S, J=J))
}

print(SJ(pars))

##          S1          S2          S3           J 
## 0.099999995 0.009999999 0.092592593 0.090909091

Evaluate the steady state concentrations and flux by varying the \(V_{\mathrm{max}}\) values. Use the manipulate package for an interactive visualization and plot the original values along.

library(manipulate)

# original steady state values
ssout0 <- SJ(pars)

# remove v0, ..., v3 from parameter vector
parsGlobal <- pars[1:5]

# call the interactive plot with sliders for v0, ..., v3
manipulate({
  pars <- c(parsGlobal, v0=v0, v1=v1, v2=v2, v3=v3)
  ssout <- SJ(pars)
  data <- rbind(data.frame(state = names(ssout), value = as.numeric(ssout), pars="adjusted"),
                data.frame(state = names(ssout0), value = as.numeric(ssout0), pars="original"))
  
  ratios <- format(ssout/ssout0, digits=2, nsmall=2)
  
  P <- ggplot(data, aes(x=state, y=value, group=pars, fill=pars)) + 
    geom_bar(stat="identity", position="dodge", color="black") + 
    scale_fill_manual(values = c("black", "white")) +
    annotate("text", y=rep(0.3, length(ratios)), x=1:4, label=ratios) +
    ylim(c(0, 0.3)) 
  print(P)
}, 
v0 = slider(0, 0.5, 0.1, step=0.05),
v1 = slider(0, 2, 1, step=0.1),
v2 = slider(0, 2, 1, step=0.1),
v3 = slider(0, 10, 5, step=1)
)

Calculate the control coefficients using the numDeriv package. Hint: have a look at outer() for calculating the outer product of two vectors.

library(numDeriv)

# calculate the control coefficients
parsV <- pars[grepl("v", names(pars))]
CkSJ <- (outer(1/SJ(pars), pars, "*"))*jacobian(SJ, pars)
print(CkSJ)

##             S          K0           K1            K2            K3
## S1 0.09999988 -0.10000003 9.999999e-01 -6.094092e-08  7.290394e-08
## S2 0.09999997 -0.09999997 0.000000e+00  1.000000e+00  0.000000e+00
## S3 0.09259259 -0.09259259 8.243829e-11 -4.127238e-10  1.000000e+00
## J  0.09090909 -0.09090909 8.106958e-11 -4.056171e-10 -1.698710e-10
##          v0            v1            v2            v3
## S1 1.100000 -1.100000e+00  6.094092e-08  0.000000e+00
## S2 1.100000  0.000000e+00 -1.100000e+00  0.000000e+00
## S3 1.018519  7.225437e-12  4.127088e-10 -1.018519e+00
## J  1.000000  7.359780e-12  4.056023e-10  1.540729e-10

# plot the control coefficients
plotCkSJ <- function(M) {
  nrow <- dim(M)[1]
  data <- data.frame(state = rownames(M), pars = rep(colnames(M), each=nrow), coefficient = as.numeric(M))
  ggplot(data, aes(x=state, y=coefficient, group=pars, fill=pars, color=pars)) + geom_bar(stat="identity", position="dodge")
}

print(plotCkSJ(CkSJ))

Dynamical Systems in Biology

Tutorial 8

Daniel Kaschek, Mirjam Fehling-Kaschek

Summer term 2017

Exercise 1: Enzyme reaction chain

Cathedral exercise: