Models with many conditions specific parameters
In this tutorial we show to create PEtabODEproblem
for a small ODE-model (≤20 states, ≤20 parameters) with many parameters to estimate (≈70), because most parameter are specific to a subset of simulation conditions. For example, in cond1 we have τcond1 and in cond2 we have τcond2 that maps to the ODE-system parameter τ respectively.
To run the code you need the Beer PEtab files which can be found here. A fully runnable example of this tutorial can be found here.
First we read the model and load necessary libraries.
using PEtab
using OrdinaryDiffEq
using Printf
pathYaml = joinpath(@__DIR__, "Beer", "Beer_MolBioSystems2014.yaml")
petabModel = readPEtabModel(pathYaml, verbose=true)
PEtabModel for model Beer. ODE-system has 4 states and 9 parameters.
Generated Julia files are at ...
Handling condition specific parameters
For a small ODE-system like Beer the most efficient gradient method is gradientMethod=:ForwardDiff
, and we can further compute the hessian via hessianMethod=:ForwardDiff
. However, there are several condition specific parameters, and as we by default we compute the gradient and hessian via a single call to ForwardDiff.jl it means that for this model we have to perform as many forward-passes (solve the ODE model) as there are model-parameters. To force one ForwardDiff.jl call per simulation condition we can use the option splitOverConditions=true
. Thus, for a model like Beer good options are:
odeSolverOptions
- Rodas5P() (works well for smaller models ≤ 15 states), and we use the defaultabstol, reltol .= 1e-8
.gradientMethod
- For small models like Beer forward mode automatic differentiation (AD) is fastest, thus below we choose:ForwardDiff
.hessianMethod
- For small models like Boehm with ≤20 parameters it is computationally feasible to compute the full Hessian via forward-mode AD. Thus, we choose:ForwardDiff
.splitOverConditions=true
- Force a call to ForwardDiff.jl per simulation condition. Most efficient for models where a majority of parameters are specific to a subset of simulation conditions.
odeSolverOptions = getODESolverOptions(Rodas5P(), abstol=1e-8, reltol=1e-8)
petabProblem = setupPEtabODEProblem(petabModel, odeSolverOptions,
gradientMethod=:ForwardDiff,
hessianMethod=:ForwardDiff,
splitOverConditions=true)
p = petabProblem.θ_nominalT
gradient = zeros(length(p))
hessian = zeros(length(p), length(p))
cost = petabProblem.computeCost(p)
petabProblem.computeGradient!(gradient, p)
petabProblem.computeHessian!(hessian, p)
@printf("Cost = %.2f\n", cost)
@printf("First element in the gradient = %.2e\n", gradient[1])
@printf("First element in the hessian = %.2f\n", hessian[1, 1])
Cost = -58622.91
First element in the gradient = 7.17e-02
First element in the hessian = 755266.33