Estimating some parameters of DDEs with callbacks

Specific Domains Modelling & Simulations
1

Hi,

I have a delay differential equation model with 3 model parameters and a drug perturbation with one parameter (alpha). I’ve set up the simulation as:

using DifferentialEquations
using DataFrames, DataFramesMeta
using CairoMakie

# Model parameters
τ = 10.0
k1 = 0.1
k2 = 0.1

# dose
α = 2.0
D = 0.0

p = [τ, k1, k2, α, D]

# Initial conditions
x0 = [1.0, 1.0]

# DDEs
function model(dx, x, xτ, p, t)
    τ, k1, k2, α, D = p
    xτ1 = xτ(p, t-τ)[2]
    dx[1] = 1/xτ1*(1 + α*D) - k1*x[1]
    dx[2] = x[1] - k2*x[2]
end

# Define dose
dose = [0.0, 1.0, 3.0]
dose_time = 100.0

# Simulation time
t_span = (0.0, 300.0)
t_step = 1.0

# Initialise delay
xτ(p, t) = ones(2)

# Define problem
prob = DDEProblem(model, x0, xτ, t_span, p; constant_lags=[τ])

# Loop over doses
sol = DataFrame()

for i in eachindex(dose)

    # Reset dose
    p[5] = 0.0

    # Callback
    function affect!(integrator)
        integrator.p[5] = dose[i]
    end
    cb_i = PresetTimeCallback(dose_time, affect!, save_positions=(false, false))

    # Solve
    sol_i = solve(prob, MethodOfSteps(Tsit5()), callback=cb_i, tstops=dose_time, saveat=t_step)

    # Store
    sol_i_df = @rtransform(DataFrame(sol_i), :dose = dose[i])

    sol = [sol; sol_i_df]
end

# Plot
function _plt()
    fig = Figure()
    ax1 = Axis(fig[1, 1])
    ax2 = Axis(fig[1, 2])

    for i in eachindex(dose)
        sol_i = @rsubset(sol, :dose == dose[i])
        lines!(ax1, sol_i.timestamp, sol_i.value1)
        lines!(ax2, sol_i.timestamp, sol_i.value2)
    end
    fig
end

_plt()

Now, I used this to generate dummy data as

# Generate dummy data
t_vec = hcat([sol[sol.dose.==i, :].timestamp for i in dose]...)'
data_mat = hcat([sol[sol.dose.==i, :].value2 .+ .5randn(size(t_vec,2)) for i in dose]...)'

Without the external forcing, I managed to estimate the parameters by writing a custom objective function along the lines of

 Custom objective function
function _optim_fcn(prob, data_t, data_obs, p0)

    function _cost_fcn(p, _)

        # Redefine problem
        new_prob = remake(prob; p=p)

        # Solve DE
        sol = solve(new_prob, MethodOfSteps(Tsit5()), saveat=data_t)

        # Calculate sum of squares
        loss = 0.0
        for i in eachindex(sol_fc)
            loss += (data_obs[i] - sol[i])^2
        end

        # Return
        return loss, sol_fc
    end

    # Optimize for the parameters that minimize the loss
    optf = OptimizationFunction(_cost_fcn, Optimization.AutoForwardDiff())
    optprob = OptimizationProblem(optf, p0, lb=zeros(length(p0)), ub=300.0 * ones(length(p0)))
    sol = solve(optprob, BFGS())

    # return the parameters we found
    return sol

end

But now, with the drug function across different doses and keeping alpha fixed, I’m struggling to implement a way to estimate the model parameters.
Any help would greatly be appreciated.

2

How so / why? Nothing you described should have issues with ForwardDiff.

3

I expected that this should not be a problem. Mainly I just don’t know how. How do I fix the alpha? How do I use the callback for the drug function in the optimisation and make sure that D which is technically defined as a parameter is not optimised? And how do I estimate the parameters based on time courses of different doses (and not just one dose at the time)?

4

You can write it so your optimization vector theta is a subset of the ODE p, i.e. p = [theta;D;alpha], and then solve etc. The D and alpha will be up-promoted to duals for type stability in the derivative calculation but that’s fine (and can be avoided by making p a tuple instead of a vector)

5

Ok. Can you give me a rewrite of the optimisation function that would achieve that?

6 

    function _cost_fcn(p, _)

        # Redefine problem
        new_prob = remake(prob; p=[p;D;alpha])

        # Solve DE
        sol = solve(new_prob, MethodOfSteps(Tsit5()), saveat=data_t)

        # Calculate sum of squares
        loss = 0.0
        for i in eachindex(sol_fc)
            loss += (data_obs[i] - sol[i])^2
        end

        # Return
        return loss, sol_fc
    end
7

Thanks. I now have

function model(dx, x, xτ, p, t)
    τ, k1, k2 = p[1]
    α, D = p[2:3]

    xτ1 = xτ(p, t - τ)[2]

    dx[1] = 1 / xτ1 * (1 + α * D) - k1 * x[1]
    dx[2] = x[1] - k2 * x[2]
end

# Generate dummy data for just one dose
t_vec = sol[sol.dose.==1.0, :].timestamp
data_mat = sol[sol.dose.==1.0, :].value2 .+ 0.5randn(size(t_vec, 2))

# Custom objective function
function _optim_fcn(prob, data_t, data_obs, p0, α, D)

    function _cost_fcn(p, _)

        # Redefine problem
        new_prob = remake(prob; p=[p; α; D])

        # Reset dose
        p[3] = 0.0

        # Callback
        function affect!(integrator)
            integrator.p[3] = D
        end
        cb_i = PresetTimeCallback(dose_time, affect!, save_positions=(false, false))

        # Solve DE
        sol = solve(new_prob, MethodOfSteps(Tsit5()), callback=cb_i, tstops=dose_time, saveat=data_t, save_idxs=2)

        # Calculate sum of squares
        loss = 0.0
        for i in eachindex(sol.u)
            loss += (data_obs[i] - sol.u[i])^2
        end

        # Return
        return loss, sol
    end

    # Optimize for the parameters that minimize the loss
    optf = OptimizationFunction(_cost_fcn, Optimization.AutoForwardDiff())
    optprob = OptimizationProblem(optf, p0, lb=zeros(length(p0)), ub=300.0 * ones(length(p0)))
    sol = solve(optprob, BFGS())

    # return the parameters we found
    return sol

end

p_estimate_i = _optim_fcn(prob, t_vec, data_mat, [10.0, 0.1, 0.1], 1.0, 1.0)

for which I receive the error:

ERROR: BoundsError: attempt to access ForwardDiff.Dual{ForwardDiff.Tag{Optimization.var"#169#186"{OptimizationFunction{true, Optimization.AutoForwardDiff{nothing}, var"#_cost_fcn#989"{DDEProblem{Vector{Float64}, Tuple{Float64, Float64}, Vector{Float64}, Tuple{}, true, Vector{Any}, DDEFunction{true, SciMLBase.FullSpecialize, typeof(model), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, typeof(xτ), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, Vector{Float64}, Vector{Float64}, Float64, Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Optimization.ReInitCache{Vector{Float64}, SciMLBase.NullParameters}}, Float64}, Float64, 3} at index [2]
8

When I try to run that code I get sol is not defined, and it’s clearly not at tvec because just a function above is defined. Can you please share a full example.

9

Full working example:

## Optimisation of parameters for DDEs with callback --------------------------
using DifferentialEquations
using DataFrames, DataFramesMeta
using CairoMakie
using Optimization, OptimizationOptimJL

# Model parameters
τ = 10.0
k1 = 0.1
k2 = 0.1

pm = [τ, k1, k2]

# dose
α = 2.0
D = 0.0

p = [pm, α, D]

# Initial conditions
x0 = [1.0, 1.0]

# DDEs
function model(dx, x, xτ, p, t)
    τ, k1, k2 = p[1]
    α, D = p[2:3]

    xτ1 = xτ(p, t - τ)[2]

    dx[1] = 1 / xτ1 * (1 + α * D) - k1 * x[1]
    dx[2] = x[1] - k2 * x[2]
end

# Define dose
dose = [0.0, 1.0, 3.0]
dose_time = 100.0

# Simulation time
t_span = (0.0, 300.0)
t_step = 1.0

# Initialise delay
xτ(p, t) = ones(2)

# Define problem
prob = DDEProblem(model, x0, xτ, t_span, p; constant_lags=[τ])

# Loop over doses
sol = DataFrame()

for i in eachindex(dose)

    # Reset dose
    p[3] = 0.0

    # Callback
    function affect!(integrator)
        integrator.p[3] = dose[i]
    end
    cb_i = PresetTimeCallback(dose_time, affect!, save_positions=(false, false))

    # Solve
    sol_i = solve(prob, MethodOfSteps(Tsit5()), callback=cb_i, tstops=dose_time, saveat=t_step)

    # Store
    sol_i_df = @rtransform(DataFrame(sol_i), :dose = dose[i])

    sol = [sol; sol_i_df]
end

# Generate dummy data
# t_vec = hcat([sol[sol.dose.==i, :].timestamp for i in dose]...)'
# data_mat = hcat([sol[sol.dose.==i, :].value2 .+ .5randn(size(t_vec,2)) for i in dose]...)'

# Generate dummy data for just one dose
t_vec = sol[sol.dose.==1.0, :].timestamp
data_mat = sol[sol.dose.==1.0, :].value2 .+ 0.5randn(size(t_vec, 2))

# Custom objective function
function _optim_fcn(prob, data_t, data_obs, p0, α, D)

    function _cost_fcn(p, _)

        # Redefine problem
        new_prob = remake(prob; p=[p; α; D])

        # Reset dose
        p[3] = 0.0

        # Callback
        function affect!(integrator)
            integrator.p[3] = D
        end
        cb_i = PresetTimeCallback(dose_time, affect!, save_positions=(false, false))

        # Solve DE
        sol = solve(new_prob, MethodOfSteps(Tsit5()), callback=cb_i, tstops=dose_time, saveat=data_t, save_idxs=2)

        # Calculate sum of squares
        loss = 0.0
        for i in eachindex(sol.u)
            loss += (data_obs[i] - sol.u[i])^2
        end

        # Return
        return loss, sol
    end

    # Optimize for the parameters that minimize the loss
    optf = OptimizationFunction(_cost_fcn, Optimization.AutoForwardDiff())
    optprob = OptimizationProblem(optf, p0, lb=zeros(length(p0)), ub=300.0 * ones(length(p0)))
    sol = solve(optprob, BFGS())

    # return the parameters we found
    return sol

end

p_estimate_i = _optim_fcn(prob, t_vec, data_mat, [10.0, 0.1, 0.1], 1.0, 1.0)
10 
## Optimisation of parameters for DDEs with callback --------------------------
using DifferentialEquations
using DataFrames, DataFramesMeta
using Optimization, OptimizationOptimJL

# Model parameters
τ = 10.0
k1 = 0.1
k2 = 0.1

pm = [τ, k1, k2]

# dose
α = 2.0
D = 0.0

p = [pm, α, D]

# Initial conditions
x0 = [1.0, 1.0]

# DDEs
function model(dx, x, xτ, p, t)
    τ, k1, k2 = p[1]
    α, D = p[2:3]

    xτ1 = xτ(p, t - τ)[2]

    dx[1] = 1 / xτ1 * (1 + α * D) - k1 * x[1]
    dx[2] = x[1] - k2 * x[2]
end

# Define dose
dose = [0.0, 1.0, 3.0]
dose_time = 100.0

# Simulation time
t_span = (0.0, 300.0)
t_step = 1.0

# Initialise delay
xτ(p, t) = ones(2)

# Define problem
prob = DDEProblem(model, x0, xτ, t_span, p; constant_lags=[τ])

# Loop over doses
sol = DataFrame()

for i in eachindex(dose)

    # Reset dose
    p[3] = 0.0

    # Callback
    function affect!(integrator)
        integrator.p[3] = dose[i]
    end
    cb_i = PresetTimeCallback(dose_time, affect!, save_positions=(false, false))

    # Solve
    sol_i = solve(prob, MethodOfSteps(Tsit5()), callback=cb_i, tstops=dose_time, saveat=t_step)

    # Store
    sol_i_df = @rtransform(DataFrame(sol_i), :dose = dose[i])

    sol = [sol; sol_i_df]
end

# Generate dummy data
# t_vec = hcat([sol[sol.dose.==i, :].timestamp for i in dose]...)'
# data_mat = hcat([sol[sol.dose.==i, :].value2 .+ .5randn(size(t_vec,2)) for i in dose]...)'

# Generate dummy data for just one dose
t_vec = sol[sol.dose.==1.0, :].timestamp
data_mat = sol[sol.dose.==1.0, :].value2 .+ 0.5randn(size(t_vec, 2))

# Custom objective function
p0, α, D = [10.0, 0.1, 0.1], 1.0, 1.0

function _cost_fcn(p, _)

       # Redefine problem
       new_prob = remake(prob; p=[p, α, D])

       # Reset dose
       p[3] = 0.0

       # Callback
       function affect!(integrator)
       integrator.p[3] = D
       end
       cb_i = PresetTimeCallback(dose_time, affect!, save_positions=(false, false))

       # Solve DE
       sol = solve(new_prob, MethodOfSteps(Tsit5()), callback=cb_i, tstops=dose_time, saveat=t_vec, save_idxs=2)

       # Calculate sum of squares
       loss = 0.0
       for i in eachindex(sol.u)
       loss += (data_mat[i] - sol.u[i])^2
       end

       # Return
       return loss, sol
end

# Optimize for the parameters that minimize the loss
optf = OptimizationFunction(_cost_fcn, Optimization.AutoForwardDiff())
optprob = OptimizationProblem(optf, p0, lb=zeros(length(p0)), ub=300.0 * ones(length(p0)))
sol = solve(optprob, BFGS(initial_stepnorm = 0.01))

Works. Basically the issue is just that you don’t make your parameters the same. remake(prob; p=[p; α; D]) should be changed to remake(prob; p=[p, α, D]).

11

D’oh. Thanks a lot. Could you also suggest a way to not only fit the data for one dose but for all doses simultaneously, i.e. using

t_vec = hcat([sol[sol.dose.==i, :].timestamp for i in dose]...)'
data_mat = hcat([sol[sol.dose.==i, :].value2 .+ .5randn(size(t_vec,2)) for i in dose]...)'
12

What’s the issue with doing that?

13

Figured it out. Thanks anyway.
Maybe as a last question: I regularly get the error:

AssertionError: isfinite(phi_d) && isfinite(gphi)

I read that I could get rid of it using BFGS(linesearch=LineSearches.BackTracking()) or LBFGS() but I don’t quite understand why and what the difference between these solutions is or which should be preferred. Could you give some insight or link to a good resource?

14

That happens when the integrator takes a step to where the parameters are infinite, i.e. where your model starts to diverge. Optim.jl is particularly prone to this with its default BFGS so I typically don’t recommend using Optim for the optimization and instead recommend using one of the NLopt BFGS implementations which seems more robust.

1 Like