Simulation (2018-03-09)

model.R

################################################################################
# A simulation of the following system:
#
# dx1/dt =          - k1 * x1 + f1(x9)
# dx2/dt = mu2 * x1 - k2 * x2
# dx3/dt = mu3 * x2 - k3 * x3
# dx4/dt = mu4 * x3 - k4 * x4
# dx5/dt =          - k5 * x5 + f5(x4)
# dx6/dt = mu6 * x5 - k6 * x6
# dx7/dt = mu7 * x6 - k7 * x7
# dx8/dt =          - k8 * x8 + f8(x7)
# dx9/dt = mu9 * x8 - k9 * x9
#
# Parameters (1):
#   k1 = 1.1, k2 = 1.2, k3 = 1.3, k4 = 1.4, k5 = 1.5, k6 = 1.7, k7 = 1.7, k8 = 1.8, k9 = 1.9
#   mu2 = 1.15, mu3 = 1.25, mu4 = 1.35, mu6 = 1.55, mu7 = 1.65, mu9 = 1.85
# Parameters (2):
#   k1 = 0.8, k2 = 0.8, k3 = 0.8, k4 = 0.8, k5 = 0.6, k6 = 0.6, k7 = 0.6, k8 = 0.6, k9 = 0.6
#   mu2 = 0.8, mu3 = 0.8, mu4 = 0.8, mu6 = 0.6, mu7 = 0.6, mu9 = 0.6
# Functions:
#  f1(p) = 30 * exp(-3 * p)
#  f5(p) = 40 / (1 + p^2)
#  f8(p) = 30 * exp(-4 * p)
#
# Author: Andrey Akinshin
################################################################################

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

f1 <- function(p) 30 * exp(-3 * p)
f5 <- function(p) 40 / (1 + p^2)
f8 <- function(p) 30 * exp(-4 * p)

model <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dx1 <-          - k1 * x1 + f1(x9)
    dx2 <- mu2 * x1 - k2 * x2
    dx3 <- mu3 * x2 - k3 * x3
    dx4 <- mu4 * x3 - k4 * x4
    dx5 <-          - k5 * x5 + f5(x4)
    dx6 <- mu6 * x5 - k6 * x6
    dx7 <- mu7 * x6 - k7 * x7
    dx8 <-          - k8 * x8 + f8(x7)
    dx9 <- mu9 * x8 - k9 * x9
    
    list(c(dx1, dx2, dx3, dx4, dx5, dx6, dx7, dx8, dx9))
  })
}

parameters1 <- c(k1 = 1.1, k2 = 1.2, k3 = 1.3, k4 = 1.4, k5 = 1.5, k6 = 1.7, k7 = 1.7, k8 = 1.8, k9 = 1.9,
                 mu2 = 1.15, mu3 = 1.25, mu4 = 1.35, mu6 = 1.55, mu7 = 1.65, mu9 = 1.85)
parameters2 <- c(k1 = 0.8, k2 = 0.8, k3 = 0.8, k4 = 0.8, k5 = 0.6, k6 = 0.6, k7 = 0.6, k8 = 0.6, k9 = 0.6,
                 mu2 = 0.8, mu3 = 0.8, mu4 = 0.8, mu6 = 0.6, mu7 = 0.6, mu9 = 0.6)
state       <- c(x1 = 0, x2 = 0, x3 = 0, x4 = 0, x5 = 0, x6 = 0, x7 = 0, x8 = 0, x9 = 0)
times1      <- seq(0, 40, by = 0.01)
times2      <- seq(0, 100, by = 0.01)

out1 <- ode(y = state, times = times1, func = model, parms = parameters1)
out2 <- ode(y = state, times = times2, func = model, parms = parameters2)

traj1 <- data.frame(
  time = out1[,"time"],
  x1 = out1[,"x1"],
  x2 = out1[,"x2"],
  x3 = out1[,"x3"],
  x4 = out1[,"x4"],
  x5 = out1[,"x5"],
  x6 = out1[,"x6"],
  x7 = out1[,"x7"],
  x8 = out1[,"x8"],
  x9 = out1[,"x9"]
)
traj2 <- data.frame(
  time = out2[,"time"],
  x1 = out2[,"x1"],
  x2 = out2[,"x2"],
  x3 = out2[,"x3"],
  x4 = out2[,"x4"],
  x5 = out2[,"x5"],
  x6 = out2[,"x6"],
  x7 = out2[,"x7"],
  x8 = out2[,"x8"],
  x9 = out2[,"x9"]
)

simPlot1 <-
  traj1 %>%
  gather("var", "value", 2:10) %>%
  ggplot(aes(x=time, y=value, group=var, colour=var)) +
  geom_line() +
  theme_bw() +
  ggtitle("Plot ov vector-function x_i[t]") +
  xlab("Time") +
  ylab(expression(x[i]))
simPlot1
ggsave("system1-plot2d.png")

simPlot2 <-
  traj2 %>%
  gather("var", "value", 2:10) %>%
  ggplot(aes(x=time, y=value, group=var, colour=var)) +
  geom_line() +
  theme_bw() +
  ggtitle("Plot ov vector-function x_i[t]") +
  xlab("Time") +
  ylab(expression(x[i]))
simPlot2
ggsave("system2-plot2d.png")

plot3d(traj1[,c(10,3,5)], axes=T, type="l", col="red")
plot3d(traj1[,c(9,2,6)], axes=T, type="l", col="red")

plot3d(traj2[,c(10,3,5)], axes=T, type="l", col="red")
plot3d(traj2[,c(9,2,6)], axes=T, type="l", col="red")

System 1

System 2