Use simultaneous ODE in STAN with a ridiculously simple example

Here is another example for simultaneous ODE is involved in the model. In this example, we assume the model with 1st order absorption. Carefully check the difference between this one and the previous example. You could get an overall view of what are necessary to use ODE in STAN.

1-compartment model for 1st order absorption

Consider the following 1-compartment model for 1st order absorption with log normal error. Let\(\hat{C}(t)\) be the concentration at time \(t\), \(\displaystyle k\left(:=\frac{CL}{V}\right)\) be the elimination constant, \(G(t)\) be the amount of drug in the absorption compartment, and \(k_a\) be the 1st order absorption constant, then we have the following simultaneous ODE: \begin{align} \left\{ \begin{array}{l} \displaystyle \frac{\partial G(t)}{\partial t}=-k_aG(t) \\ \displaystyle \frac{\partial \hat{C}(t)}{\partial t}=k_aG(t)-k\hat{C}(t) \end{array} \right. \label{eq:odeabs} \end{align}

Data generation

Analytical solution for Eq.(\ref{eq:odeabs}) is: \begin{align} \hat{C}(t)=\frac{D}{V}\frac{k_a}{k_a-k} \left( \exp \left( -kt\right) - \exp \left( -k_at\right) \right), \label{eq:anasolabs} \end{align} where \(D\) is amount of dose (bioavailability is not considered here.) and \(V\) is volume of distribution. We can generate simulation dataset using Eq.(\ref{eq:anasolabs}) as follows: \begin{align} C(t)\sim \mbox{lognormal} \left( \log(\hat{C}(t)), \sigma \right), \label{eq:geneabs} \end{align} where \(\sigma\) is standard deviation of the measurement error. Suppose we measured from \(t=1.0\) to \(24.0\) by \(0.5\) hour and \(D=20\), \(V=30\), \(CL=6\), \(k_a=0.9\) and \(\sigma=0.05\).

library(rstan)

D = 20
theta=c(30,6,0.9)
t = seq(1,24,0.5)
sigma = 0.05

FOSD = function(theta,t){
  V  = theta[1]
  k  = theta[2]/theta[1]
  ka = theta[3]
  
  res = D/V*ka/(ka-k)*(exp(-k*t)-exp(-ka*t))
  return(res)
}

SCFOSD = function(theta,t,sigma){
  yhat = FOSD(theta,t)
  res = yhat*exp(rnorm(length(yhat),0,sigma))
  return(res)
}

res = c()
for (i in 1:length(t)){
  res = c(res,SCFOSD(theta,t[i],sigma))
}

p = ggplot(data.frame(t, C=res))
p + geom_point(aes(x = t, y = C))

The generated dataset is shown as below.

Implementation in STAN

The following STAN code is using the function "integrate_ode_rk45". The functions \(G(t)\) and \(\hat{C}(t)\) are assigned to \(y[1]\) and \(y[2]\) respectively. However, only \(y[2]\) is used in the model section. The parameter K in transformed parameters section is set in order to avoid flip-flop kinetics, i.e., \(k_a>k\).

model_code_1_comp_abs='
functions {
real[] absorp(real t,
              real[] y,
              real[] theta,
              real[] x_r,
              int[] x_i) {
  real dydt[2];
  real k;
  real ka;
  k = theta[2]/theta[1];
  ka = theta[3];
  dydt[1] = -ka*y[1];
  dydt[2] = ka*y[1]-k*y[2];
  return dydt;
}

}
data {
  int T;
  real t0;
  real ts[T];
  real y[T];
}

transformed data {
  real x_r[0];
  int x_i[0];
}

parameters{
  real sigma;
  real theta[3]; // V, CL, Ka
}

transformed parameters{
  positive_ordered[2] K;
  real y0[2];
  real y_hat[T, 2];

  K[1]=theta[2]/theta[1];
  K[2]=theta[3];
  
  y0[1] = 20/theta[1];
  y0[2] = 0;

  y_hat = integrate_ode_rk45(absorp, y0, t0, ts, theta, x_r, x_i);
}

model {
  for (i in 1:T)
    y[i] ~ lognormal(log(y_hat[i,2]),sigma);
}
'

model_1_comp_abs <- stan_model(model_code = model_code_1_comp_abs)

Then, we can do sampling by the following script. The model mathematically has 2 solutions (ill-identified condition) and "init_fun" is necessary to converged to correct solution and make it faster.

init_fun <- function(){list(theta=c(runif(1,10,50),runif(1,2,10),runif(1,0.5,2)))}

fit <- sampling(model_1_comp_abs,
                 data = list(T = length(t),
                             ts = t,
                             t0 = 0,
                             y = res),
                 chains = 1, iter = 1000, init=init_fun
)

We get the following estimates and all the parameters are almost accurately recovered.

> print(fit)
Inference for Stan model: c9107cb9548c51f0c9f9926746ed55b0.
1 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=500.

              mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
sigma         0.05    0.00 0.01   0.04   0.05   0.05   0.06   0.06   219    1
theta[1]     30.64    0.03 0.45  29.77  30.32  30.64  30.94  31.46   225    1
theta[2]      6.09    0.00 0.06   5.97   6.04   6.09   6.13   6.23   250    1
theta[3]      0.88    0.00 0.04   0.80   0.85   0.88   0.91   0.98   302    1
K[1]          0.20    0.00 0.00   0.20   0.20   0.20   0.20   0.20   309    1
K[2]          0.88    0.00 0.04   0.80   0.85   0.88   0.91   0.98   302    1
y0[1]         0.65    0.00 0.01   0.64   0.65   0.65   0.66   0.67   227    1
y0[2]         0.00    0.00 0.00   0.00   0.00   0.00   0.00   0.00   500  NaN
y_hat[1,1]    0.27    0.00 0.01   0.24   0.26   0.27   0.28   0.29   296    1
y_hat[1,2]    0.34    0.00 0.01   0.32   0.33   0.34   0.35   0.36   302    1
y_hat[2,1]    0.17    0.00 0.01   0.15   0.17   0.17   0.18   0.20   304    1
y_hat[2,2]    0.40    0.00 0.01   0.38   0.39   0.40   0.41   0.42   292    1
y_hat[3,1]    0.11    0.00 0.01   0.09   0.11   0.11   0.12   0.13   308    1
y_hat[3,2]    0.42    0.00 0.01   0.41   0.42   0.42   0.43   0.44   277    1
                                 (snip)
y_hat[45,1]   0.00    0.00 0.00   0.00   0.00   0.00   0.00   0.00   254    1
y_hat[45,2]   0.01    0.00 0.00   0.01   0.01   0.01   0.01   0.01   500    1
y_hat[46,1]   0.00    0.00 0.00   0.00   0.00   0.00   0.00   0.00   283    1
y_hat[46,2]   0.01    0.00 0.00   0.01   0.01   0.01   0.01   0.01   491    1
y_hat[47,1]   0.00    0.00 0.00   0.00   0.00   0.00   0.00   0.00   298    1
y_hat[47,2]   0.01    0.00 0.00   0.01   0.01   0.01   0.01   0.01   478    1
lp__        118.65    0.13 1.48 115.07 117.85 118.99 119.76 120.51   135    1

Not to mentioned, you can include individual variability with this ODE model with simple modifications.