# Estimating some parameters of DDEs with callbacks

**URL:** <https://discourse.julialang.org/t/estimating-some-parameters-of-ddes-with-callbacks/97614>\
**Category:** Modelling & Simulations\
**Tags:** question\
**Created:** [April 18, 2023, 12:39pm UTC](https://discourse.julialang.org/t/estimating-some-parameters-of-ddes-with-callbacks/97614 "2023-04-18T12:39:31Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![pascal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pascal/32/207922_2.png) [@pascal](https://discourse.julialang.org/u/pascal)\
**Post date:** [April 18, 2023, 12:39pm UTC](https://discourse.julialang.org/t/estimating-some-parameters-of-ddes-with-callbacks/97614/1 "2023-04-18T12:39:31Z")

</div>

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:

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

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

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 19, 2023, 7:29am UTC](https://discourse.julialang.org/t/estimating-some-parameters-of-ddes-with-callbacks/97614/2 "2023-04-19T07:29:44Z")

</div>

> [@pascal](#):
>
> 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.

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

---

<div class="post-metadata">

**Author:** ![pascal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pascal/32/207922_2.png) [@pascal](https://discourse.julialang.org/u/pascal)\
**Post date:** [April 19, 2023, 7:50am UTC](https://discourse.julialang.org/t/estimating-some-parameters-of-ddes-with-callbacks/97614/3 "2023-04-19T07:50:37Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 19, 2023, 7:55am UTC](https://discourse.julialang.org/t/estimating-some-parameters-of-ddes-with-callbacks/97614/4 "2023-04-19T07:55:37Z")

</div>

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)

---

<div class="post-metadata">

**Author:** ![pascal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pascal/32/207922_2.png) [@pascal](https://discourse.julialang.org/u/pascal)\
**Post date:** [April 19, 2023, 9:12am UTC](https://discourse.julialang.org/t/estimating-some-parameters-of-ddes-with-callbacks/97614/5 "2023-04-19T09:12:21Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 19, 2023, 1:46pm UTC](https://discourse.julialang.org/t/estimating-some-parameters-of-ddes-with-callbacks/97614/6 "2023-04-19T13:46:31Z")

</div>

```julia

    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

```

---

<div class="post-metadata">

**Author:** ![pascal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pascal/32/207922_2.png) [@pascal](https://discourse.julialang.org/u/pascal)\
**Post date:** [April 20, 2023, 8:42am UTC](https://discourse.julialang.org/t/estimating-some-parameters-of-ddes-with-callbacks/97614/7 "2023-04-20T08:42:39Z")

</div>

Thanks. I now have

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

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 22, 2023, 11:39am UTC](https://discourse.julialang.org/t/estimating-some-parameters-of-ddes-with-callbacks/97614/8 "2023-04-22T11:39:11Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![pascal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pascal/32/207922_2.png) [@pascal](https://discourse.julialang.org/u/pascal)\
**Post date:** [April 24, 2023, 9:03am UTC](https://discourse.julialang.org/t/estimating-some-parameters-of-ddes-with-callbacks/97614/9 "2023-04-24T09:03:24Z")

</div>

Full working example:

```julia
## 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)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 24, 2023, 1:47pm UTC](https://discourse.julialang.org/t/estimating-some-parameters-of-ddes-with-callbacks/97614/10 "2023-04-24T13:47:55Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![pascal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pascal/32/207922_2.png) [@pascal](https://discourse.julialang.org/u/pascal)\
**Post date:** [April 24, 2023, 2:19pm UTC](https://discourse.julialang.org/t/estimating-some-parameters-of-ddes-with-callbacks/97614/11 "2023-04-24T14:19:05Z")

</div>

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

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 24, 2023, 2:50pm UTC](https://discourse.julialang.org/t/estimating-some-parameters-of-ddes-with-callbacks/97614/12 "2023-04-24T14:50:58Z")

</div>

What’s the issue with doing that?

---

<div class="post-metadata">

**Author:** ![pascal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pascal/32/207922_2.png) [@pascal](https://discourse.julialang.org/u/pascal)\
**Post date:** [April 25, 2023, 7:56am UTC](https://discourse.julialang.org/t/estimating-some-parameters-of-ddes-with-callbacks/97614/13 "2023-04-25T07:56:52Z")

</div>

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

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 25, 2023, 12:00pm UTC](https://discourse.julialang.org/t/estimating-some-parameters-of-ddes-with-callbacks/97614/14 "2023-04-25T12:00:28Z")

</div>

> [@pascal](#):
>
> 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.