Simulation (18-10-16)

model.R

################################################################################
# A simulation of the following system:
#
# dm1/dt = -k1 * m1 + f1(p3)   dp1/dt = mu1 * (m1 - p1)
# dm2/dt = -k2 * m2 + f2(p1)   dp2/dt = mu2 * (m2 - p2)
# dm3/dt = -k3 * m3 + f3(p2)   dp3/dt = mu3 * (m3 - p3)
#
# Parameters: 
#   k1  = 1.5   k2  = 1.7   k3  = 2.0
#   mu1 = 3.0   mu2 = 3.5   mu3 = 4.0
# Functions:
#  f1(p) = 10*e^(-p)  + 0.5
#  f2(p) = 20/(1+p^3) + 0.8
#  f3(p) = 25/(1+p^2) + 0.6
#
# Author: Andrey Akinshin
################################################################################

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

f1 <- function(p) 10 * exp(-p)    + 0.5
f2 <- function(p) 20 / (1 + p^3) + 0.8
f3 <- function(p) 25 / (1 + p^2) + 0.6

model <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dm1 <- -k1 * m1 + f1(p3)
    dm2 <- -k2 * m2 + f2(p1)
    dm3 <- -k3 * m3 + f3(p2)
    dp1 <- mu1 * (m1 - p1)
    dp2 <- mu2 * (m2 - p2)
    dp3 <- mu3 * (m3 - p3)
    list(c(dm1, dm2, dm3, dp1, dp2, dp3))
  })
}

parameters <- c(k1 = 1.5, k2 = 1.7, k3 = 2.0, mu1 = 3.0, mu2 = 3.5, mu3 = 4.0)
state      <- c(m1 = 0, m2 = 0, m3 = 0, p1 = 2, p2 = 2, p3 = 2)
times      <- seq(0, 30, by = 0.01)

out <- ode(y = state, times = times, func = model, parms = parameters)
df <- data.frame(
  time = out[,"time"],
  p1 = out[,"p1"],
  p2 = out[,"p2"],
  p3 = out[,"p3"]
)

simPlot <-
  df %>%
  gather("protein", "value", 2:4) %>%
  ggplot(aes(x=time, y=value, group=protein, colour=protein)) +
  geom_line() +
  theme_bw() +
  ggtitle("Simulation") +
  xlab("Time (minutes)") +
  ylab("Proteins per cell")
simPlot

plot3d(df[,2:4], axes=T, type="l")

plot2d.png

plot3d.png