Using optimizers directly

A PEtabODEProblem contains all information needed to use an optimizer directly for parameter estimation. While PEtab.jl provides high-level functions for single-start estimation (calibrate) and multi-start estimation (calibrate_multistart), direct use can be useful for unsupported optimization packages and/or custom workflows.

This tutorial demonstrates how to use the IPNewton method from Optim.jl directly via a PEtabODEProblem. As a running example, the Michaelis–Menten model from the starting tutorial is used.

julia

using Catalyst, PEtab
rn = @reaction_network begin
    @parameters S0 c3=3.0
    @species begin
        S(t) = S0
        E(t) = 50.0
        SE(t) = 0.0
        P(t) = 0.0
    end
    @observables begin
        obs1 ~ S + E
        obs2 ~ P
    end
    c1, S + E --> SE
    c2, SE --> S + E
    c3, SE --> P + E
end

# Observables
@parameters sigma
petab_obs1 = PEtabObservable(:petab_obs1, :obs1, 3.0)
petab_obs2 = PEtabObservable(:petab_obs2, :obs2, sigma)
observables = [petab_obs1, petab_obs2]

# Parameters to estimate
p_c1 = PEtabParameter(:c1)
p_c2 = PEtabParameter(:c2)
p_S0 = PEtabParameter(:S0)
p_sigma = PEtabParameter(:sigma)
pest = [p_c1, p_c2, p_S0, p_sigma]

# Measurements; simulate with 'true' parameters
using DataFrames, OrdinaryDiffEqRosenbrock
ps = [:c1 => 1.0, :c2 => 10.0, :c3 => 1.0, :S0 => 100.0]
u0 = [:E => 50.0, :SE => 0.0, :P => 0.0]
tspan = (0.0, 10.0)
oprob = ODEProblem(rn, u0, tspan, ps)
sol = solve(oprob, Rodas5P(); saveat = 0:0.5:10.0)
obs1 = (sol[:S] + sol[:E]) .+ randn(length(sol[:E]))
obs2   = sol[:P] .+ randn(length(sol[:P]))
df1 = DataFrame(obs_id = "petab_obs1", time = sol.t, measurement = obs1)
df2 = DataFrame(obs_id = "petab_obs2", time = sol.t, measurement = obs2)
measurements = vcat(df1, df2)

# Create the PEtabODEProblem
petab_model = PEtabModel(rn, observables, measurements, pest)
petab_prob = PEtabODEProblem(petab_model)

PEtabODEProblem ReactionSystemModel: 4 parameters to estimate
(for more statistics, call `describe(petab_prob)`)

Extracting relevant functions from a PEtabODEProblem

A numerical optimizer requires an objective function, often a gradient function, and sometimes a Hessian function. PEtab.jl works with likelihoods, so the objective corresponds to the negative log-likelihood (nllh):

julia

x = get_x(petab_prob)
nllh = petab_prob.nllh(x; prior = true)

21508.98163982404

The keyword argument prior = true (default) includes parameter priors (if defined). In addition, the PEtabODEProblem provides both in-place and out-of-place gradients:

julia

g_inplace = similar(x)
petab_prob.grad!(g_inplace, x)
g_outplace = petab_prob.grad(x)

ComponentVector{Float64}(log10_c1 = -174.85991323811464, log10_c2 = 173.5889319511079, log10_S0 = 168613.00346875747, log10_sigma = 24.715040093970234)

and Hessians:

julia

h_inplace = zeros(length(x), length(x))
petab_prob.hess!(h_inplace, x)
h_outplace = petab_prob.hess(x)

4×4 Matrix{Float64}:
   484.914     -481.687     -1005.76        -0.137188
  -481.687      481.387      1008.98         0.136463
 -1005.76      1008.98          1.0961e6  -117.111
    -0.137188     0.136463   -117.111      108.863

The input x is typically a ComponentVector, but Vector inputs are also supported. More details on what is available in a PEtabODEProblem can be found in the API documentation

Lastly, for ODE models, parameter bounds are often important: without bounds, the optimizer often explores regions where the ODE solver fails, increasing runtime [2]. The bounds in petab_probcan be accessed with:

julia

lb, ub = petab_prob.lower_bounds, petab_prob.upper_bounds

lb and ub are ComponentVectors. If an optimizer does not support ComponentVector, convert them to Vectors with collect (see below).

Wrapping Optim.jl `IPNewton`

Optim.jl’s IPNewton expects an objective, gradient, Hessian, and box constraints (bounds provided as plain vectors). Using the functions in a PEtabODEProblem (see above), this can be set up as:

julia

using Optim
x0 = collect(get_x(petab_prob))
df  = TwiceDifferentiable(
    petab_prob.nllh, petab_prob.grad!, petab_prob.hess!, x0
)
dfc = TwiceDifferentiableConstraints(
    collect(lb), collect(ub)
)

Here, collect converts ComponentVectors to Vectors for Optim.jl. The optimization can then be run from x0:

julia

opt_res = Optim.optimize(df, dfc, x0, IPNewton())

 * Status: success

 * Candidate solution
    Final objective value:     1.570987e+02

 * Found with
    Algorithm:     Interior Point Newton

 * Convergence measures
    |x - x'|               = 0.00e+00 ≤ 0.0e+00
    |x - x'|/|x'|          = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
    |g(x)|                 = 5.12e-03 ≰ 1.0e-08

 * Work counters
    Seconds run:   1  (vs limit Inf)
    Iterations:    39
    f(x) calls:    87
    ∇f(x) calls:   87
    ∇f(x)ᵀv calls: 0
    ∇²f(x) calls:  39
    ∇²f(x)v calls: 0

References

F. Fröhlich and P. K. Sorger. Fides: Reliable trust-region optimization for parameter estimation of ordinary differential equation models. PLoS computational biology 18, e1010322 (2022).

Using optimizers directly ​

Extracting relevant functions from a PEtabODEProblem ​

Wrapping Optim.jl IPNewton ​

References ​

Using optimizers directly

Extracting relevant functions from a PEtabODEProblem

Wrapping Optim.jl `IPNewton`

References