# Estimating some parameters of DDEs with callbacks

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.

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

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)?

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)

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

``````
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
``````

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]
``````

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.

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)
``````
``````## 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])`.

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]...)'
``````

What’s the issue with doing that?

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?

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