PEtab fails on simple (?) problem. Unable to compute cost function for true value

[TorkelE]

I have a relatively simple case of parameter fitting where I try to use PEtab, but unfortunately, it fails. It might be that I have not used a good noise formula or similar. However, there are some warnings when I try to evaluate the cost function for the true parameter set which makes me suspicious.

Code

The model is from a catalysis chain network used for e.g. chemical synthesis. Parameter values and initial conditions are at least somewhat realistic.

Generate data.

using Catalyst
# Create model.
rn = @reaction_network begin
    k1, LPd + ArBr --> LPdArBr
    k2, LPdArBr + RB --> LPdArR + BBr
    k3, LPdArR --> LPd + ArR
end

# Prepare data generation simulation.
u0 = [:LPd => 0.008, :ArBr => 0.08, :LPdArBr => 0.0, :RB => 0.2, :LPdArR => 0.0, :BBr => 0.0, :ArR => 0.0]
p_true = [:k1 => 50000.0, :k2 => 1000.0, :k3 => 12000.0] 

# Generate synthetic data.
using OrdinaryDiffEq
oprob_true = ODEProblem(rn, u0, (0.0, 0.1), p_true)
true_sol = solve(oprob_true, Tsit5())
data_sol = solve(oprob_true, Tsit5(); saveat=0.01)
data_ts = data_sol.t[2:end]
data_vals = (0.8 .+ 0.4*rand(10)) .* data_sol[:ArR][2:end]

# Plots the true solutions and the (synthetic data) measurements.
using Plots
plot(true_sol; idxs=:ArR, label="True solution", lw=8)
plot!(data_ts, data_vals; label="Measurements", seriestype=:scatter, ms=6, color=:blue)

Fit model

using PEtab
@unpack ArR = rn
obs_ArR = PEtabObservable(ArR, ArR*0.5)
observables = Dict("obs_ArR" => obs_ArR)

# Parameters.
par_k1 = PEtabParameter(:k1)
par_k2 = PEtabParameter(:k2)
par_k3 = PEtabParameter(:k3)
params = [par_k1, par_k2, par_k3]

# Measurements.
using DataFrames
measurements = DataFrame(obs_id="obs_ArR", time=data_ts, measurement=data_vals)

# Create PEtab problem.
petab_model = PEtabModel(rn, observables, measurements, params; state_map=u0)
petab_problem = PEtabODEProblem(petab_model)

# Fit model to data.
using Optim
@time res_ms = calibrate_model_multistart(petab_problem, IPNewton(), 50, "OptimisationRuns"; save_trace=true)

# Plot the fit.
plot(res_ms, petab_problem; observable_ids=["obs_ArR"])

plot(res_ms; plot_type=:best_objective)

Weird stuff

Here, this code

petab_problem.compute_cost(last.(p_true)) # Returns Inf.

generates these warnings:

┌ Warning: Automatic dt set the starting dt as NaN, causing instability. Exiting.
└ @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/2nLli/src/solve.jl:571
┌ Warning: NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.
└ @ SciMLBase ~/.julia/packages/SciMLBase/LBy2M/src/integrator_interface.jl:577
┌ Warning: Failed to solve ODE model
└ @ PEtab ~/.julia/packages/PEtab/fBBay/src/Compute_cost/Compute_cost.jl:79

Which seems weird. Surely the original parameter set should at least be simulatable? It seems PEtab is trying to determine an initial step size for this parameter set, but fails, resulting in a failed simulation, and infinite cost.

[sebapersson]

Two things:

I think the reason for the bad performance in the parameter estimation is the bounds, by default they are [1e-3, 1e3] which you true parameter values exceed (typically for most sysbio models we stay withing these bounds hence they are the default), bounds can be set in PEtabParameter

Your second problem is because by default a PEtabParameter is assumed to be on the log10-scale. The following should work:

petab_problem.compute_cost(log10.(last.(p_true)))

Hopefully these two things fix the problem!

[TorkelE]

You are correct, on both. Thanks a lot! :D