Introduction

The purpose of this package is to allow users to simulate and estimate the parameters from a system of ordinary differential equations (ODE) in Stan while not having to write Stan code. Events are supported by providing an event data frame.

Unfortunately, due to the user-defined nature ODE system it is not possible for the Stan model to be precompiled. Thus some compile time will be involve when simulating/fitting the model specified by the user.

NOTE: estimating parameters from the ODE (i.e. sampling = TRUE) is not yet supported.

Libraries

Load the stanode package and the deSolve package (for comparison).

library(rstanode)
library(deSolve)
rstan::rstan_options(auto_write = TRUE)

Simple Harmonic Oscillator

The mathematical representation of a simple harmonic oscillator ODE is, $$ \begin{align} dy_1 &= y_2 \ dy_2 &= -y_1 -\theta \cdot y_2 \end{align} $$

In deSolve we can fit a simple harmonic oscillator as follows,

# specify the ODE
sho <- function(t, y, p) {
  with(as.list(c(y,p)), {
    dy1 = y2
    dy2 = - y1 - theta * y2
    return(list(c(dy1 = dy1, dy2 = dy2)))
  })
}
# declare the parameters, initial conditions, and time steps
pars <- c("theta" = 0.15)
yini <- c("y1" = 1, "y2" = 0)
time_steps <- seq(1, 100, by = 0.001)
sho_deSolve <- ode(func = sho, y = yini, times = time_steps, parms = pars)

In stanode we can use the same function arguments to simulate from the same model,

sho_stanode <- stan_ode(sho, state = yini,
                  pars = pars,
                  times = time_steps,
                  integrator = "rk45",
                  sampling = FALSE)
sims <- sho_stanode$simulations
plot(sho_deSolve[,2], sho_deSolve[,3], col = "#808080", lwd = 3,
     type = "l",
     main = "Simple Harmonic Oscillator", xlab = "y1", ylab = "y2")
lines(sims[,2], sims[,3], col = "#FF6688")
legend("bottomleft", c("deSolve","stanode"), col = c("#808080", "#FF6688"), lwd = c(3, 1))

PKPD Model

The stanode package also supports events. Specifically, the user provides an event table to which describes the time the event takes place, the value that each state should take on at the specfied state, and the method to use (currently "add" and "multiply" are supported).

Below we fit the two compartment PKPD model in deSolve

two_cpt <- function(t, y, parms) {
  with(as.list(c(y, parms)), {
    dy_gut = -ka * y_gut
    dy_cent = ka * y_gut - (CL/V_cent + Q/V_cent) * y_cent + (Q/V_peri) * y_peri
    dy_peri = (Q/V_cent) * y_cent - (Q/V_peri) * y_peri
    return(list(c(dy_gut=dy_gut, dy_cent=dy_cent, dy_peri=dy_peri)))
  })
}
dosing_table <- data.frame(var = "y_gut",
                           time = seq(10, 70, by=10),
                           value = 5,
                           method = "add")
pars <- c("CL" = 10, Q = 13, "V_cent" = 20, "V_peri" = 73, ka = 3)
yini <- c("y_gut" = 0, "y_cent" = 0, "y_peri" = 0)
time_steps <- seq(0, 150, by = 0.005)
two_cpt_deSolve <- deSolve::ode(y = yini, times = time_steps, func = two_cpt, parms = pars,
                                events = list(data = dosing_table))

To fit the model in stanode we need to specify the event table,

dosing_table_stan <- data.frame(time = dosing_table$time,
                                "y_gut" = 5, "y_cent" = 0, "y_peri" = 0,
                                method = "add")
print(dosing_table_stan)

Now we can simulate the model,

two_cpt_stanode <- stan_ode(two_cpt, state = yini,
                            pars = pars,
                            times = time_steps,
                            events = dosing_table_stan,
                            integrator = "rk45",
                            sampling = FALSE)
sims <- two_cpt_stanode$simulations

Below we plot the simulations for each compartment from both packages.

par(mfrow = c(1,3))

plot(two_cpt_deSolve[,1], two_cpt_deSolve[,2], col = "#808080", lwd = 3,
     type = "l",
     main = "Gut Compartment", xlab = "Time", ylab = "Concentration")
lines(sims[,1], sims[,2], col = "#FF6688")
legend("topright", c("deSolve","stanode"), col = c("#808080", "#FF6688"), lwd = c(3, 1))

plot(two_cpt_deSolve[,1], two_cpt_deSolve[,3], col = "#808080", lwd = 3,
     type = "l",
     main = "Central Compartment", xlab = "Time", ylab = "Concentration")
lines(sims[,1], sims[,3], col = "#FF6688")
legend("topright", c("deSolve","stanode"), col = c("#808080", "#FF6688"), lwd = c(3, 1))

plot(two_cpt_deSolve[,1], two_cpt_deSolve[,4], col = "#808080", lwd = 3,
     type = "l",
     main = "Peripheral Compartment", xlab = "Time", ylab = "Concentration")
lines(sims[,1], sims[,4], col = "#FF6688")
legend("topright", c("deSolve","stanode"), col = c("#808080", "#FF6688"), lwd = c(3, 1))

Arenstorf Orbits

Arenstorf <- function(t, y, p) {
  with(as.list(c(y,p)), {
    D1 <- ((y[1] + mu1)^2 + y[2]^2)^(3/2)
    D2 <- ((y[1] - mu2)^2 + y[2]^2)^(3/2)
    dy1 <- y[3]
    dy2 <- y[4]
    dy3 <- y[1] + 2*y[4] - mu2*(y[1]+mu1)/D1 - mu1*(y[1]-mu2)/D2
    dy4 <- y[2] - 2*y[3] - mu2*y[2]/D1 - mu1*y[2]/D2
    return(list( c(dy1, dy2, dy3, dy4) ))
  })
}

Simulate in deSolve.

mu1 <- 0.012277471
pars <- c(mu1 = mu1, mu2 = 1 - mu1)
yini <- c(y1 = 0.994, y2 = 0,
          y3 = 0, y4 = -2.00158510637908252240537862224)
time_steps <- seq(from = 0, to = 18, by = 0.01)
orbit_deSolve <- ode(func = Arenstorf, y = yini, parms = pars, times = time_steps)

Simulate in rstanode.

orbit_rstanode <- stan_ode(Arenstorf,
                           yini,
                           pars,
                           time_steps,
                           integrator = "rk45")
sims <- orbit_rstanode$simulations

Overlay both simulations.

plot(orbit_deSolve[,2], orbit_deSolve[,3], type = "l", col = "#808080", lwd = 3)
lines(sims[,2], sims[,3], col = "#FF6688")
legend("bottomleft", c("deSolve","stanode"), col = c("#808080", "#FF6688"), lwd = c(3, 1))


imadmali/stanode documentation built on May 3, 2019, 11:48 p.m.