Distributed Delay Equation conversion to ODEs

[ChrisRackauckas]

https://link.springer.com/article/10.1007/s10928-018-9570-4

[korsbo]

I have been thinking about this in the context of DiffEqBiological.jl but I quickly put that on the shelf since it seemed impractical. How should this be implemented here?

A naive way of achieving this result would be something like:

using ModelingToolkit
using DifferentialEquations
using Plots

function distributed_delay!(eqs, input, output, n, r)
    internals = [Variable(:internal, i)(t) for i in 1:n-1]
    vars = vcat(input, internals) 
    append!(eqs, [D(vars[i]) ~ r * (vars[i-1] - vars[i]) for i in 2:n])
    push!(eqs, D(output) ~ r * ( vars[end] - output))
    return internals
end

@parameters t r
n = 12
@derivatives D'~t
@variables x(t) z(t)

eqs = [D(x) ~ - x]
internals = distributed_delay!(eqs, x, z, n, r)

de = ODESystem(eqs)
f = ODEFunction(de, vcat(x, internals, z), [r])


tspan = (0., 10.)
p = [3]
u0 = fill(0., length(de.dvs))
u0[1] = 1

prob = ODEProblem(f, u0, tspan, p)


### Hide all the internal variables that we used
mask =  @. !occursin("internal", String(Symbol(de.dvs)))
sol = solve(prob; save_idxs = mask)

plot(sol; label = ["input", "output"])

This may be better that nothing, but it still has a problem in that the n parameter must be set before creating the ODESystem and cannot be used by passing it to the ODEProblem. This would, for example, be irksome when performing parameter optimisation. I imagine that allowing n to be a normal parameter would either require the introduction of loops into our generated DE functions or that the model is instantiated only after the parameter vector has been passed to the DEProblem.

I consider this implementation to have four issues which should be addressed:

• the inelegant n parameter.
• the definition of the u0 vector (size is dependent on n). Also, this distributed delay formulation is not correct for all initial conditions.
• the cleaning up of internal variables in the solve function (easy).
• uniquely defining the internal variables so that one can have several distributed delays in the same model.

Thoughts?

Edit:
Re-instantiating the model in order to change the n parameter comes at a very real compilation cost.

[ChrisRackauckas]

I was thinking of it completely differently. Instead, a distributed delay is like

D(x(t)) ~ x(t-Gamma(lambda))

so the delay time is a random variable, and we need to interpret that and transform it to the set of ODEs.

[korsbo]

I recently wrote a paper on how gamma-distributed delays are often good for simplifying arbitrary linear signalling pathways and I now find myself craving an easy way to use them in more complicated models. I'm thus returning to this issue, initially to gauge how much work the implementation would take and whether it is wise to implement this during the little time I have left of my PhD.

I just spent some time familiarizing myself with ModelingToolkit but I can't claim to grok it yet.

I like the proposed notation except that the gamma distribution takes two parameters

D(x(t)) ~ x(t - Gamma(α, β))

At first glance, I don't think that it looks too hard to accommodate this API, what I'm more worried about are the next steps.

I would much prefer the parameter α to be an actual parameter which is passed to the DEProblems. However, since this parameter changes the number of variables of the system of DEs, this could get rather fiddly.

For an ODE, I would like

D(x(t)) ~ x(t - Gamma(α, β)) - d_x * x(t)

to generate a function:

(du, u, p, t)->begin
    ## Sort out the parameters first since they are used in the next step
    let (α, β, d_x) = (Int(p[1]), p[2], p[3])

        ## Make nice handles for derivatives and variables
        let (dx, dx_internals, x, x_internals) = (@view(du[1]), @view(du[2:1+α]), u[1], @view(u[2:1+α])) 

            ## The gamma distrubuted delay (with some bounds checks lest we error for low α values)
            if length(dx_internals) >= 1
                dx_internals[1] = β * (x - x_internals[1])
                if length(dx_internals) >= 2
                    @. dx_internals[2:end] = β * (x_internals[1:end-1] - x_internals[2:end])
                end
            else
                x_internals = [x]
            end

            ## The "normal" mathematics
            dx[1] = x_internals[end] - d_x * x
        end
    end
end

Simulating this (after assigning the function to the variable model) would look something like:

using DifferentialEquations
using Plots

p = [10, 2, 1.1]
tspan = (0., 20.)
u0 = vcat(1., fill(0., Int(p[1])))

prob = ODEProblem(model, u0, tspan, p)
sol = solve(prob)
plot(sol, alpha=0.2, label="")
plot!(sol, vars=[1], color=1)

The interface could be tidied up a bit by allowing automatic expansion of the u0 array in the ODEProblem overload, tweaking the solution object or the plot recipe to ignore the internal variables, and substituting the x_internals variable names with something that won't collide with user-defined variables or parameters.

However, I currently don't see how to nicely implement this since I don't see a simple way of making this work with the calculation of the Jacobian, etc.
Am I missing anything here?

A simpler, but less appealing alternative

A less appealing alternative would be to not allow α to be a free parameter but ensure that it is statically set before the construction of a new array of Equations. The transformation could then be something which took

eqs = [D(x(t)) ~ x(t - Gamma(3, β)) - d_x * x(t)]

and created some internal variables and then spat out

eqs = [
    D(x(t)) ~ x_internal[3] - d_x * x(t),
    D(x_internal[1]) ~ β * (x(t) - x_internal[1]),
    D(x_internal[2]) ~ β * (x_internal[1] - x_internal[2]),
    D(x_internal[3]) ~ β * (x_internal[2] - x_internal[3]),
]

Here, all the downstream calculations would be given for free, but one would have to create a whole new object every time one wished to change the α parameter.

Does anyone have a better idea of how to implement this? And if not, which of the two would be preferred?

[ChrisRackauckas]

I kind of like the second, since then it would work in general and all acceleration would work. But I understand the issue there. Maybe there can be DistributedDelaySystem for defining these equations, and the transformation function can then take in the parameters required to convert to an ODESystem, or if you compile it directly you get your first version?

[sdwfrost]

A distributed delay approach would be really handy right now - is there any progress on this?

[ChrisRackauckas]

Not on my end. But I like what @korsbo created here, so I would be happy to get it added to the library.

[sdwfrost]

Are 'standard' fixed delays (x(t-\tau)) implemented (or on the timeline)? For distributed delays, I'll see what I can do with @korsbo's suggestion. Have you seen this discrete time implementation that does bookkeeping of time spent in a certain state?

[ChrisRackauckas]

Not quite yet: we need to make a DDESystem that lowers properly. It's planned to be done, but not finished yet. There's a few tricks that need to be worked out in the function generation. That said, everything on the DifferentialEquations.jl side is already ready of course: it's just an MTK thing.

[sdwfrost]

Hi, just following up - not on the distributed delay, which I understand needs more work on the symbolic side, but on the fixed delay side (x(t-\tau)). Feel free to assign me something if I can help - I just need a pointer.

Here's a simple example DDE model:

https://github.com/epirecipes/sir-julia/blob/master/markdown/dde/dde.md

[ChrisRackauckas]

DDEs are fairly difficult. To illustrate this, I wrote down exactly what needs to be done:

#955

If you want to take a crack at it, see if you can get a good automatic detector for delay terms in the RHS, along with an @rule application for swapping to an h(p,t) term. We can take that discussion all to that issue. I don't think this one will take a day to crack, but a good dedicated 2 weeks.