PEtab.jl API

PEtabModel

A PEtabModel for model fitting can be imported from the PEtab standard format, or defined directly in Julia. In the Julia interface, observables that link model outputs to data are specified with PEtabObservable, estimated parameters with PEtabParameter, simulation conditions with PEtabConditions, and events/callbacks with PEtabEvent.

PEtab.PEtabObservable Type

PEtabObservable(observable_id, observable_formula, noise_formula; distribution = Normal)

Observation model linking model output (observable_formula) to measurement data via a likelihood defined by distribution with noise/scale given by noise_formula.

For examples, see the online package documentation.

Arguments

• observable_id::Union{String, Symbol}: Observable identifier. Measurements are linked to this observable via the column obs_id in the measurement table.

• observable_formula: Observable expression. Two supported forms:

  • Model-observable: A Symbol identifier matching an observable defined in a Catalyst ReactionSystem @observables block, or a non-differential variable defined in a ModelingToolkit ODESystem @variables block.

  • expression: A String, :Symbol, Real, or a Symbolics expression (Num). Can include standard Julia functions (e.g. exp, log, sin, cos). Variables may reference model species, model parameters, or PEtabParameters. Can include time-point-specific parameters (see documentation for examples).

• noise_formula: Noise/scale expression (String or Symbolics equation), same rules as observable_formula. May include time-point-specific parameters (see documentation).

Keyword Arguments

• distribution: Measurement noise distribution. Valid options are Normal (default), Laplace, LogNormal, Log2Normal, Log10Normal and LogLaplace. See below for mathematical definition.

Mathematical description

For a measurement m, model output y = observable_formula, and a noise parameter σ = noise_formula, PEtabObservable defines the likelihood linking the model output to the measurement data: $\pi (m\mid y,\sigma )$.

For distribution = Normal, the measurement is assumed to be normally distributed with $m\sim \mathcal{N}(y,{\sigma }^2)$. The likelihood formula is:

$$\pi (m\mid y,\sigma )=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp (-\frac{(m-y{)}^2}{2{\sigma }^2})$$

If $\sigma =1$, this likelihood reduces to the least-squares objective function.

For distribution = Laplace, the measurement is assumed to be Laplace distributed with $m\sim \mathcal{L}(y,\sigma )$. The likelihood formula is:

$$\pi (m\mid y,\sigma )=\frac{1}{2\sigma }\exp (-\frac{|m-y|}{\sigma })$$

For distribution = LogNormal, the log of the measurement is assumed to be Normal distributed with $\log (m)\sim \mathcal{N}(\log (y),{\sigma }^2)$ (requires $m>0$ and $y>0$). The likelihood formula is:

$$\pi (m\mid y,\sigma )=\frac{1}{\sqrt{2\pi {\sigma }^2}m}\exp (-\frac{(\log (m)-\log (y){)}^2}{2{\sigma }^2})$$

For distribution = Log2Normal|Log10Normal, similar to the LogNormal, ${\log }_2(m)$ and ${\log }_{10}(m)$ are assumed to be normally distributed.

For distribution = LogLaplace, the log of the measurement is assumed to be Laplace distributed with $\log (m)\sim \mathcal{L}(\log (y),\sigma )$ (requires $m>0$ and $y>0$). The likelihood formula is:

$$\pi (m\mid y,\sigma )=\frac{1}{2\sigma m}\exp (-\frac{|\log (m)-\log (y)|}{\sigma })$$

For numerical stability, PEtab.jl works with the log-likelihood in practice.

source

PEtab.PEtabParameter Type

PEtabParameter(parameter_id; kwargs...)

Parameter-estimation data for parameter parameter_id (bounds, scale, prior, and whether to estimate).

All parameters estimated in a PEtabODEProblem must be declared as PEtabParameter. parameter_id must correspond to a model parameter and/or a parameter appearing in an observable_formula or noise_formula of a PEtabObservable.

Keyword Arguments

• scale::Symbol = :log10: Scale the parameter is estimated on. One of :log10 (default), :log2, :log, or :lin. Estimating on a log scale often improves performance and is recommended.

• lb: Lower bound, specified on the linear scale; e.g. with scale = :log10, pass lb = 1e-3, not log10(1e-3). Defaults to 1e-3 without a prior, otherwise to the lower bound of the prior support.

• ub: Upper bound, same convention as lb. Defaults to 1e3 without a prior, otherwise to the upper bound of the prior support.

• prior = nothing: Optional prior distribution acting on the linear parameter scale (i.e. even if scale = :log10, the prior is on x, not on log10(x)). If the prior’s support extends beyond provided lb/ub bounds, it is truncated by [lb, ub]. Any continuous univariate distribution from Distributions.jl is supported, including truncated distributions.

• estimate::Bool = true: Whether the parameter is estimated (default true) or treated as a constant (false).

• value = nothing: Value used when estimate = false, and the value returned by get_x. Defaults to the midpoint of [lb, ub].

Priors and Parameter Estimation

If at least one parameter in a PEtabODEProblem has a prior specified, parameter estimation uses a maximum-a-posteriori (MAP) objective:

$$\underset{\mathbf{x}}{min}-\ell (\mathbf{x})-\sum _{i\in \mathcal{I}}\log {\pi }_i(x_i)$$

where $\mathcal{I}$ indexes parameters with an explicit prior density ${\pi }_i$. If no parameter has a prior, parameter estimation reduces to the maximum-likelihood (ML) objective $-\ell (\mathbf{x})$.

source

PEtab.PEtabCondition Type

PEtabCondition(condition_id, assignments::Pair...; t0 = 0.0)

Simulation condition that overrides model entities according to assignments under condition_id.

Used to set initial values and/or model parameters for different experimental conditions. For examples, see the online package documentation.

Arguments

• condition_id::Union{String, Symbol}: Simulation condition identifier. Measurements are linked to this condition via the column simulation_id in the measurement table.

• assignments: One or more assignments of the form target_id => target_value.

  • target_id: Entity to assign (Symbol, String, or Symbolics Num). Can be a model state id or model parameter id for a parameter that is not estimated.

  • target_value: Value/expression assigned to target_id. A String, Real, or a Symbolics expression (Num) which can use standard Julia functions (e.g. exp, log, sin, cos). Any variables referenced must be model parameters or PEtabParameters (model state variables are not allowed).

Keyword Arguments

• t0 = 0.0: Simulation start time for the condition.

source

PEtab.PEtabEvent Type

PEtabEvent(condition, assignments::Pair...; condition_ids = [:all])

Model event triggered when condition transitions from false to true, applying the updates in assignments.

For examples, see the online package documentation.

Arguments

• condition: Boolean expression that triggers the event on a false → true transition. For example, t == 3.0 triggers at t = 3.0; S > 2.0 triggers when species S crosses 2.0 from below.

• assignments: One or more updates of the form target_id => target_value.

  • target_id: Entity to update (Symbol, String, or Symbolics Num). May be a model state id or model parameter id.

  • target_value: Value/expression assigned to target_id (Real, String, or Symbolics Num). May use standard Julia functions (e.g. exp, log, sin, cos) and reference model states/parameters.

Keyword Arguments

• condition_ids: Simulation condition identifiers (as declared by PEtabCondition) for which the event applies. If [:all] (default), the event applies to all conditions.

Event evaluation order

target_value expressions are evaluated at the condition trigger point using pre-event model values, meaning all assignments are applied simultaneously (updates do not see each other’s new values). If a time-triggered event fires at the same time as a measurement, the model observable is evaluated after applying the event.

source

Given a dynamic model (as a ReactionSystem or ODESystem), measurement data as a DataFrame, and optional simulation conditions, a PEtabModel can be created with:

PEtab.PEtabModel Type

PEtabModel(path_yaml; kwargs...)

Import a PEtab problem in the standard (YAML + tables) format from path_yaml as a PEtabModel for parameter estimation.

Keyword Arguments

• ifelse_to_callback::Bool = true: Rewrite ifelse (SBML piecewise) expressions as callbacks. Typically improves simulation performance.

• write_to_file::Bool = false: Write generated Julia functions to dirname(path_yaml)/Julia_model_files/ (useful for debugging).

• verbose::Bool = false: Print progress while building the model.

source

PEtabModel(sys, observables, measurements::DataFrame, parameters; kwargs...)

Create a PEtabModel for parameter estimation from model system sys, observables linking model output to measurements, and parameters to estimate.

If there are multiple observables, parameters, simulation_conditions, and/or events, pass them as a Vector (e.g. Vector{PEtabObservable}).

For examples, see the online package documentation.

Arguments

• sys: Model system (ReactionSystem or ODESystem).

• observables: PEtabObservable(s) linking model output to measurements.

• measurements: Measurement table (see documentation for required columns).

• parameters: PEtabParameter(s) to estimate.

Keyword Arguments

• simulation_conditions = nothing: Optional PEtabCondition(s) specifying condition-specific overrides (initial values and/or model parameters).

• events = nothing: Optional PEtabEvent(s) defining model events/callbacks.

• speciemap: Optional vector of pairs [:state_id => value, ...] setting default initial values for model states/species. Only needed if values are not already defined in the model system (recommended) or provided elsewhere in the PEtab problem.

• parametermap: Like speciemap, but for model parameters.

• verbose::Bool = false: Print progress while building the model.

See also: PEtabObservable, PEtabParameter, PEtabCondition, and PEtabEvent.

source

PEtabODEProblem

From a PEtabModel, a PEtabODEProblem can be created with:

PEtab.PEtabODEProblem Type

PEtabODEProblem(model::PEtabModel; kwargs...) -> PEtabODEProblem

Create a PEtabODEProblem from model for parameter estimation and Bayesian inference.

If no keyword arguments are provided, PEtab.jl chooses defaults for the ODESolver, gradient method, and Hessian method (see the online documentation on default options).

Keyword Arguments

• odesolver::ODESolver: ODE solver options used when evaluating the objective (negative log-likelihood).

• odesolver_gradient::ODESolver: ODE solver options used when evaluating gradients. Defaults to odesolver.

• ss_solver::SteadyStateSolver: Steady-state solver options used when evaluating the objective. Only applicable for models with steady-state (pre-equilibration) simulations.

• ss_solver_gradient::SteadyStateSolver: Steady-state solver options used when evaluating gradients. Defaults to ss_solver.

• gradient_method::Symbol: Method used to compute gradients. Available options and defaults are described in the online documentation.

• hessian_method::Symbol: Method used to compute Hessians. Available options and defaults are described in the online documentation.

• FIM_method = nothing: Method used to compute the empirical Fisher Information Matrix (FIM). Accepts the same options as hessian_method. If nothing, the default is used.

• sparse_jacobian::Bool = false: Whether to use a sparse Jacobian when solving the ODE. Can substantially improve performance for large models.

• sensealg = nothing: Sensitivity algorithm used for gradient computations. Available and default options depend on gradient_method. See the documentation for details.

• chunksize = nothing: Chunk size used by ForwardDiff.jl for forward-mode automatic differentiation (gradients and Hessians). If not provided, a default is chosen. Tuning chunksize can improve performance, but the optimal value is model-dependent.

• split_over_conditions::Bool = false: Whether to split ForwardDiff-based derivative computations across simulation conditions. Can improve performance for models with many condition-specific parameters, but adds overhead otherwise.

• reuse_sensitivities::Bool = false: Whether to reuse forward sensitivities from gradient computations when forming a Gauss-Newton Hessian approximation. Only applies when hessian_method = :GaussNewton and gradient_method = :ForwardEquations. This can substantially improve performance when the optimizer evaluates gradient and Hessian together.

• verbose::Bool = false: Whether to print progress while building the PEtabODEProblem.

Returns

A PEtabODEProblem, which provides everything needed for wrapping an optmization algorithm, key fields are:

• nllh: Negative log-likelihood function (or negative log-posterior if priors are included); nllh(x).

• grad!: In-place gradient function; grad!(g, x).

• grad: Out-of-place gradient function; g = grad(x).

• nllh_grad: Function to compute nllh and grad simultanesouly: nllh, g = nllh_grad(x). More efficient than calling nllh and grad separately since quantities from computing grad are resued for nllh.

• hess!: In-place Hessian function; hess!(H, x).

• hess: Out-of-place Hessian function; H = hess(x).

• FIM: Empirical Fisher Information Matrix function; F = FIM(x).

• chi2: Chi-squared test statistic function; χ² = chi2(x).

• residuals: Residual vector function; r = residuals(x).

• simulated_values: Model-simulated values corresponding to each measurement row, in the same order as the measurement table.

• lower_bounds, upper_bounds: Parameter bounds used for optimization/inference.

Unless otherwise noted, the input x can be a Vector or a ComponentArray (the output type matches the input type where applicable). x must be provided in the order expected by a PEtabODEProblem, see get_x

See also: PEtabModel, ODESolver, SteadyStateSolver.

Mathemathical description

Following the PEtab standard, the objective function for a PEtabODEProblem is a likelihood (or posterior, when priors are included). For numerical stability, PEtabODEProblem works with the negative log-likelihood; so nllh equals:

$$-\ell (\mathbf{x})=-\sum _{i=1}^N{\ell }_i(\mathbf{x}),$$

where ${\ell }_i$ is the log-likelihood contribution from measurement i$.The‘gradien{t}_{m}ethod‘and‘hessia{n}_{m}ethod‘determinehow$-\nabla \ell(\mathbf{x})$and$-\nabla^2 \ell(\mathbf{x})`` are computed.

In addition to $-\ell$ and its derivatives, PEtabODEProblem provides diagnostics and quantities useful for model assessment. The chi-squared value is computed (per measurement) as [1]:

$${\chi }^2=\frac{(y-h{)}^2}{{\sigma }^2},$$

where y is the measurement, h is the observable value (obs_formula), and σ is the noise standard deviation (noise_formula). Residuals are computed as

$$r=\frac{(y-h)}{\sigma }.$$

The empirical Fisher Information Matrix (FIM) can be used for identifiability analysis. When feasible, it should be computed with an exact Hessian method. The inverse of the FIM provides a lower bound on the parameter covariance matrix. In practice, profile-likelihood methods often provide more reliable results [2].

1. Cedersund et al., The FEBS Journal, pp. 903–922 (2009).

2. Raue et al., Bioinformatics, pp. 1923–1929 (2009).

source

A detailed overview of problem size and configuration is available via:

DataAPI.describe Method

describe(prob::PEtabODEProblem)

Print summary and configuration statistics for prob

source

A PEtabODEProblem has many configurable options. Two important options are the ODESolver and, for problems with steady-state simulations, the SteadyStateSolver:

PEtab.ODESolver Type

ODESolver(solver; kwargs...)

ODE solver configuration (solver, tolerances, etc.) used to simulate the model in a PEtabODEProblem.

Any ODE solver from OrdinaryDiffEq.jl or Sundials.jl is supported. Solver recommendations and the default configuration (when ODESolver is not provided) are described in the online documentation.

Keyword Arguments

• abstol = 1e-8: Absolute tolerance for the forward ODE solve. For gradient-based estimation, values larger than 1e-6 are generally not recommended, as they may yield inaccurate gradients.

• reltol = 1e-8: Relative tolerance for the forward ODE solve. For gradient-based estimation, values larger than 1e-6 are generally not recommended, as they may yield inaccurate gradients.

• dtmin = nothing: Minimum step size for the forward ODE solve.

• maxiters = 10_000: Maximum number of solver iterations. Increasing this can substantially increase runtime during parameter estimation.

• verbose::Bool = true: Whether to print warnings if the solver terminates early. Keeping this enabled is recommended to detect problematic solver configurations.

• adj_solver = solver: ODE solver used for the adjoint solve. Defaults to solver. Only relevant when gradient_method = :Adjoint.

• abstol_adj = abstol: Absolute tolerance for the adjoint solve (only relevant when gradient_method = :Adjoint).

• reltol_adj = reltol: Relative tolerance for the adjoint solve (only relevant when gradient_method = :Adjoint).

source

PEtab.SteadyStateSolver Type

SteadyStateSolver(method::Symbol; kwargs...)

Options for computing a steady state, i.e. a state u where the ODE right-hand side f(u, p, t) satisfies du = f(u, p, t) ≈ 0.

Two approaches are supported:

• method = :Simulate: Simulate the ODE forward in time until a steady-state termination criterion is met. The simulation uses the ODESolver settings provided to the PEtabODEProblem. This is the recommended and most robust approach.

• method = :Rootfinding: Solve f(u, p, t) = 0 directly using a root-finding method from NonlinearSolve.jl. This can be faster, but is generally less reliable (see below).

Keyword Arguments

• termination_check = :wrms: Termination criterion used for method = :Simulate. Options:

  • :wrms: Weighted root-mean-square of du. Terminate when $\sqrt{\frac{1}{N}\sum _{i=1}^N{\left(\frac{d{u}_i}{\mathrm{reltol}{u}_i+\mathrm{abstol}}\right)}^2}<1$.

  • :Newton: Terminate when the Newton step Δu is sufficiently small: $\sqrt{\frac{1}{N}\sum _{i=1}^N{\left(\frac{\mathrm{\Delta }{u}_i}{\mathrm{reltol}{u}_i+\mathrm{abstol}}\right)}^2}<1$. This requires solving a linear system involving the Jacobian of f. If the Jacobian is singular, a pseudo-inverse is used when pseudoinverse = true; otherwise the criterion falls back to :wrms. :wrms is recommended.

• pseudoinverse = true: Use a pseudo-inverse when the Jacobian inversion fails for termination_check = :Newton.

• rootfinding_alg = nothing: NonlinearSolve.jl algorithm used for method = :Rootfinding. If nothing, the NonlinearSolve.jl default is used.

• abstol: Absolute tolerance used both in the termination criterion and (for method = :Rootfinding) the NonlinearSolve.jl solve. Defaults to 100 * abstol_ode, where abstol_ode is taken from the ODESolver in the PEtabODEProblem.

• reltol: Relative tolerance used both in the termination criterion and (for method = :Rootfinding) the NonlinearSolve.jl solve. Defaults to 100 * reltol_ode, where reltol_ode is taken from the ODESolver in the PEtabODEProblem.

• maxiters: Maximum number of iterations for method = :Rootfinding.

Mathemathical description

For an ODE model

$$\frac{\mathrm{d}u}{\mathrm{d}t}=f(u,p,t).$$

A steady state is a solution $u^{\ast }$ such that

$$f(u^{\ast },p,t)=0.$$

The steady state can be obtained either by simulation (:Simulate) or by direct root-finding (:Rootfinding). Although root-finding can be more computationally efficient, it may converge to an unintended root. For example, mass-action models with positive initial conditions often have a physically meaningful steady state with u ≥ 0. Simulation tends to preserve such invariants (except in case of solver failure), while a generic root-finder has no such guarantees and may return a different root that still satisfies f(u, p, t) = 0. In line with this, benchmarks show that :Simulate is more robust [1]. Lastly, if the steady state can be solved for symbolically, it is the best approach [1].

1. Fiedler et al., BMC Systems Biology (2016), pp. 1–19.

source

Utility functions for interacting with a PEtabODEProblem (including efficient subsetting via remake) are:

PEtab.get_x Function

get_x(prob::PEtabODEProblem; linear_scale = false)::ComponentArray

Get the nominal parameter vector with parameters in the correct order expected by prob for parameter estimation/inference. Nominal values can optionally be specified when creating a PEtabParameter, or in the parameters table if the problem is provided in the PEtab standard format.

For ease of interaction (e.g., changing values), the parameter vector is returned as a ComponentArray. For how to interact with a ComponentArray, see the documentation and the ComponentArrays.jl documentation.

See also PEtabParameter.

Keyword argument

• linear_scale: Whether to return parameters on the linear scale. By default, parameters are returned on the scale they are estimated, which by default is log10 as this often improves parameter estimation performance.

source

SciMLBase.remake Method

remake(prob::PEtabODEProblem; conditions=Symbol[], parameters=nothing) -> PEtabODEProblem

Create a new PEtabODEProblem from prob, restricted to a subset of simulation conditions and/or with a subset of parameters to estimate in prob fixed to constant values.

Intended for efficient subsetting (e.g. evaluating nllh/grad!/hess! on a subset of conditions, or with a reduced set of parameters to estimate). Typically faster than constructing a new PEtabODEProblem, since compiled functions from prob are reused (avoids recompilation).

Keyword arguments

• conditions: Simulation conditions to keep. If empty (default), all conditions are kept. Format depends on whether the model has pre-equilibration:

  • No pre-equilibration: Provide Vector{Symbol} of simulation condition ids (e.g. [:cond1, :cond2]).

  • With pre-equilibration: Provide Vector{Pair} of pre_eq_id => simulation_id (e.g. [:pre1 => :cond1, :pre1 => :cond2]).

• experiments: Experimental time course ids to keep, as Vector{Symbol}. Only applicable for problems in PEtab v2 standard format.

• parameters: Parameters to fix to constant values, as a vector of pairs [:p1 => val1, :p2 => val2, ...]. Only parameters that are estimated in prob can be fixed. Values are given on the linear scale; e.g. if a parameter is estimated on :log10, pass val (not log10(val)).

Examples

# Keep only simulation conditions :cond1 and :cond3
prob_sub = remake(prob; conditions = [:cond1, :cond3])

# Fix parameters k1 and k2
prob_sub = remake(prob; parameters = [:k1 => 3.0, :k2 => 4.0])

source

Utilities for interacting with the underlying dynamic model (ODEProblem) are:

PEtab.get_ps Function

get_ps(res, prob::PEtabODEProblem; kwargs...)

Return parameter values for simulating the ODE model in prob.

res can be a parameter-estimation result (e.g. PEtabMultistartResult) or a parameter vector in the internal order expected by prob (see get_x).

Keyword arguments

• condition: Simulation condition to retrieve parameters for. Used for PEtab v1 problems and problems defined via the Julia interface. If the model has pre-equilibration, pass pre_eq_id => simulation_id; otherwise pass simulation_id. IDs may be String or Symbol.

• experiment: Experiment time course to retrieve parameters for (String or Symbol). Only applicable for problem in the PEtab v2 standard format.

• retmap = true: If true, return a vector of pairs p = [k1 => v1, ...] suitable for ODEProblem(; p). If false, return a parameter vector.

Note

condition and experiment are mutually exclusive (provide at most one).

See also: get_u0, get_odeproblem, get_odesol.

source

PEtab.get_u0 Function

get_u0(res, prob::PEtabODEProblem; kwargs...)

Retrieve the ODEProblem initial values for simulating the ODE model in prob. res can be a parameter estimation result (e.g., PEtabMultistartResult) or a Vector with parameters in the order expected by prob (see get_x).

For information on keyword arguments see get_ps.

See also get_odeproblem and get_odesol.

source

PEtab.get_system Function

get_system(res, prob::PEtabODEProblem; kwargs...) -> (sys, p, u0, callbacks)

Retrieve the dynamic model system, parameter map (p), initial species map (u0), and callbacks (CallbackSet) for the model in prob. The argument res can be a parameter estimation result (e.g., PEtabMultistartResult) or a Vector of parameters in the order expected by prob (see get_x).

The system type returned depends on the input to PEtabModel. If the model is provided as a ReactionSystem, a ReactionSystem is returned. The same applies for an ODESystem. If the model is provided via an SBML file, a ReactionSystem is returned.

For information on keyword arguments, see get_ps.

See also: get_u0 and get_odesol.

source

PEtab.get_odeproblem Function

get_odeproblem(res, prob::PEtabODEProblem; kwargs...) -> (ode_prob, cbs)

Retrieve the ODEProblem and cbs (CallbackSet) for simulating the ODE model in prob. res can be a parameter estimation result (e.g., PEtabMultistartResult) or a Vector with parameters in the order expected by prob (see get_x).

For information on keyword arguments see get_ps.

See also: get_u0 and get_odesol.

source

PEtab.get_odesol Function

get_odesol(res, prob::PEtabODEProblem; kwargs...)::ODESolution

Retrieve the ODESolution from simulating the ODE model in prob. res can be a parameter estimation result (e.g., PEtabMultistartResult) or a Vector with parameters in the order expected by prob (see get_x).

For information on keyword arguments see get_ps.

See also: get_u0 and get_odeproblem.

source

PEtab.solve_all_conditions Function

solve_all_conditions(x, prob::PEtabODEProblem, solver; kwargs)

Solve the ODE model in prob for all simulation conditions with the provided ODE-solver.

x should be a Vector or ComponentArray with parameters in the order expected by prob (see get_x).

Keyword Arguments

• abstol=1e-8: Absolute tolerance for the ODE solver.

• reltol=1e-8: Relative tolerance for the ODE solver.

• maxiters=1e4: Maximum iterations for the ODE solver.

• ntimepoints_save=0: The number of time points at which to save the ODE solution for each condition. A value of 0 means the solution is saved at the solvers default time points.

• save_observed_t=false: When set to true, this option overrides ntimepoints_save and saves the ODE solution only at the time points where measurement data is available.

Returns

• odesols: A dictionary containing the ODESolution for each simulation condition.

source

Parameter estimation

A PEtabODEProblem contains everything needed to use an optimizer directly, but PEtab.jl also provides convenience wrappers. For single-start parameter estimation:

PEtab.calibrate Function

calibrate(prob::PEtabODEProblem, x0, alg; kwargs...)::PEtabOptimisationResult

From starting point x0 using optimization algorithm alg, estimate unknown model parameters for prob, and get results as a PEtabOptimisationResult.

x0 can be a Vector or a ComponentArray, where the individual parameters must be in the order expected by prob. To get a vector in the correct order, see get_x.

A list of available and recommended optimization algorithms (alg) can be found in the documentation. Briefly, supported algorithms are from:

• Optim.jl: LBFGS(), BFGS(), or IPNewton() methods.

• Ipopt.jl: IpoptOptimizer() interior-point Newton method.

• Fides.jl: Newton trust region for box-constrained problem. Fides can either use the Hessian in prob (alg = Fides.CustomHessian()), or any of its built in Hessian approximations (e.g. alg = Fides.BFGS()). A full list of Hessian approximations can be found in the Fides documentation

Different ways to visualize the parameter estimation result can be found in the documentation.

See also PEtabOptimisationResult and IpoptOptimizer

Keyword Arguments

• save_trace::Bool = false: Whether to save the optimization trace of the objective function and parameter vector. Only applicable for some algorithms; see the documentation for details.

• options = DEFAULT_OPTIONS: Configurable options for alg. The type and available options depend on which package alg belongs to. For example, if alg = IPNewton() from Optim.jl, options should be provided as an Optim.Options() struct. A list of configurable options can be found in the documentation.

source

PEtab.PEtabOptimisationResult Type

PEtabOptimisationResult

Parameter estimation statistics from single-start optimization with calibrate.

See also: calibrate

Fields

• xmin: Minimizing parameter vector found by the optimization.

• fmin: Minimum objective value found by the optimization.

• x0: Starting parameter vector.

• alg: Parameter estimation algorithm used.

• niterations: Number of iterations for the optimization.

• runtime: Runtime in seconds for the optimization.

• xtrace: Parameter vector optimization trace. Empty if save_trace = false was provided to calibrate.

• ftrace: Objective function optimization trace. Empty if save_trace = false was provided to calibrate.

• converged: Convergence flag from alg.

• original: Original result struct returned by alg. For example, if alg = IPNewton() from Optim.jl, original is the Optim return struct.

source

For multi-start parameter estimation:

PEtab.calibrate_multistart Function

calibrate_multistart([rng::AbstractRng], prob::PEtabODEProblem, alg, nmultistarts::Integer;
                     nprocs = 1, dirsave=nothing, kwargs...)::PEtabMultistartResult

Perform nmultistarts parameter estimation runs from randomly sampled starting points using the optimization algorithm alg to estimate the unknown model parameters in prob.

A list of available and recommended optimisation algorithms (alg) can be found in the package documentation and in the calibrate documentation.

As with get_startguesses, the rng controlling the generation of starting points is optional; if omitted, Random.default_rng() is used. For reproducible starting points, pass a seeded rng (e.g., MersenneTwister(42)).

If nprocs > 1, the parameter estimation runs are performed in parallel using the pmap function from Distributed.jl with nprocs processes. If parameter estimation on a single process (nprocs = 1) takes longer than 5 minutes, we strongly recommend setting nprocs > 1, as this can greatly reduce runtime. Note that nprocs should not be larger than the number of cores on the computer.

If dirsave is provided, intermediate results for each run are saved in dirsave. It is strongly recommended to provide dirsave for larger models, as parameter estimation can take hours (or even days!),and without dirsave, all intermediate results will be lost if something goes wrong.

Different ways to visualize the parameter estimation result can be found in the documentation.

See also PEtabMultistartResult, get_startguesses, and calibrate.

Keyword Arguments

• sampling_method = LatinHypercubeSample(): Method for sampling a diverse (spread out) set of starting points. See the documentation for get_startguesses, which is the function used for generating starting points.

• sample_prior::Bool = true: See the documentation for get_startguesses.

• options = DEFAULT_OPTIONS: Configurable options for alg. See the documentation for calibrate.

source

PEtab.get_startguesses Function

get_startguesses([rng::AbstractRNG], prob::PEtabODEProblem, n::Integer; kwargs...)

Generate n random parameter vectors within the parameter bounds in prob.

rng is optional and if omitted defaults to Random.default_rng(). If n = 1, a single random vector is returned. For n > 1, a vector of random parameter vectors is returned. In both cases, parameter vectors are returned as a ComponentArray. For details on how to interact with a ComponentArray, see the documentation and the ComponentArrays.jl documentation.

See also calibrate and calibrate_multistart.

Keyword Arguments

• sampling_method = LatinHypercubeSample(): Method for sampling a diverse (spread out) set of parameter vectors. Any algorithm from QuasiMonteCarlo is allowed, but the default LatinHypercubeSample is recommended as it usually performs well.

• sample_prior::Bool = true: Whether to sample random parameter values from the prior distribution if a parameter has a prior.

• allow_inf::Bool = false: Whether to return parameter vectors for which the likelihood cannot be computed (typically happens because the ODEProblem cannot be solved). Often it only makes sense to use starting points with a computable likelihood for parameter estimation, hence it typically does not make sense to change this option.

source

PEtab.PEtabMultistartResult Type

PEtabMultistartResult

Parameter estimation statistics from multi-start optimization with calibrate_multistart.

See also calibrate_multistart and PEtabOptimisationResult.

Fields

• xmin: Best minimizer across all runs.

• fmin: Best minimum across all runs.

• alg: Parameter estimation algorithm.

• nmultistarts: Number of parameter estimation runs.

• sampling_method: Sampling method used for generating starting points.

• dirsave: Path of directory where parameter estimation run statistics are saved if dirsave was provided to calibrate_multistart.

• runs: Vector of PEtabOptimisationResult with the parameter estimation results for each run.

PEtabMultistartResult(dirres::String; which_run::String="1")

Import multistart parameter estimation results saved at dirres.

Each time a new optimization run is performed, results are saved with unique numerical endings. Results from a specific run can be retrieved by specifying the numerical ending with which_run.

source

And finally model selection (PEtab-Select interface):

PEtab.petab_select Function

petab_select(path_yaml, alg; nmultistarts = 100, kwargs...) -> path_res

For a PEtab-select problem, perform model selection with the method specified in the PEtab select problem files. Returns the path (path_res) to a YAML-file with model selection results.

The general model selection (e.g. to use forward-search) options are specified in the PEtab-select files. For details on how to set this up, see the PEtab-select documentation.

For each round of model selection, the candidate models are parameter estimated using multi-start parameter estimation, with nmultistarts performed for each model. The objective values obtained from parameter estimation are then used for the next round of model evaluation.

A list of available and recommended optimization algorithms (alg) can be found in the package documentation and calibrate documentation.

See also calibrate_multistart.

Keyword Arguments

• kwargs: The same keywords accepted by PEtabODEProblem and calibrate.

source

For calibration, optimizers from Optim.jl, Ipopt.jl, and Fides.jl. Ipopt-specific configuration is available via:

PEtab.IpoptOptimizer Type

IpoptOptimizer(LBFGS::Bool)

Setup the Ipopt Interior-point Newton method optmizer for parameter estimation.

Ipopt can be configured to use either the Hessian method from the PEtabODEProblem (LBFGS=false) or an LBFGS scheme (LBFGS=true). For setting other Ipopt options, see IpoptOptions.

See also calibrate and calibrate_multistart.

Description

Ipopt is an Interior-point Newton method for constrained non-linear optimization problems. More information on the algorithm can be found in [1].

1. Wächter and Biegler, Mathematical programming, pp 25-57 (2006)

source

PEtab.IpoptOptions Type

IpoptOptions(; kwargs...)

Options for parameter estimation with IpoptOptimizer.

More details on the options can be found in the Ipopt documentation.

See also IpoptOptimizer, calibrate, and calibrate_multistart.

Keyword Arguments

• print_level = 0: Output verbosity level (valid values are 0 ≤ print_level ≤ 12)

• max_iter = 1000: Maximum number of iterations

• tol = 1e-8: Relative convergence tolerance

• acceptable_tol = 1e-6: Acceptable relative convergence tolerance

• max_wall_time 1e20: Maximum wall time optimization is allowed to run

• acceptable_obj_change_tol 1e20: Acceptance stopping criterion based on objective function change.

source

Parameter-estimation results can be visualized using the plotting recipes described on Plotting parameter estimation results, and with:

PEtab.get_obs_comparison_plots Function

get_obs_comparison_plots(res, prob::PEtabODEProblem; kwargs...)::Dict

Plot the model fit against data for all simulation conditions and all observable ids for a PEtab.jl parameter estimation result (res).

Each entry in the returned Dict corresponds to plot(res, prob; obsids=[obsid], cid=cid, kwargs...) for all possible cid and all obsid.

source

As an alternative workflow, a PEtabODEProblem can be converted to an OptimizationProblem to access solvers via Optimization.jl:

PEtab.OptimizationProblem Function

OptimizationProblem(prob::PEtabODEProblem; box_constraints::Bool = true)

Create an Optimization.jl OptimizationProblem from prob.

To use algorithms not compatible with box constraints (e.g., Optim.jl NewtonTrustRegion), set box_constraints = false. Note that with this option, optimizers may move outside the bounds, which can negatively impact performance. More information on how to use an OptimizationProblem can be found in the Optimization.jl documentation.

source

Parameter-estimation results for problems in the PEtab standard format can be exported with:

PEtab.export_petab Function

export_petab(dir, prob::PEtabODEProblem, res) -> yaml_path

Export in the PEtab standard format prob to the directory dir, using the parameters in res to populate the PEtab parameter table. res may be a parameter-estimation result (e.g. PEtabMultistartResult) or a parameter vector in the order expected by prob (see get_x).

Returns the path to the exported PEtab problem YAML file.

The exporter currently supports only problems that were provided in the PEtab standard format (PEtab tables), and exported tables keep the same filenames as in the original PEtab problem. Problems constructed via the Julia interface are not yet exportable.

source

Bayesian inference

PEtab.jl supports Bayesian inference by exposing a PEtabLogDensity compatible with the LogDensityProblems.jl interface, which can be sampled with AdvancedHMC.jl (including NUTS) or AdaptiveMCMC.jl. This functionality is planned to move into a separate package, so the API will change in the future.

PEtab.PEtabLogDensity Type

PEtabLogDensity(prob::PEtabODEProblem)

Create a LogDensityProblem using the posterior and gradient functions from prob.

This LogDensityProblem interface defines everything needed to perform Bayesian inference with packages such as AdvancedHMC.jl (which includes algorithms like NUTS, used by Turing.jl), and AdaptiveMCMC.jl.

source

PEtab.to_prior_scale Function

to_prior_scale(xpetab, target::PEtabLogDensity)

Transforms parameter x from the PEtab problem scale to the prior scale.

This conversion is needed for Bayesian inference, as in PEtab.jl Bayesian inference is performed on the prior scale.

Note

To use this function, the Bijectors, LogDensityProblems, and LogDensityProblemsAD packages must be loaded: using Bijectors, LogDensityProblems, LogDensityProblemsAD

source

PEtab.to_chains Function

to_chains(res, target::PEtabLogDensity; kwargs...)::MCMCChains

Converts Bayesian inference results obtained with PEtabLogDensity into an MCMCChains.

res can be the inference results from AdvancedHMC.jl or AdaptiveMCMC.jl. The returned chain has the parameters on the prior scale.

Keyword Arguments

• start_time: Optional starting time for the inference, obtained with now().

• end_time: Optional ending time for the inference, obtained with now().

Note

To use this function, the MCMCChains package must be loaded: using MCMCChains

source

In addition to priors from Distributions.jl, PEtab.jl also supports a LogLaplace prior:

PEtab.LogLaplace Type

LogLaplace(μ,θ)

The log-Laplace distribution with location μ and scale θ has probability density function

$$f(x;\mu ,\theta )=\frac{1}{2\theta x}\exp (-\frac{|lnx-\mu |}{\theta })$$

External links

• Log-Laplace distribution on Wikipedia

Implementation note

LogLaplace does not yet have support in Distributions.jl, so to support it for Bayesian inference PEtab.jl implements support its pdf, logpdf cdf, logcdf, median and sampling

source