String together multiple ODE initializations simultaneously?

Here’s a possibly odd question: does anyone know how to supply more than a single set of initial conditions – both for time t0 and the state of the system R0 – for the integrate_ode_rk45() function?

Context: I’m fitting ODEs to empirically observed experimental data, where we model each treatment as a separate ODE; each ODE is assumed to have the same parameter values *, which we want to estimate. We know initial conditions for both state and time, and we observed the system change through time. Easy enough so far.

Here’s where it gets difficult: we (intentionally) “re-initialized” our experimental conditions twice after first beginning the experiment, but we did this without shuffling treatments. In essence, we have a sequence of three temporally correlated experimental trials. The key for why we want to “re-initialize” the ODEs in Stan, and why I put the * asterisk above: we clearly have time-dependency in the parameters we’re wanting to estimate (i.e., the parameter values change across the three experimental trials), and this time-dependency is especially of interest for the biology we’re trying to model.

Here is simplified R code that uses an ODE to simulate data, add noise, and then set up a list to pass to Stan. This first example does NOT include multiple ODE initializations, and this first example works well:

## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## ~~~~~~~~~~~ Simulate an ODE to fit in STAN ~~~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~





## open ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
rm(list=ls())
library(rstan)
library(StanHeaders)
library(dplyr)
library(deSolve)
library(ggplot2)
library(reshape)
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~





## Simulate ODE ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## params 
params <- c(k=.025)

## create a list of initial conditions
R0 <- seq(50, 300, length.out=5)

## specify ODE function 
resourceLoss <- function(t, R0, params) {
  with(as.list(c(R0, params)),{
    dR <- R0 * -k
    list(c(dR))
  })
}

## set vector of times 
t.list <- seq(1, 50, length.out=3)


#t.list <- c(0, 1, 100)

## run ODE 
out <- ode(y=c(R0), times=t.list, func=resourceLoss, parms=params)

## save as data frame
loss <- round(as.data.frame(out),3)

## save list without Times for later use 
dRdt.list <- round(out[,-1],3)

## create and add prefix for col names
loss.dat <- as.data.frame(dRdt.list)
colnames(loss.dat) <- paste("T", colnames(loss.dat), sep="")

## plot all dRdt 
all_dRdt <- melt(loss, id.vars="time")

tplot <- ggplot(all_dRdt, aes(time, value, col=variable,group=variable)) +
  geom_point() + geom_line(color="black") +
  xlab("Time") + ylab("Density") + ggtitle("Remaining density = R0 - dRdT")

print(tplot)
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~





## add noise ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## set dims for matrix and set sigma 
rws=length(t.list)  ## number of rows ie data points to simulate
cls=ncol(dRdt.list)  ## number of columns ie separate FR runs 
sd=5

set.seed(10)

## apply rgamma to deterministic data 
mult_ODE <- function(n, dat, noise){
  noise.out <- rnorm(n=length(dRdt.list), mean=dat, sd=sd)
  return(matrix(noise.out, nrow=rws, ncol=cls, byrow=FALSE))
}

noisy.dat <- round(mult_ODE(n, dat=dRdt.list, noise=sig),digits=3)

## replace noisy row 1 with deterministic R0 
noisy.dat[1,] <- dRdt.list[1,]

## create data frame for plotting, add time as a column 
noise.dat.df <- as.data.frame(noisy.dat)
noise.dat.df <- cbind(loss$time, noise.dat.df)
names(noise.dat.df)[1]<-"time"

## plot 
melt_noise <- melt(noise.dat.df, id.vars="time")

noise_plot <- ggplot(melt_noise, aes(time, value, col=variable, group=variable)) +
  geom_point() + geom_line(data=all_dRdt, aes(time, value, group=variable), col="black") +
  xlab("Time") + ylab("Density") + ggtitle("Density remaining = R0 - dRdT + noise")

print(noise_plot)
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~





## configure data list for STAN ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## extract deterministic list of initial conditions 
R0.list <- as.numeric(loss.dat[1,])


## list for Stan
## standard time
resourceLoss.dat <- list(Tm = nrow(loss.dat)-1,
                 t0 = t.list[1],
                 th = t.list[-1],
                 nLev = ncol(noisy.dat),
                 R0_obs = R0.list,
                 R_obs = noisy.dat[-1,])


## run STAN
fit1 <- stan(
  file = "ExpDecay_StanExample.stan",
  data = resourceLoss.dat, 
  chains = 1, 
  warmup = 1000,
  iter = 2000,
  cores = 1,
  refresh = 10)
## END standard time ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

And here’s the associated code in Stan, and again, this too works fine:

functions {
   real[] resourceLoss(
       real t,              // time
       real[] R,            // state 
       real[] theta,        // params
       real[] x_r, 
       int[] x_i) {
		real dRdt[x_i[1]];
		real k; 
		k = theta[1]; 
				
		for (i in 1:x_i[1]) { 
		    dRdt[i] = R[i] * -k;
				     }    
                return dRdt;
	 	    }
}
data {
   int <lower=1> Tm;                // Tm = length(t.list - initial t point)
   int <lower=1> nLev;              // number of different R0, i.e. # of ODEs
   real R0_obs[nLev];               // Intial conditions R0 list 
   real <lower=1> t0;               // starting time t = 1 
   real <lower=t0> th[Tm];          // array of times for ODE, length Tm  
   real R_obs[Tm, nLev];            // simulated dat 
}
transformed data{
  real x_r[0];                
  int x_i[1];
  x_i[1] = nLev;      
}
parameters {
   real <lower=0> k;        
   real <lower=0> sigma;      
}
transformed parameters {
   real theta[1];
   theta[1] = k;
}
model { 
   real R[Tm, nLev]; 
 
   R = integrate_ode_rk45(resourceLoss, R0_obs, t0, th, theta, x_r, x_i);

   k ~ normal(0,5);            
   sigma ~ normal(0,5);

   for (t in 1:Tm) {
     for (i in 1:nLev) {
       R_obs[t,i] ~ normal(R[t,i], sigma);
                       }
                   }  
}

Here’s where I attempt to pass a vector for t0. Note that this code does not even include the “multiple re-initializations” of the ODE . . . I’m first simply trying to get Stan to accept a vector for t0, after which I’d then expand to include multiple re-initializations (this code replaces the end of the code above, starting at ## configure data list for STAN ~~~)

## treating intial time as a vector
t.list <- seq(1, 50, length.out=3)
t0.list <- rep(t.list[1], ncol(loss.dat))
t.obs_1 <- rep(t.list[2], ncol(loss.dat))
t.obs_2 <- rep(t.list[3], ncol(loss.dat))                   

## all time points -- partitioned and passed to Stan below 
full.t.list <- rbind(t0.list, t.obs_1, t.obs_2)

## extract deterministic list of initial conditions 
R0.list <- as.numeric(loss.dat[1,])

## intialize list for Stan 
resourceLoss.dat <- list(Tm = nrow(loss.dat)-1,
                         t0 = full.t.list[1,],
                         th = full.t.list[-1,],
                         nLev = ncol(loss.dat),
                         R0_obs = R0.list,
                         R_obs = array(loss.dat[-1,]))


## run STAN
fit1 <- stan(
  file = "ExpDecay_StanExample.stan",
  data = resourceLoss.dat, 
  chains = 1, 
  warmup = 1000,
  iter = 2000,
  cores = 1,
  refresh = 10)

Unsurprisingly, I ran into:

Wrong signature for function integrate_ode_rk45; first argument must be the name of a function with signature (real, real[ ], real[ ], data real[ ], data int[ ]): real[ ].

In case it’s at all useful, here’s the conceptual layout of what the initial state + observed values would look like, Where R0_X are the experimentally supplied initial conditions for the state of the system, and DX represent observed data points as the state of the system changes through time; we supply all these values; this represents four treatments and thus 4 ODEs “re-initialized” twice, i.e., three total temporal ODE “runs”:

R0_1, R0_1, R0_1, R0_1
D2, D2, D2, D2,
D3, D3, D3, D3,
R0_2, R0_2, R0_2, R0_2,
D4, D4, D4, D4,
D5, D5, D5, D5,
R0_3, R0_3, R0_3, R0_3
D6, D6, D6, D6,
D7, D7, D7, D7,

And where the first initializing time point is t0_1 = 1, the second initialization of the system is t0_2 = 4, and the third is t0_3 = 7. (what I haven’t included here is a vector of data with “cumulative time”, and it is this vector that will be the predictor for the parameter’s temporal dependency).

t0_1, t0_1, t0_1, t0_1,
2, 2, 2, 2,
3, 3, 3, 3,
t0_2, t0_2, t0_2, t0_2,
5, 5, 5, 5,
6, 6, 6, 6,
t0_3, t0_3, t0_3, t0_3,
8, 8, 8, 8,
9, 9, 9, 9

In essence, we want to string together three separate initializations of an ODE system. And it seems one way to do this is to pass a vector to t0 in Stan. Is this at all possible? If it’s not possible in the context of integrate_ode_rk45(), is there some other context that might work?

Thank you all in advance!!

I don’t think that something like this is possible.

Can you not just do it serially? Or even in parallel, things don’t depend on each other as I understood it.

1 Like

In case anyone else comes across this problem (and is unfamiliar with more advanced ODE principals, as I am), I believe what I attempted to describe in my post here is an example of forced system, for which there are excellent resources and links discussed here.

I’ve since reformulated my model such that this forced approach is not necessary, and instead, we model the effects of time by capturing the change in a parameter when a single ODE system is simultaneously fit to all the data, and where binary indicator variables enable different versions of the parameter to initialize for the different time periods of data. In short, @Funko_Unko, I am indeed following your suggestion of running in parallel :-)

Thanks, folks – the resources accessible on this forum are really helpful!!

Like this? Forcing function (differential equations) - Wikipedia

Just with the forcing function being a series of Dirac deltas?

I think the forcing function page you linked to is an example of what’s utilized to analyze this type of problem.

Here’s another page with more info, from StanCon 2019. Example two (starting on page 21) talks about the pharmacokinetics (PK) and pharmacodynamics (PD) problem that’s addressed with a two compartment model. This model also uses an “event schedule”, where an exterior force alters the state of the system. In the PK/PD problem, this simulates reapplying a dose of a drug at certain time points (check out page 27, where the jumps in the state of the system are the “forced” applications of a treatment, with the subsequent dynamics being the absorption, metabolism, etc). This event schedule appears to use a Torsten function.

Exerting this external force upon system state (via the event schedule) is precisely what I was originally attempting to describe. Full disclosure, I’m officially beyond the edge of my competency here, so take my interpretation with a grain of salt. As mentioned with my previous post, I’ve backed away from this type of modeling approach, and I’m now pursuing a simpler solution that is still process orientated.

Without looking into the code, let’s get straight with the terms, as I’m not sure what “initialization” means here.

Are you trying to do the following?

Solve for u_i(t;\theta), i=1,2,3,4 as treatment index, in t_0<t<t_{\text{end}}, s.t.

u_i' = f_i(u_i, t; \theta(t)), \quad u_i(t_0) = u_i^0, t_0\le t < t_1\\ u_i' = f_i(u_i, t; \theta(t)), \quad u_i(t_1) = u_i^1, t_1\le t < t_2\\ u_i' = f_i(u_i, t; \theta(t)), \quad u_i(t_2) = u_i^2, t_2\le t < t_{\text{end}}
2 Likes