Dynamical Systems in Biology

Tutorial 9

Exercise 1: Positive cooperative feedback

Consider a protein \(A\) being phosphorylated by a stimulus \(S\) with rate \(k_1\) and dephosphorylated with rate \(k_{-1}\). The phosphorylation is cooperatively enhanced by a positiv feedback with hill coefficient \(k=4\) and rate \(k_2\):

\[ \begin{align} \frac{\mathrm{d}}{\mathrm{d}t} pA= k_{1} \cdot S \cdot A + k_{2} \cdot A \cdot \frac{pA^4}{K_{m}^4 + pA^4} - k_{-1}\cdot pA. \end{align} \]

Further assume that the total amount of protein equals 1 (\(pA + A = A_{\mathrm{total}} = 1\)):

\[ \begin{align*} \frac{\mathrm{d}}{\mathrm{d}t} pA = \left( k_{1} \cdot S + k_{2} \cdot \frac{pA^4}{K_m^4 + pA^4}\right) \cdot (1-pA) - k_{ -1}\cdot pA. \end{align*} \]

The parameters are given by \(k_{1} = 0.1\), \(k_{-1} = 1\), \(k_{2} = 2\) and \(K_{m} = 0.3\).

• Implement the system and plot the solution in config-space for \(S=1\). Compare the solution with the solution without a positive feedback.

## phosphorylation with positive feedback 
fn <- function(time, state, pars) {
    with(as.list(c(state, pars)), {
      dpA <- (k1*S+k2*pA^4/(Km^4 + pA^4))*(1-pA) - k_1*pA
      return(list(c(dpA), A=1-pA))
    })
}

pars <- c(S = 1, k1 = 0.1, k_1 = 1, k2 = 2, Km = 0.3)
yini <- c(pA = 0)
times <- seq(0, 20, len=500)

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

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

parsNoF <- c(S = 1, k1 = 0.1, k_1 = 1, k2 = 0, Km = 0.3)

resultNoF <- ode(yini, times, fn, parsNoF)

outNoF <- cbind(gather(as.data.frame(resultNoF), key = "species", value = "value", -time), feedback="no")

ggplot(rbind(out, outNoF), aes(x = time, y = value, color = species, lty = feedback)) + geom_line()

• Stimulate the system with \(S\in[0,2]\) starting with \(S = 0\), increasing \(S\) in small steps until \(S = 2\) and then reducing \(S\) in small steps until \(S = 0\) again. In each step, simulate or compute the steady state of \(pA\) for each value of \(S\) either by simulating the system for long times or by rootSolve::steady() (see last exercise) and take its value as initial value for the next \(S\) value. Start with \(pA=0\) for \(S=0\).
  • How is the observed phenomenon called in physics?
  • What is its meaning in biological systems?

times <- seq(0, 100, length.out = 500)

SSvsS <- function(S = seq(0, 2, by = 0.01), yini = c(pA=0)) {
    SS <- NULL
    for(s in S) {
    
        p <- c(S = s, pars[-1])
        result <- ode(yini, times, fn, p)
    yini <- result[length(times),"pA"]
    
        SS <- c(SS, yini)
    
    }
    return(data.frame(S=S, pA = SS))
}

outforth <- SSvsS()
outforth$way <- "forth"

yini <- c(pA = tail(outforth,1)$pA)
outback <- SSvsS(S = seq(2, 0, by = -0.01), yini = yini)
outback$way <- "back"

SSvsSplot <- ggplot(rbind(outforth, outback), aes(x = S, y = pA, color = way)) + geom_path()
print(SSvsSplot)

• Generate a rate balance plot to illustrate the mode of action: plot the build and decay fluxes of \(pA\) as a function of \(pA\).
  • How can you extract information on the qualitative dynamics, especially the stability of the fixed points, from this plot?

## flux calculations 
flux <- function(state, pars) {
    with(as.list(c(state, pars)), {
      build <- (k1*S+k2*pA^4/(Km^4 + pA^4))*(1-pA)
      decay <- k_1*pA
      dpA <- build - decay
      return(c(dpA = dpA, build = build, decay = decay, S=S, pA = pA))
    })
}

S <- seq(0, 2, len=9)
pA <- seq(0, 1, len=100)

SpA <- expand.grid(S=S, pA=pA)

out <- do.call(rbind, lapply(1:nrow(SpA), function(i){
  pars <- c(S = SpA$S[i], pars[-1])
  state <- c(pA = SpA$pA[i])
  flux(state, pars)
})
)

fluxes <- cbind(gather(as.data.frame(out), key = "flux", value = "value", -c(S,pA,dpA)))

ggplot(fluxes, aes(x = pA, y = value, color = flux, group = flux)) + facet_wrap(~as.character(S)) + geom_path() 

• Evaluate the system for smaller values of \(k_2\) and discover another, qualitatively different behavior of the system. How can this behavior be explained?

# change k2 to 0.5
pars["k2"] <- 0.5

outforth <- SSvsS()
outforth$way <- "forth"

yini <- c(pA = tail(outforth,1)$pA)
outback <- SSvsS(S = seq(2, 0, by = -0.01), yini = yini)
outback$way <- "back"

SSvsSplot <- ggplot(rbind(outforth, outback), aes(x = S, y = pA, color = way)) + geom_path()
print(SSvsSplot)

## flux calculations for varying k2 values
flux <- function(state, pars) {
    with(as.list(c(state, pars)), {
      build <- (k1*S+k2*pA^4/(Km^4 + pA^4))*(1-pA)
      decay <- k_1*pA
      dpA <- build - decay
      return(c(dpA = dpA, build = build, decay = decay, k2=k2, pA = pA))
    })
}

k <- seq(0, 2, len=9)
pA <- seq(0, 1, len=100)

kpA <- expand.grid(k=k, pA=pA)

out <- do.call(rbind, lapply(1:nrow(SpA), function(i){
  pars["k2"] <- kpA$k[i]
  state <- c(pA = SpA$pA[i])
  flux(state, pars)
})
)

fluxes <- cbind(gather(as.data.frame(out), key = "flux", value = "value", -c(k2,pA,dpA)))

ggplot(fluxes, aes(x = pA, y = value, color = flux, group = flux)) + facet_wrap(~as.character(k2)) + geom_path() 

Cathedral exercise:

What is the story of the most famous waterspout?