Drosophila Model

DrosophilaModel.R

library(ggplot2)
library(deSolve)
library(ggpubr)
library(latex2exp)

'%=%' <- function(x, y) {
  x <- as.character(substitute(x)[-1])
  if (length(y) < length(x))
    y <- rep(y, ceiling(length(x) / length(y)))
  if (length(y) > length(x))
    y <- y[1:length(x)]
  mapply(assign, x, y, MoreArgs = list(envir = parent.frame()))
  invisible()
}

h <- function(v) v > 0

dt_x <- 0
dt_y <- 0
dt_z <- 0
dt_u <- 0
dt_w <- 0
dt_p <- 0
Gro <- 1
EMC <- 1
D <- 2.05
L <- 1
UB <- 1.17
SINA <- 5.6
C <- 1
n1 <- 2
n3 <- 2
n5 <- 2
m1 <- 1.77
m2 <- 1
m3 <- 1
m4 <- 1
m5 <- 1
m6 <- 1
m7 <- 1
kx <- 1
ky <- 1
kz <- 1
ku <- 1
kw <- 1
kp <- 1
a3 <- 3.61
b3 <- 4.96
c3 <- 1.35
a4 <- 4.43
b4 <- 6.09
c4 <- 1.66
a5 <- 8.09
b5 <- 11.13
c5 <- 3.03
a6 <- 2.67
b6 <- 3.67
c6 <- 1
d1 <- 7.46
d3 <- 2.77
d5 <- 1.24
x0 <- 0.56
y0 <- 1.59
z0 <- 0.15
u0 <- 0
w0 <- 0
p0 <- 0

sigma1 <- function (v) {
  h(v) * d1 * v ^ n1 / (1 + v ^ n1)
}
sigma3 <- function (v) {
  h(v) * d3 * v ^ n3 / (1 + v ^ n3)
}
sigma5 <- function (v) {
  h(v) * d5 * v ^ n5 / (1 + v ^ n5)
}
s3 <- function (v) {
  h(v) * a3 * exp((v - b3) / c3) / (1 + exp((v - b3) / c3))
}
s4 <- function (v) {
  h(v) * a4 * exp((v - b4) / c4) / (1 + exp((v - b4) / c4))
}
s5 <- function (v) {
  h(v) * a5 * exp((v - b5) / c5) / (1 + exp((v - b5) / c5))
}
s6 <- function (v) {
  h(v) * a6 * exp((v - b6) / c6) / (1 + exp((v - b6) / c6))
}


# Helper for x'
x_helper <- function(t, lagPoint, currentPoint = lagPoint) {
  c(x, y, z, u, w, p) %=% lagPoint
  c(x_curr, y_curr, z_curr, u_curr, w_curr, p_curr) %=% currentPoint
  kx * (sigma1(D * x) + sigma3(z) + sigma5(w)) / ((1 + Gro * y) * (1 + EMC * x)) - m1 * x * (1 + p_curr * UB * SINA)
}

# Modelling for x'
x_model <- function(t, dt, startValue, point, args, ...) {
  if (t < 0) { startValue }
  else if (t < dt || dt < 0.01) {
    lagPoint <- x_helper(t, point)
  } else {
    lagPoint <- lagvalue(t - dt)
    x_helper(t, lagPoint, point)
  }
}

# Helper for y'
y_helper <- function(t, lagPoint, currentPoint = lagPoint) {
  c(x, y, z, u, w, p) %=% lagPoint
  c(x_curr, y_curr, z_curr, u_curr, w_curr, p_curr) %=% currentPoint
  ky * C / (1 + u ^ m7) - m2 * y_curr
}

# Modelling for y'
y_model <- function(t, dt, startValue, point, args, ...) {
  if (t < 0) { startValue }
  else if (t < dt || dt < 0.01) {
    lagPoint <- y_helper(t, point)
  } else {
    lagPoint <- lagvalue(t - dt)
    y_helper(t, lagPoint, point)
  }
}

# Helper for z'
z_helper <- function(t, lagPoint, currentPoint = lagPoint) {
  c(x, y, z, u, w, p) %=% lagPoint
  c(x_curr, y_curr, z_curr, u_curr, w_curr, p_curr) %=% currentPoint
  kz * s3(D * x) - m3 * z_curr
}

# Modelling for z'
z_model <- function(t, dt, startValue, point, args, ...) {
  if (t < 0) { startValue }
  else if (t < dt || dt < 0.01) {
    lagPoint <- z_helper(t, point)
  } else {
    lagPoint <- lagvalue(t - dt)
    z_helper(t, lagPoint, point)
  }
}

# Helper for u'
u_helper <- function(t, lagPoint, currentPoint = lagPoint) {
  c(x, y, z, u, w, p) %=% lagPoint
  c(x_curr, y_curr, z_curr, u_curr, w_curr, p_curr) %=% currentPoint
  ku * s4(D * x) - m4 * u_curr
}

# Modelling for u'
u_model <- function(t, dt, startValue, point, args, ...) {
  if (t < 0) { startValue }
  else if (t < dt || dt < 0.01) {
    lagPoint <- u_helper(t, point)
  } else {
    lagPoint <- lagvalue(t - dt)
    u_helper(t, lagPoint, point)
  }
}

# Helper for w'
w_helper <- function(t, lagPoint, currentPoint = lagPoint) {
  c(x, y, z, u, w, p) %=% lagPoint
  c(x_curr, y_curr, z_curr, u_curr, w_curr, p_curr) %=% currentPoint
  kw * s5(D * x) - m5 * w_curr
}

# Modelling for w'
w_model <- function(t, dt, startValue, point, args, ...) {
  if (t < 0) { startValue }
  else if (t < dt || dt < 0.01) {
    lagPoint <- w_helper(t, point)
  } else {
    lagPoint <- lagvalue(t - dt)
    w_helper(t, lagPoint, point)
  }
}

# Helper for p'
p_helper <- function(t, lagPoint, currentPoint = lagPoint) {
  c(x, y, z, u, w, p) %=% lagPoint
  c(x_curr, y_curr, z_curr, u_curr, w_curr, p_curr) %=% currentPoint
  h(t - 12) * (kp * s6(D * x) * (t - 12) ^ 2 / (L + (t - 12) ^ 2) - m6 * p_curr)
}

# Modelling for p'
p_model <- function(t, dt, startValue, point, args, ...) {
  if (t < 0) { startValue }
  else if (t < dt || dt < 0.01) {
    lagPoint <- p_helper(t, point)
  } else {
    lagPoint <- lagvalue(t - dt)
    p_helper(t, lagPoint, point)
  }
}

# Modelling stuff
simulate <- function() {
  model <- function(t, p, parms, ...) {
    list(c(x_model(t, dt_x, x0, p), y_model(t, dt_y, y0, p), z_model(t, dt_z, z0, p), u_model(t, dt_u, u0, p), w_model(t, dt_w, w0, p), p_model(t, dt_p, p0, p)))
  }
  t <- seq(0, 16, by = 0.1)
  start <- c(x0 * kx, y0 * ky, z0 * kz, u0 * ku, w0 * kw, p0 * kp)
  traj <- as.data.frame(dede(start, t, func = model))
  
  plot_x <- traj[, 2]
  plot_y <- traj[, 3]
  plot_z <- traj[, 4]
  plot_u <- traj[, 5]
  plot_w <- traj[, 6]
  plot_p <- traj[, 7]
  
  df <- rbind(
    data.frame(t, value = plot_x, f = "AS-C: x(t)"),
    data.frame(t, value = plot_y, f = "Hairy: y(t)"),
    data.frame(t, value = plot_z, f = "Senseless: z(t)"),
    data.frame(t, value = plot_u, f = "Scratch: u(t)"),
    data.frame(t, value = plot_w, f = "Charlatan: w(t)"),
    data.frame(t, value = plot_p, f = "Phyllopod: p(t)")
  )
  return(df)
}

draw_plot <- function(df, filename) {
  p <- ggplot(df, aes(x = t, y = value, color = f, group = f)) +
    geom_line(aes(linetype = f)) +
    ylim(0, max(df$value)) +
    ylab("") +
    labs(color = "Function", linetype = "Function") + 
    ggtitle(TeX(paste0("$\\k_x=", kx, "$"))) +
    theme_bw()
  ggsave(paste0(filename, ".png"), p, dpi = 1200, width = 6, height = 4)
  show(p)
}

simulate_and_draw <- function(filename, kx.value) {
  kx <<- kx.value
  df <- simulate()
  draw_plot(df, filename)
}

for (i in 640:660) {
  simulate_and_draw(paste0("figure_k0", i), i * 0.001)
}