Simulation (2017-07-04)

model.R

################################################################################
# A simulation of the following system:
#
# dm1/dt = -k1 * m1 + f1(p9)   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 = -k1 * m1 + f1(p3)   dp4/dt = mu1 * m1 - nu1 * p4
# dm5/dt = -k2 * m2 + f2(p4)   dp5/dt = mu2 * m2 - nu2 * p5
# dm6/dt = -k3 * m3 + f3(p5)   dp6/dt = mu3 * m3 - nu3 * p6
# dm7/dt = -k1 * m1 + f1(p6)   dp7/dt = mu1 * m1 - nu1 * p7
# dm8/dt = -k2 * m2 + f2(p7)   dp8/dt = mu2 * m2 - nu2 * p8
# dm9/dt = -k3 * m3 + f3(p8)   dp9/dt = mu3 * m3 - nu3 * p9
#
# Parameters: 
#   k1 = 0.91, k2 = 0.92, k3 = 0.93
#   mu1 = 1,   mu2 = 1.1, mu3 = 1.2
#   nu1 = 1.4, nu2 = 1.5, nu3 = 1.6
# Functions:
#  f1(p) = 40 / (1 + p^2)
#  f2(p) = 40 * exp(-3 * p)
#  f3(p) = 40 / (1 + p^4)
#
# Author: Andrey Akinshin
################################################################################

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

f1 <- function(p) 40 / (1 + p^2)
f2 <- function(p) 40 * exp(-3 * p)
f3 <- function(p) 40 / (1 + p^4)

model <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dm1 <- -k1 * m1 + f1(p9)
    dm2 <- -k2 * m2 + f2(p1)
    dm3 <- -k3 * m3 + f3(p2)
    dm4 <- -k1 * m1 + f1(p3)
    dm5 <- -k2 * m2 + f2(p4)
    dm6 <- -k3 * m3 + f3(p5)
    dm7 <- -k1 * m1 + f1(p6)
    dm8 <- -k2 * m2 + f2(p7)
    dm9 <- -k3 * m3 + f3(p8)
    
    dp1 <- mu1 * m1 - nu1 * p1
    dp2 <- mu2 * m2 - nu2 * p2
    dp3 <- mu3 * m3 - nu3 * p3
    dp4 <- mu1 * m1 - nu1 * p4
    dp5 <- mu2 * m2 - nu2 * p5
    dp6 <- mu3 * m3 - nu3 * p6
    dp7 <- mu1 * m1 - nu1 * p7
    dp8 <- mu2 * m2 - nu2 * p8
    dp9 <- mu3 * m3 - nu3 * p9
    
    list(c(dm1, dp1, dm2, dp2, dm3, dp3, dm4, dp4, dm5, dp5, dm6, dp6, dm7, dp7, dm8, dp8, dm9, dp9))
  })
}

parameters <- c(k1 = 0.91, k2 = 0.92, k3 = 0.93,
                mu1 = 1,   mu2 = 1.1, mu3 = 1.2,
                nu1 = 1.4, nu2 = 1.5, nu3 = 1.6)
state1     <- 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, m6 = 0.25, p6 = 0.2, m7 = 36, p7 = 35, m8 = 0.2, p8 = 0.1, m9 = 35, p9 = 36)
state2     <- c(m1 = 36, p1 = 35, m2 = 0.1, p2 = 0.2, m3 = 35, p3 = 36, m4 = 36,  p4 = 35,  m5 = 0.1, p5 = 0.2, m6 = 35,  p6 = 36,  m7 = 36, p7 = 35, m8 = 0.1, p8 = 0.2, m9 = 35, p9 = 36)
times      <- seq(0, 50, by = 0.01)

stateLabel <- "point1" # "point2"
state      <- state1   # state2

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

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

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

dfM9 <- data.frame(
  time = out[,"time"],
  m1 = out[,"m1"],
  m2 = out[,"m2"],
  m3 = out[,"m3"],
  m4 = out[,"m4"],
  m5 = out[,"m5"],
  m6 = out[,"m6"],
  m7 = out[,"m7"],
  m8 = out[,"m8"],
  m9 = out[,"m9"]
)

dfP9 <- data.frame(
  time = out[,"time"],
  p1 = out[,"p1"],
  p2 = out[,"p2"],
  p3 = out[,"p3"],
  p4 = out[,"p4"],
  p5 = out[,"p5"],
  p6 = out[,"p6"],
  p7 = out[,"p7"],
  p8 = out[,"p8"],
  p9 = out[,"p9"]
)


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(paste0(stateLabel,  "-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(paste0(stateLabel, "-plot2d-time-p1-p2-p3.png"))

simPlotM9 <-
  dfM9 %>%
  gather("mRNA", "value", 2:10) %>%
  ggplot(aes(x=time, y=value, group=mRNA, colour=mRNA)) +
  geom_line() +
  theme_bw() +
  ggtitle("Simulation (m1,m2,m3,m4,m5,m6,m7,m8,m9)") +
  xlab("Time (minutes)") +
  ylab("mRNA per cell")
simPlotM9
ggsave(paste0(stateLabel,  "-plot2d-time-m1-m2-m3-m4-m5-m6-m7-m8-m9.png"))

simPlotP9 <-
  dfP9 %>%
  gather("protein", "value", 2:10) %>%
  ggplot(aes(x=time, y=value, group=protein, colour=protein)) +
  geom_line() +
  theme_bw() +
  ggtitle("Simulation (p1,p2,p3,p4,p5,p6,p7,p8,p9)") +
  xlab("Time (minutes)") +
  ylab("Proteins per cell")
simPlotP9
ggsave(paste0(stateLabel, "-plot2d-time-p1-p2-p3-p4-p5-p6-p7-p8-p9.png"))

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

point1-plot2d-time-m1-m2-m3-m4-m5-m6-m7-m8-m9.png

point1-plot2d-time-m1-m2-m3.png

point1-plot2d-time-p1-p2-p3-p4-p5-p6-p7-p8-p9.png

point1-plot2d-time-p1-p2-p3.png

point2-plot2d-time-m1-m2-m3-m4-m5-m6-m7-m8-m9.png

point2-plot2d-time-m1-m2-m3.png

point2-plot2d-time-p1-p2-p3-p4-p5-p6-p7-p8-p9.png

point2-plot2d-time-p1-p2-p3.png