Is it possible to solve delay differential equation?

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.

Thanks for your help and advice.

If you can reduce it to an ODE with unknown initial conditions and parameters, then yes. I’m not familiar with DDEs myself—@wds15 should know.

Not sure what you mean by different Bayesian inference solution. Different than MCMC?

Sorry, I am not familiar with these either.

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.

1 Like

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.

This is on my to-do list.

What specifically?

One use case from my angle are forcing functions such as dosing inputs to PK/PD systems or PK functions as inputs to PD models, etc.

I guess you plan something general rather than a specific pharma thing and I am curious what your plans are.

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.

Not twice, I’d hope. I meant something like

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 had dropped the ball on this earlier this year but I want to pick it up again.

@maedoc could explain would you mean by history_func? How would I integrate that into an ODE function which I input into the ODE solver?

The history func here is just to provide the state of the system for t < 0, it could be just zero to keep things simple

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:

\frac{dy_1}{dt} = f_1(\mathbf{y}) \\ \frac{dy_2}{dt} = \rho (y_1 - y_2) \\ \frac{dy_3}{dt} = \rho (y_2 - y_3) \\ ... \\ \frac{dy_n}{dt} = \rho y_{n-1} + f(\mathbf{y}) \\

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.

Thanks both. I will have a go at this and report back if I get stuck or something is unclear

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 have a C++ binding for dde_solver. But to make it work with Stan we need

  • Add fortran compiler.
  • Add complex-number support.

The complex support is at the corner, “which is nice”.


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?