Simulation (2018-12-09)

model.R

################################################################################
# A simulation of the following system:
#
# dm1/dt = L1(p3) - k1 * m1
# dm2/dt = L2(p1) - k2 * m2
# dm3/dt = L3(p2) - k3 * m3
# dp1/dt = G1(m1) - l1 * p1
# dp2/dt = G2(m2) - l2 * p2
# dp3/dt = G3(m3) - l3 * p3
#
# Parameters:
#   k1 = 1.2   a1 = 1.3  l1 = 1.4   b1 = 1.5
#   k2 = 2.1   a2 = 2.3  l2 = 2.4   b2 = 2.5  
#   k3 = 3.1   a3 = 3.3  l3 = 3.4   b3 = 3.5 
# Functions:
#  L1(p) = | k1 * a1 for p <= 1
#          | 0       for p >  1
#  L2(p) = | k2 * a2 for p <= 1
#          | 0       for p >  1
#  L3(p) = | k3 * a3 for p <= 1
#          | 0       for p >  1
#  G1(m) = | l1 * b1 for m >  1
#          | 0       for m <= 1
#  G2(m) = | l2 * b2 for m >  1
#          | 0       for m <= 1
#  G3(m) = | l3 * b3 for m >  1
#          | 0       for m <= 1
#
# Author: Andrey Akinshin
################################################################################

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

L1 <- function(p, k1, a1) ifelse(p <= 1, k1 * a1, 0)
L2 <- function(p, k2, a2) ifelse(p <= 1, k2 * a2, 0)
L3 <- function(p, k3, a3) ifelse(p <= 1, k3 * a3, 0)
G1 <- function(m, l1, b1) ifelse(m >  1, l1 * b1, 0)
G2 <- function(m, l2, b2) ifelse(m >  1, l2 * b2, 0)
G3 <- function(m, l3, b3) ifelse(m >  1, l3 * b3, 0)

model <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {

    dm1 <- L1(p3, k1, a1) - k1 * m1
    dm2 <- L2(p1, k2, a2) - k2 * m2
    dm3 <- L3(p2, k3, a3) - k3 * m3
    dp1 <- G1(m1, l1, b1) - l1 * p1
    dp2 <- G2(m2, l2, b2) - l2 * p2
    dp3 <- G3(m3, l3, b3) - l3 * p3
    list(c(dm1, dm2, dm3, dp1, dp2, dp3))
  })
}

parameters <- c(k1 = 1.2, a1 = 1.3, l1 = 1.4, b1 = 1.5,
                k2 = 2.1, a2 = 2.3, l2 = 2.4, b2 = 2.5,  
                k3 = 3.1, a3 = 3.3, l3 = 3.4, b3 = 3.5)
state      <- c(m1 = 0, m2 = 0, m3 = 0, p1 = 0, p2 = 0, p3 = 0)
times      <- seq(0, 40, by = 0.01)

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

traj.m <- data.frame(
  time = out[,"time"],
  m1   = out[,"m1"],
  m2   = out[,"m2"],
  m3   = out[,"m3"]
)
traj.p <- data.frame(
  time = out[,"time"],
  p1   = out[,"p1"],
  p2   = out[,"p2"],
  p3   = out[,"p3"]
)

simPlot.m <-
  traj.m %>%
  gather("var", "value", 2:4) %>%
  ggplot(aes(x=time, y=value, group=var, colour=var)) +
  geom_line() +
  theme_bw() +
  ggtitle("Plot of vector-function m_i[t]") +
  xlab("Time") +
  ylab(expression(m[i]))
simPlot.m
ggsave("plot.m.2d.png")

simPlot.p <-
  traj.p %>%
  gather("var", "value", 2:4) %>%
  ggplot(aes(x=time, y=value, group=var, colour=var)) +
  geom_line() +
  theme_bw() +
  ggtitle("Plot of vector-function p_i[t]") +
  xlab("Time") +
  ylab(expression(p[i]))
simPlot.p
ggsave("plot.p.2d.png")

plot3d(traj.m[,c(2,3,4)], axes=T, type="l", col="red")
plot3d(traj.p[,c(2,3,4)], axes=T, type="l", col="red")