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\).

## 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()

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)

## 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() 

# 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?