Simulation (2017-06-15)

model.R

################################################################################
# A simulation of the following system:
#
# dm1/dt = -k1 * m1 + f1(p5)   dp1/dt = mu1 * m1 - nu1 * p1
# dm2/dt = -k2 * m2 + f2(p1)   dp2/dt = mu2 * m2 - nu2 * p2
# dm3/dt = -k3 * m3 + f3(p2)   dp3/dt = mu3 * m3 - nu3 * p3
# dm4/dt = -k4 * m4 + f4(p3)   dp4/dt = mu4 * m4 - nu4 * p4
# dm5/dt = -k5 * m5 + f5(p4)   dp5/dt = mu5 * m5 - nu5 * p5
#
# Parameters: 
#   k1 = 0.91, k2 = 0.92, k3 = 0.93, k4 = 0.94, k5 = 0.95
#   mu1 = 1,   mu2 = 1.1, mu3 = 1.2, mu4 = 1.3, mu5 = 1.4
#   nu1 = 1.4, nu2 = 1.5, nu3 = 1.6, nu4 = 1.7, nu5 = 1.8
# Functions:
#  f1(p5) = 40 / (1 + p5^2)
#  f2(p1) = 40 * exp(-2 * p1)
#  f3(p2) = 40 / (1 + p2^3)
#  f4(p3) = 40 * exp(-3 * p3)
#  f5(p4) = 40 / (1 + p4^4)
#
# Author: Andrey Akinshin
################################################################################

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

f1 <- function(p5) 40 / (1 + p5^2)
f2 <- function(p1) 40 * exp(-2 * p1)
f3 <- function(p2) 40 / (1 + p2^3)
f4 <- function(p3) 40 * exp(-3 * p3)
f5 <- function(p4) 40 / (1 + p4^4)

model <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dm1 <- -k1 * m1 + f1(p5)
    dm2 <- -k2 * m2 + f2(p1)
    dm3 <- -k3 * m3 + f3(p2)
    dm4 <- -k4 * m4 + f4(p3)
    dm5 <- -k5 * m5 + f5(p4)
    dp1 <- mu1 * m1 - nu1 * p1
    dp2 <- mu2 * m2 - nu2 * p2
    dp3 <- mu3 * m3 - nu3 * p3
    dp4 <- mu4 * m4 - nu4 * p4
    dp5 <- mu5 * m5 - nu5 * p5
    list(c(dm1, dp1, dm2, dp2, dm3, dp3, dm4, dp4, dm5, dp5))
  })
}

parameters <- c(k1 = 0.91, k2 = 0.92, k3 = 0.93, k4 = 0.94, k5 = 0.95,
                mu1 = 1,   mu2 = 1.1, mu3 = 1.2, mu4 = 1.3, mu5 = 1.4,
                nu1 = 1.4, nu2 = 1.5, nu3 = 1.6, nu4 = 1.7, nu5 = 1.8)
state      <- c(m1 = 36, p1 = 35, m2 = 0.1, p2 = 0.2, m3 = 35, p3 = 36, m4 = 0.2, p4 = 0.1, m5 = 34, p5 = 33)
times      <- seq(0, 50, by = 0.01)

out <- ode(y = state, times = times, func = model, parms = parameters)

dfM <- data.frame(
  time = out[,"time"],
  m1 = out[,"m1"],
  m2 = out[,"m2"],
  m3 = out[,"m3"]
)

dfP <- data.frame(
  time = out[,"time"],
  p1 = out[,"p1"],
  p2 = out[,"p2"],
  p3 = out[,"p3"]
)

simPlotM <-
  dfM %>%
  gather("mRNA", "value", 2:4) %>%
  ggplot(aes(x=time, y=value, group=mRNA, colour=mRNA)) +
  geom_line() +
  theme_bw() +
  ggtitle("Simulation (m1,m2,m3)") +
  xlab("Time (minutes)") +
  ylab("mRNA per cell")
simPlotM
ggsave("plot2d-time-m1-m2-m3.png")

simPlotP <-
  dfP %>%
  gather("protein", "value", 2:4) %>%
  ggplot(aes(x=time, y=value, group=protein, colour=protein)) +
  geom_line() +
  theme_bw() +
  ggtitle("Simulation (p1,p2,p3)") +
  xlab("Time (minutes)") +
  ylab("Proteins per cell")
simPlotP
ggsave("plot2d-time-p1-p2-p3.png")

plot3d(dfM[,2:4], axes=T, type="l", col="red")
plot3d(dfP[,2:4], axes=T, type="l", col="red")

plot2d-time-m1-m2-m3.png

plot2d-time-p1-p2-p3.png

plot3d-m1-m2-m3.png

plot3d-p1-p2-p3.png