Here we will create a PEtabODEproblem
for the Brannmark model, which is a model with requires pre-equilibration before comparing the model against data. In other words, to emulate a perturbation experiment the model must first reach a steady state where
To run the code, you will need the Brannmark PEtab files which can be found here. A fully runnable example of this tutorial is available here.
First, we load the necessary libraries and read the model.
using PEtab
using OrdinaryDiffEq
using Printf
path_yaml = joinpath(@__DIR__, "Brannmark", "Brannmark_JBC2010.yaml")
petab_model = PEtabModel(path_yaml, verbose=true)
PEtabModel for model Brannmark. ODE-system has 9 states and 23 parameters.
Generated Julia files are at ...
When dealing with pre-equilibration models, we must first reach a steady state :Rootfinding
, where we use any algorithm from NonlinearSolve.jl to find the roots of :Simulate
, where we simulate the model from the initial condition until it reaches a steady state. The latter method is more stable and often performs better.
When creating a PEtabODEProblem
, we can set steady-state solver options via SteadyStateSolver
. The first argument is the method to use, either :Rootfinding
or :Simulate
(recommended). For :Simulate
, we can choose how to terminate the steady-state simulation using the check_simulation_steady_state
argument, which accepts two options:
-
:wrms
: the weighted root-mean square$\sqrt{\sum_{i=1}^n \bigg( \frac{du[i]}{\mathrm{reltol}*u[i] + \mathrm{abstol}} \bigg) \frac{1}{n}} \leq 1$ , where$n$ is the number of ODEs. -
:Newton
: if the Newton-step$\Delta u$ is sufficiently small$\sqrt{\sum_{i=1}^n \bigg( \frac{\Delta u[i]}{\mathrm{reltol}*u[i] + \mathrm{abstol}} \bigg) \frac{1}{n}} \leq 1$ .
Newton often performs better, but it requires an invertible Jacobian. If the Jacobian is non-invertible, the code automatically switches to :wrms
. The default values for abstol
and reltol
are the tolerances of the ODE solver divided by 100.
In the example below, we use :Simulate
with :wrms
termination:
petab_problem = PEtabODEProblem(petab_model,
ode_solver=ODESolver(Rodas5P()),
ss_solver=SteadyStateSolver(:Simulate,
check_simulation_steady_state=:wrms),
gradient_method=:ForwardDiff)
p = petab_problem.θ_nominalT
gradient = zeros(length(p))
cost = petab_problem.compute_cost(p)
petab_problem.compute_gradient!(gradient, p)
@printf("Cost= %.2f\n", cost)
@printf("First element in the gradient = %.2e\n", gradient[1])
Cost = 141.89
First element in the gradient = 2.70e-03
Some useful notes regarding the steady-state solver:
- If you do not specify a
SteadyStateSolverOption
, the default option is:Simulate
with:wrms
. - You can also set a separate steady-state solver option for the gradient using
ss_solver_gradient
. - All gradient and hessian options are compatible with
:Simulate
. However,:Rootfinding
is only compatible with approaches that use forward-mode automatic differentiation.