Hi,
I am currently working on parameter inference of ODE models in Stan and was wondering if it is also possible to use delay differential equations in Stan. I know there is no explicit solver for DDE in Stan but is it feasible to implement it manually?

Furthermore, I was wondering if people would consider different Bayesian inference solution if the bulk of one work will be based on time series data and ODE models.

I know that I am in a stan forum but since I am a beginner in the field of Bayesian statistics/inference and I am happy about any advice, I thought I would just ask.

From glancing over wikipedia what DDEs are, I would say that you are very limited in implementing these â€śmanuallyâ€ť. The only case I can think of is whenever you use forcing functions. These forcing functions can be evaluated at any given time point within your ODE functor. However, I am not sure if people would call this a DDE anymore.

Hmmâ€¦ maybe (really maybe) â€¦ you could emulate a DDE by using a nested ODE integration scheme? Should be numerically a nightmare to do and maybe this is non-sense.

If you can formulate it such that the delay is an integer number of steps of the time step of a, say, Euler or Heun discretisation, and you donâ€™t need complicated stuff like event handling, discontinuity propagation (e.g. dde23 in MATLAB), then it works OK. Itâ€™s numerically similar to an AR model with non-zero coefficients with large time lags.

Historically these were the first implementations of DDE solvers until users need fancy things like discontinuity propagation.

If the equations arenâ€™t efficiently integrated by low order scheme, it can be useful to use a nice ODE integrator like Stanâ€™s to do the time stepping, passing the matrix of states over time as a parameter, and reference the state history from the ode func.

The main difficulty in Stan is that interpolation of history isnâ€™t going to work (well).

True. It works and I have done it, but it is a hack and we need better solutions here. Maybe a self-made splines approach could workâ€¦ not sure, but this was a thought of mine recently which I havenâ€™t followed up yet.

I think weâ€™ll have some free cycles over the next year or two for various kinds of pharma-related extensions. As usual, we need to start with some kind of functional spec for what is needed.

Most of the PKPD systems enjoy the simplicity of constant delays, so we can just do the explicit RK + method of steps. Thatâ€™ll be the first thing to do. Non-constant delays can be done in a similar but more messy way.

For state-dependent delays we need to calculate the location of discontinuity on the fly. So the cost is much higher as each discontinuity would incur a solution if an algebraic system. I havenâ€™t thought through the details yet.

If I understand it correctly, one basically solves the ODEs twice and you use the result from the first run to determine the values of the state variables at lagged timepoint when solving the ODEs the second time?
Has anyone tried that before to get an idea how exactly that would look like?

For now, I have luckily found a way to drop the delays and I will be able to apply ODEs to the problem.

I will give @maedoc suggestion a try though since I will certainly be working with DDE in the future.

I meant it in the sense, that if one works mainly with ODEs and potentially SDEs in the future, if you or others would recommend a different environment due to potential limitations of the Stan environment when it comes to dynamical systems?

Regardless, I am pretty happy with Stan and I like how intuitive it is to use and how easy it is for me to further and further expand my model.

real f(real y, real y_delayed);
real dt;
int delay_steps; // delay in units of dt
real y[T + delay_steps];
y[1:delay_steps] = history_func(dt, delay_steps);
for (t in (delay_steps:T+delay_steps))
{
y[t + 1] = y[t] + dt * f(y[t], y[t - delay_steps]);
}

Iâ€™m assuming you need y(t-\tau), not just (t-\tau), as in a time-dependent forcing term, e.g. \frac{dy}{dt}=f(y(t-\tau)) = sin(t-\tau) g(y), since that is normally feasible with ODE solvers.

I think the solution above is the best if your model can be discretized and solved with these simple differences instead of a more sophisticated solver (I often have to implement discrete stochastic versions of ODE models and compare their output to the mean field approximations and I rarely have problems, but of course that is model and parameterization-dependent). If that holds it would be actually what you are asking for, just being a â€śdelay difference equationâ€ť instead.

As an alternative, you may be able to formulate your model with a series of compartments before the component y_n(t-\tau) of the state that should have a delay, something like:

that would delay the â€śavailabilityâ€ť of y_n by replacing an exponential waiting time with a gamma-distributed waiting time resulting form the sum of sequential exponential times, with mean waiting time dependent only on rate \rho and the number of compartments n (although it will increase the number of components of the state \mathbf{y}, in case you need to estimate initial conditions). This would effectively introduce a delay, although Iâ€™m not sure if and how it is related to DDEs.
If that makes sense for your problem it is formulated as an ODE without any time dependence, so it should be straightforward to do it with the Stan solver.

To be clear, I didnâ€™t come up with that myself, some references where gamma-distributed rates affect the dynamics in models that can be written as ODEs can be found in Wallinga & Lipsitch 2006 and Bellan 2010, although Iâ€™m not sure they explicitly formulate that as sequential compartments with exponential waiting times, as I described above, but there should be references in there where that formulation is used.

I came back to this recently and tried to put together a working example for estimating parameters fora stochastic delay differential equation. Hereâ€™s the gist on GitHub, the first variant where the time delay is fixed, the time stepping looks like

for (t in (d + 1):nt)
{
real x1 = x_t_hat[t - 1];
real xd = x_t_hat[t - d];
real dx = x1 - x1*x1*x1/3.0 + I_hat + K_hat * xd + eta_hat[t]*eta_sd;
x_t_hat[t] = x1 + dt * dx;
}

where xd is the delayed state. It looks a lot like a high-lag AR to be honest (though here itâ€™s nonlinear). For the case of trying to estimate the delay, the prior on the delay can be evaluated across the entire time grid and the product with the state over time is the delayed state,

for (t in (nh + 1):nt)
{
real x1 = x_t[t - 1];
real xd = exp(-square((ts - (t*dt - dmu))/dsd))' * x_t; // normal prior on delay
real dx = x1 - x1*x1*x1/3.0 + I + K * xd + eta[t]*eta_sd;
x_t[t] = x1 + dt * dx;
}

The second variant is a lot more intensive than the first, so if any ideas for interpolation have came in the mean time Iâ€™d be really interested. Otherwise I hope this be helpful.

I think you also need adjoints to compute forward sensitivities. It could also make sense to facilitate calling into the Julia differential equations infrastructure where lots of solvers are available and there are plenty of sensitivity options.

@ChrisRackauckas To follow up on the other thread, the workflow which would be straightforward to implement is a user setting up their DEs in Julia, precompiling the solution and gradient calculation to C compatible functions, and generating a header file that can be included in the Stan model. The model would then invoke the Julia solvers during the model run. What do you think?