From 3a91939fe12efe5c0ef0aab98508e691c1679b6d Mon Sep 17 00:00:00 2001 From: Torkel Date: Fri, 6 Oct 2023 21:47:11 -0400 Subject: [PATCH 01/32] Init --- docs/pages.jl | 17 +- docs/src/assets/petab_best_objective.svg | 88 ++++ docs/src/assets/petab_waterfall.svg | 88 ++++ .../parameter_estimation.md | 0 .../petab_ode_param_fitting.md | 487 ++++++++++++++++++ 5 files changed, 672 insertions(+), 8 deletions(-) create mode 100644 docs/src/assets/petab_best_objective.svg create mode 100644 docs/src/assets/petab_waterfall.svg rename docs/src/{catalyst_applications => inverse_problems}/parameter_estimation.md (100%) create mode 100644 docs/src/inverse_problems/petab_ode_param_fitting.md diff --git a/docs/pages.jl b/docs/pages.jl index ca2241c5b7..2b67fd0ba5 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -1,19 +1,20 @@ pages = Any["Home" => "index.md", "Introduction to Catalyst" => Any["introduction_to_catalyst/catalyst_for_new_julia_users.md", - "introduction_to_catalyst/introduction_to_catalyst.md"], + "introduction_to_catalyst/introduction_to_catalyst.md" ], "Catalyst Functionality" => Any["catalyst_functionality/dsl_description.md", "catalyst_functionality/programmatic_CRN_construction.md", "catalyst_functionality/compositional_modeling.md", "catalyst_functionality/constraint_equations.md", "catalyst_functionality/parametric_stoichiometry.md", - "catalyst_functionality/network_analysis.md"], + "catalyst_functionality/network_analysis.md", + "Model creation examples" => Any["example_networks/basic_CRN_examples.md", + "example_networks/hodgkin_huxley_equation.md", + "example_networks/smoluchowski_coagulation_equation.md"]], "Catalyst Applications" => Any["catalyst_applications/simulation_structure_interfacing.md", "catalyst_applications/advanced_simulations.md", "catalyst_applications/homotopy_continuation.md", - "catalyst_applications/bifurcation_diagrams.md", - "catalyst_applications/parameter_estimation.md"], - "Example Networks" => Any["example_networks/basic_CRN_examples.md", - "example_networks/hodgkin_huxley_equation.md", - "example_networks/smoluchowski_coagulation_equation.md"], + "catalyst_applications/bifurcation_diagrams.md"], + "Inverse Problems" => Any["inverse_problems/parameter_estimation.md", + "inverse_problems/petab_ode_param_fitting.md"], "FAQs" => "faqs.md", - "API" => "api/catalyst_api.md"] + "API" => "api/catalyst_api.md"] \ No newline at end of file diff --git a/docs/src/assets/petab_best_objective.svg b/docs/src/assets/petab_best_objective.svg new file mode 100644 index 0000000000..1e09fed0fa --- /dev/null +++ b/docs/src/assets/petab_best_objective.svg @@ -0,0 +1,88 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/assets/petab_waterfall.svg b/docs/src/assets/petab_waterfall.svg new file mode 100644 index 0000000000..03e603bb45 --- /dev/null +++ b/docs/src/assets/petab_waterfall.svg @@ -0,0 +1,88 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/catalyst_applications/parameter_estimation.md b/docs/src/inverse_problems/parameter_estimation.md similarity index 100% rename from docs/src/catalyst_applications/parameter_estimation.md rename to docs/src/inverse_problems/parameter_estimation.md diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md new file mode 100644 index 0000000000..44d7a7b88c --- /dev/null +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -0,0 +1,487 @@ +# [Parameter Fitting for ODEs using PEtab.jl](@id petab_parameter_fitting) +The [PEtab.jl package](https://github.com/sebapersson/PEtab.jl) implements the [PEtab format](https://petab.readthedocs.io/en/latest/) for fitting the parameters of deterministic CRN models to data ([please cite the corresponding papers if you use the PEtab approach in your research](@ref petab_citations)). PEtab.jl both implements methods for creating cost functions (determining how well parameter sets fit to data), and for minimizing these cost functions. The PEtab approach covers most cases of fitting deterministic (ODE) models to data. When possible, it is recommended to use it for parameter fitting (when not, an alternative, more general, approach will have to be used). This page describes how to combine PEtab.jl and Catalyst for parameter fitting, with the PEtab.jl package providing [a more extensive documentation](https://sebapersson.github.io/PEtab.jl/stable/) (this tutorial is partially an adaptation of this documentation). + + +## Introductory example +Let us consider a simple catalysis network, where an enzyme (*E*) turns a substrate (*S*) into a product (*P*): +```petab1 +using Catalyst, PEtab + +rn = @reaction_network begin + kB, S + E --> SE + kD, SE --> S + E + kP, SE --> P + E +end +``` +From some known initial condition, and a true parameter set (which we later want to recover from the data) we generate synthetic data (on which we will demonstrate the fitting process). +```petab1 +# Define initial conditions and parameters. +u0 = [:S => 1.0, :E => 1.0, :SE => 0.0, :P => 0.0] +p_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5] + +# Generate synthetic data. +using OrdinaryDiffEq, Random +oprob = ODEProblem(rn, u0, (0.0, 10.0), p_true) +sol = solve(oprob, Tsit5(); saveat=0.1) +data = (0.8 .+ 0.4*rand(10)) .* sol[:P][10:10:end] + +# Plots the true solutions and the (synthetic data) measurements. +using Plots +plot(sol; idxs=:P, label="True solution") +plot!(sol.t[10:10:end], data; label="Measurements", seriestype=:scatter) +``` + +Generally, PEtab takes five different inputs to define an optimisation problem (the minimiser of which generates a fitted parameter set): +1. **Model**: The model which we want to fit to the data (a `ReactionSystem`). +2. **Observables**: The possible observables that can be measured (a `Dict{String,PEtabObservable}`). +3. **Estimation parameters**: The parameters which we want to fit to the data (a `Vector{PEtabParameter}`). +4. **Experimental (or simulation) conditions**: The simulations (each corresponding to a potential experiment) carried out during each step of the optimisation process (a `Dict{String,Dict{Symbol,Float64}}`). +5. **Measurements**: The measurements to which the model is fitted (a `DataFrame`). + +### Observables +The observables define the quantities that we may measure in our experiments. Typically, each corresponds to a single species, however, [more complicated observables are possible](@ref petab_observables_observables). For each observable, we also need a noise formula, defining the uncertainty in its measurements. By default, PEtab assumes normally distributed noise, with a mean equal to the true value and a standard deviation which we have to define. It is also possible to use [more advanced noise formulas](@ref petab_observables_noise_formula). + +In our example, we only have a single possible observable, the `P` species. We will assume that the noise is normally distributed with a standard deviation `0.5` (in our case this is not true, however, typically the noise distribution is unknown and a guess have to be made). We combine this information in a `PEtabObservable` struct (to access the `P` species we must use [`@unpack`](@ref simulation_structure_interfacing_symbolic_representation)). Finally, we store all our observables in a dictionary, giving each an id tag (which is later used in the measurements input). + +```petab1 +@unpack P = rn +obs_P = PEtabObservable(P, 0.5) +observables = Dict("obs_P" => obs_P) +nothing # hide +``` + +### Parameters +Each parameter of the system can either be +1. Known (described [here](@ref petab_parameters_known)). +2. Depend on experimental/simulation conditions (described [here](@ref petab_simulation_conditions)). +3. Be an unknown that we wish to fit to data. + +In our case, we wish to fit all three system parameters (*kB*, *kD*, and *kP*). For each, we create a single `PEtabParameter`, and then gather these into a single vector. +```petab1 +par_kB = PEtabParameter(:kB) +par_kD = PEtabParameter(:kD) +par_kP = PEtabParameter(:kP) +params = [par_kB, par_kD, par_kP] +nothing # hide +``` +For each parameter, it is also possible to set [a lower and/or upper bound](@ref petab_parameters_bounds), set whether to use [logarithmic or linear scale](@ref petab_parameters_scales), or add a [prior distribution of its value](@ref petab_parameters_priors). + +### Simulation conditions +Sometimes, several different experiments are performed on a system (each potentially generating several measurements). An experiment could e.g. be the time development of a system from a specific initial condition. Since each experimental condition (during the optimisation procedure, for a guess of the unknown parameters) generates a distinct simulation, these are also called simulation conditions. In our example, all data comes from a single experiment, and the simulation condition input is not required. How to define and use different experimental conditions is described [here](@ref petab_simulation_conditions). + + +### Measurements +Finally, we need to define the system measurements to which the parameters will be fitted. Each measurement combines: +1. The observable which is observed (here we use the id tag defined in the `observables` dictionary). +2. The time point of the measurement. +3. The measured value. + +For cases where several simulation conditions are given, we also need to provide: + +4. The simulation condition which generates the measurement ([here](@ref petab_simulation_conditions) is an example where this is used). + +([When [pre-equilibration](https://sebapersson.github.io/PEtab.jl/stable/Brannmark/) is used to initiate the system in a steady state, a fifth field is also required]) + +Measurements are provided in a `DataFrame` where each row corresponds to a measurement. The values are provided in the `obs_id`, `time`, `measurement`, and `simulation_id` columns. In our case we only need to fill in the first three: +```petab1 +using DataFrames +measurements = DataFrame(obs_id="obs_P", time=sol.t[10:10:end], measurement=data) +``` +Since, in our example, all measurements are of the same observable, we can set `obs_id="obs_P"`. If several different observables were measured, `obs_id` would be a vector of the same length as `time` and `measurement`. + +### Creating a PEtabModel +Finally, we combine all inputs in a single `PEtabModel`. To it, we also pass the initial conditions of our simulations (using the `state_map` argument). It is also possible to have [initial conditions with uncertainty](@ref petab_simulation_initial_conditions_uncertainty), [that vary between different simulations](@ref petab_simulation_conditions), or [that we attempt to fit to the data](@ref petab_simulation_initial_conditions_fitted). +```petab1 +petab_model = PEtabModel(rn, observables, measurements, params; state_map=u0) +nothing # hide +``` + +### Fitting parameters +We are now able to fit our model to the data. First, we create a `PEtabODEProblem`. Here, we use `petab_model` as the only input, but it is also possible to set various [numeric solver and automatic differentiation option](@ref petab_simulation_conditions) (such as method or tolerance). +```petab1 +petab_problem = PEtabODEProblem(petab_model) +nothing # hide +``` +To fit a parameter set we use the `calibrate_model` function. In addition to our `PEtabODEProblem`, we must also provide an initial guess (which can be generated with the `generate_startguesses` function) and an optimisation algorithm (which needs to be imported specifically). PEtab.jl supports various [optimisation methods and options](@ref petab_optimisation_optimisers). +```petab1 +using Optim +p0 = generate_startguesses(petab_problem, 1) +res = calibrate_model(petab_problem, p0, IPNewton()) +``` + +We can now simulate our model for the fitted parameter set, and compare the result to the measurements and true solution. +```petab1 +oprob = ODEProblem(rn, u0, (0.0, 10.0), get_ps(res, petab_problem)) + +sol = solve(oprob, Tsit5(); saveat=0.1) +plot!(sol; idxs=4, label="Fitted solution") +``` +Here we use the `get_ps` function to retrieve a full parameter set using the optimal parameters (alternatively, the `ODEProblem` or fitted simulation can be retrieved directly using the `get_odeproblem` or `get_odesol` [functions](https://sebapersson.github.io/PEtab.jl/dev/API_choosen/#PEtab.get_odeproblem), respectively). The calibration result can also be found in `res.xmin`, however, note that PEtab automatically ([unless a linear scale is selected](petab_parameters_scales)) converts parameters to logarithmic scale, so typically `10 .^res.xmin` are the values of interest. If you investigate the result from this example you might note that the found parameter set (while it fits the data well) does not correspond to the true parameter set. This phenomenon is related to *identifiability*, and is very important for parameter fitting. + +### Final notes +PEtab.jl also supports [multistart optimisation](@ref petab_multistart_optimisation), [automatic equilibration before simulations](https://sebapersson.github.io/PEtab.jl/stable/Brannmark/), and [events](@ref petab_events). Various [plot recipes](@ref petab_plotting) exist for investigating the optimisation process. Please read the [PETab.jl documentation](https://sebapersson.github.io/PEtab.jl/stable/) for a more complete description of the package's features. Below follows additional details of various options and features (generally, PEtab is able to find good default values for most options not specified). + +## [Additional features: Observables](@id petab_observables) + +### [Defining non-trivial observables](@id petab_observables_observables) +It is possible for observables to be any algebraic expression of species concentrations and parameters. E.g. in this example the total amount of `X` in the system is an observable: +```petab2 +using Catalyst, PEtab # hide +rn = @reaction_network begin + (k1,k2), X1 <--> X2 +end +@unpack X1, X2 = rn +obs_X = PEtabObservable(X1 + X2, 0.5) +``` + +### [Advanced observables noise formulas](@id petab_observables_noise_formula) +In our basic example we assumed that the normally distributed noise had a standard deviation of `0.5`. However, this value may be a parameter (or indeed any algebraic expression). E.g, we could set +```petab1 +@parameters σ +obs_P = PEtabObservable(P, σ) +``` +and then let the parameter *σ* vary between different [simulation conditions](@ref petab_simulation_conditions). Alternatively we could let the noise scale linearly with the species' concentration: +```petab1 +obs_P = PEtabObservable(P, 0.05P) +``` + +PEtab.jl assumes that noise is normally distributed (with a standard deviation equal to the second argument provided to `PEtabObservable`). The only other (currently implemented) noise distribution is log-normally distributed noise, which is designated through the `transformation=:log` argument: +```petab1 +obs_P = PEtabObservable(P, σ; transformation=:log) +``` + +## [Additional features: Parameters](@id petab_parameters) + +### [Known parameters](@id petab_parameters_known) +In our previous example, all parameters were unknowns that we wished to fit to the data. If any parameters have known values, it is possible to provide these to `PEtabModel` through the `parameter_map` argument. E.g if we had known that *kB = 1.0*, then we would only define *kD* and *kP* as parameters we wish to fit: +```petab1 +par_kD = PEtabParameter(:kD) +par_kP = PEtabParameter(:kP) +params = [par_kD, par_kP] +nothing # hide +``` +We then provide `parameter_map=[:kB => 1.0]` when we assembly our model: +```petab1 +petab_model = PEtabModel(rn, observables, measurements, params; state_map=u0, parameter_map=[:kB => 1.0]) +nothing # hide +``` + +### [Parameter bounds](@id petab_parameters_bounds) +By default, when fitted, potential parameter values are assumed to be in the interval *[1e-3, 1e3]*. When declaring a `PEtabParameter` it is possible to change these values through the `lb` and `ub` arguments. E.g. we could use +```petab1 +par_kB = PEtabParameter(:kB; lb=1e-2, ub=1e-2) +``` +to achieve the more conservative bound *[1e-2, 1e2]* for the parameter *kB*. + +### [Parameter scales](@id petab_parameters_scales) + +Internally, parameters that are fitted are converted to a logarithmic scale (generally, this is considered advantageous). To prevent this conversion, the `scale=:lin` argument can be used. Here we change the scale of *kB* to linear: + +```petab1 +par_kB = PEtabParameter(:kB; scale=:lin) +``` + +### [Parameter priors](@id petab_parameters_priors) +If we has prior knowledge about the distribution of a parameter, it is possible to incorporate this into the model. The prior can be any continuous, univariate, distribution from the [Distributions.jl package](https://github.com/JuliaStats/Distributions.jl). E.g we can use: + +```petab1 +using Distributions +par_kB = PEtabParameter(:kB; prior=Normal(1.0,0.2)) +``` +to set a normally distributed prior (with mean `1.0` and standard deviation `0.2`) on the value of *kB*. By default, the prior is assumed to be on the linear scale of the parameter (before any potential log transform). To specify that the prior is on the logarithmic scale, the `prior_on_linear_scale=false` argument can be provided: +```petab1 +par_kB = PEtabParameter(:kB; prior=Normal(1.0,0.2), prior_on_linear_scale=false) +``` + +## [Simulation conditions](@id petab_simulation_conditions) +Sometimes, we have data from several different experimental conditions. Here, when a potential parameter set is evaluated during the fitting process, each experimental condition corresponds to one simulation condition (which produces one simulation). To account for this, PEtab permits the user to define several different simulation conditions, with each condition being defined by specific values for some initial conditions and/or parameters. + +If, for our previous catalysis example, we had measured the system for two different initial values of *S* (*S(0)=1.0* and *S(0)=0.5*), these would correspond to two different simulation conditions. For each condition we define a `Dict` mapping the species to their initial condition (here, *S* is the only species in each `Dict`): +```petab1 +c1 = Dict(:S => 1.0) +c2 = Dict(:S => 0.5) +nothing # hide +``` +Similarly as for observables, we then gather the conditions in another `Dict`, giving each an id tag: +```petab1 +simulation_conditions = Dict("c1" => c1, "c2" => c2) +nothing # hide +``` +Again (like for observables), each measurement in the measurements `DataFrame` needs to be associated with a simulation condition id tag (describing which condition that measurements were taken from). Parameters, just like initial conditions, may vary between different conditions. If an initial condition (or parameter) occurs in one condition, it must occur in all of them. + +Here follows a complete version of our basic example, but with measurements both for *S(0)=1.0* and *S(0)=0.5*. +```petab3 +using Catalyst, PEtab + +rn = @reaction_network begin + kB, S + E --> SE + kD, SE --> S + E + kP, SE --> P + E +end + +u0 = [:E => 1.0, :SE => 0.0, :P => 0.0] +p_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5] + +# Simulate data. +using OrdinaryDiffEq, Random +d1, t1 = let + oprob = ODEProblem(rn, [:S=>1.0; u0], (0.0, 10.0), p_true) + sol = solve(oprob, Tsit5(); saveat=0.1) + ((0.8 .+ 0.4*rand(10)) .* sol[:P][10:10:end], sol.t[10:10:end]) +end +d2, t2 = let + oprob = ODEProblem(rn, [:S=>0.5; u0], (0.0, 10.0), p_true) + sol = solve(oprob, Tsit5(); saveat=0.1) + ((0.8 .+ 0.4*rand(10)) .* sol[:P][10:10:end], sol.t[10:10:end]) +end + +@unpack P = rn +obs_P = PEtabObservable(P, 0.5) +observables = Dict("obs_P" => obs_P) + +par_kB = PEtabParameter(:kB) +par_kD = PEtabParameter(:kD) +par_kP = PEtabParameter(:kP) +params = [par_kB, par_kD, par_kP] + +c1 = Dict(:kB => 1.0) +c2 = Dict(:kB => 0.5) +simulation_conditions = Dict("c1" => c1, "c2" => c2) + +using DataFrames +m1 = DataFrame(simulation_id="c1", obs_id="obs_P", time=t1, measurement=d1) +m2 = DataFrame(simulation_id="c2", obs_id="obs_P", time=t2, measurement=d2) +measurements = vcat(m1,m2) + +petab_model = PEtabModel(rn, simulation_conditions, observables, measurements, params; state_map=u0) +nothing # hide +``` +Note that the `u0` we pass into `PEtabModel` through the `state_map` argument no longer contains the value of *S* (as it is provided by the conditions). + +### [Varying parameters between different simulation conditions](@id petab_simulation_conditions) +Sometimes, the parameters that are used vary between the different conditions. Consider our catalysis example, if we had performed the experiment twice, using two different enzymes with different catalytic properties, this could have generated such conditions. The two enzymes could e.g. yield different rates (*kP₁* and *kP₂*) for the `SE --> P + E` reaction, but otherwise be identical. Here, the parameters *kP₁* and *kP₂* are unique to their respective conditions. PEtab.jl provides support for cases such as this, and [its documentation](https://sebapersson.github.io/PEtab.jl/dev/Beer/) provided instructions of how to handle them. + +## [Additional features: Initial conditions](@id petab_simulation_initial_conditions) + +### [Fitting initial conditions](@id petab_simulation_initial_conditions_fitted) +Sometimes, initial conditions are uncertain quantities which we wish to fit to the data. This is possible by [defining an initial condition as a parameter](dsl_description_parametric_initial_conditions): +```petab4 +using Catalyst, PEtab # hide +rn = @reaction_network begin + @parameters E0 + @species E(t)=E0 + kB, S + E --> SE + kD, SE --> S + E + kP, SE --> P + E +end +nothing # hide +``` +Here, the initial value of `E` is equal to the parameter `E0`. We modify our `u0` vector by removing `E` (which is no longer known): +```petab4 +u0 = [:S => 1.0, :SE => 0.0, :P => 0.0] +nothing # hide +``` +Next, we add `E0` to the parameters we wish to fit: +```petab4 +par_kB = PEtabParameter(:kB) +par_kD = PEtabParameter(:kD) +par_kP = PEtabParameter(:kP) +par_E0 = PEtabParameter(:E0) +params = [par_kB, par_kD, par_kP, par_E0] +nothing # hide +``` +and we can use our updated `rn`, `u0`, and `params` as input to our `PEtabModel`. + +### [Uncertain initial conditions](@id petab_simulation_initial_conditions_uncertainty) +Often, while an initial condition has been reported for an experiment, its exact value is uncertain. This can be modelled by making the initial condition a [parameter that is fitted to the data](@id petab_simulation_initial_conditions_fitted) and attach a prior to it corresponding to our certainty about its value. + +Let us consider our initial example, but where we want to add uncertainty to the initial conditions of `S` and `E`. We will add priors on these, assuming normal distributions with mean `1.0` and standard deviation `0.1`. For the synthetic measured data we will use the true values *S(0)=E(0)=1.0*. +```petab5 + +using Catalyst, PEtab + +rn = @reaction_network begin + @parameters S0 E0 + @species S(t)=S0 E(t)=E0 + kB, S + E --> SE + kD, SE --> S + E + kP, SE --> P + E +end + +u0 = [:SE => 0.0, :P => 0.0] +p_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5, :S0=>1.0, :E0 => 1.0] + +using OrdinaryDiffEq, Random +oprob = ODEProblem(rn, [:S=>1.0; :E => 1.0; u0], (0.0, 10.0), p_true) +sol = solve(oprob, Tsit5(); saveat=0.1) +data = (0.8 .+ 0.4*rand(10)) .* sol[:P][10:10:end] + +@unpack P = rn +obs_P = PEtabObservable(P, 0.5) +observables = Dict("obs_P" => obs_P) + +par_kB = PEtabParameter(:kB) +par_kD = PEtabParameter(:kD) +par_kP = PEtabParameter(:kP) +par_S0 = PEtabParameter(:S0) +par_E0 = PEtabParameter(:E0) +params = [par_kB, par_kD, par_kP, par_S0, par_E0] + +using DataFrames +measurements = DataFrame(obs_id="obs_P", time=sol.t[10:10:end], measurement=data) + +petab_model = PEtabModel(rn, observables, measurements, params; state_map=u0) +nothing # hide +``` +Here, when we fit our data we will also gain values for `S0` and `E0`, however, unless we are explicitly interested in these, they can be ignored. + +## [Additional features: Simulation options](@id petab_simulation_options) + +While in our basic example, we do not provide any additional information to our `PEtabODEProblem`, this is an opportunity to specify how the model should be simulated, and what automatic differentiation techniques to use for the optimisation procedure (if none are provided, appropriate defaults are selected). + +Here is an example, taken from the [more detailed PEtab.jl documentation](https://sebapersson.github.io/PEtab.jl/dev/Boehm/#Creating-a-PEtabODEProblem) +```petab1 +petab_problem = PEtabODEProblem(petab_model, + ode_solver=ODESolver(Rodas5P(), abstol=1e-8, reltol=1e-8), + gradient_method=:ForwardDiff, + hessian_method=:ForwardDiff) +nothing # hide +``` +where we simulate our ODE model using the `Rodas5p` method (with absolute and relative tolerance both equal `1e-8`) and use [forward automatic differentiation](https://github.com/JuliaDiff/ForwardDiff.jl) for both gradient and hessian computation. More details on available ODE solver options can be found in [the PEtab.jl documentation](https://sebapersson.github.io/PEtab.jl/dev/API_choosen/#PEtab.ODESolver). + +## [Additional features: Optimisation](@id petab_optimisation) + +### [Optimisation methods and options](@id petab_optimisation_optimisers) +For our examples, we have used the `Optim.IPNewton` optimisation method. PEtab.jl supports [several additional optimisation methods](https://sebapersson.github.io/PEtab.jl/dev/Avaible_optimisers/#options_optimisers). Furthermore, `calibrate_model`'s `options` argument permits the customisation of the options for any used optimiser. E.g. to designate the maximum number of iterations of the `Optim.IPNewton` method we would use: +```petab1 +res = calibrate_model(petab_problem, p0, IPNewton(); options=Optim.Options(iterations = 10000)) +nothing # hide +``` + +Please read the [PEtab.jl documentation](https://sebapersson.github.io/PEtab.jl/dev/Avaible_optimisers/#options_optimisers) to learn how to customise the various optimisers' properties. + +### [Optimisation path recording](@id petab_optimisation_path_recording) + +To record all the parameter sets evaluated (and their objective values) during the optimisation procedure, provide the `save_trace=true` argument to `calibrate_model` (or `calibrate_model_multistart`). This is required for the various [optimisation evaluation plots](@ref petab_plotting) provided by PEtab.jl. If desired, this information can be accessed in the calibration output's `.xtrace` and `.ftrace` fields. + + +## Objective function extraction +While PEtab.jl provides various tools for analysing the objective function generated by `PEtabODEProblem`, it is also possible to extract this function for customised analysis. Given a `PEtabODEProblem` +```petab1 +petab_problem = PEtabODEProblem(petab_model) +nothing # hide +``` +We can find the: +1. Objective function in the `petab_problem.compute_cost`. It takes a single argument (`p`) and returns the objective value. +2. Gradient in the `petab_problem.compute_gradient!` field. It takes two arguments (`g` and `p`) with the updated gradient values being written to `g`. +3. Hessian in the `petab_problem.compute_hessian` field. It takes two arguments (`H` and `p`) with the updated hessian values being written to `H`. + + +## [Multi-start optimisation](@id petab_multistart_optimisation) +To avoid the optimisation process returning a local minimum, it is often advised to run it multiple times, using different initial guesses. PEtab.jl supports this through the `calibrate_model_multistart` function. This is identical to the `calibrate_model` function, but takes two additional arguments: + +1. `n_multistart`: The number of runs to perform. +2. `dir_save`: A location to which the output is saved. If `dir_save=nothing`, no output is saved. + +And one additional optional argument: + +3. `sampling_method`: Selects the sampling method with which to select the initial guesses (`QuasiMonteCarlo.LatinHypercubeSample()` used by default). + +Because `calibrate_model_multistart` handles initial guess sampling, unlike for `calibrate_model`, no initial guess has to be provided. + +Here, we fit parameters through 10 independent optimisation runs, using QuasiMonteCarlo's `SobolSample` method, and save the result to the OptimizationRuns folder: +``` +using Optim +using QuasiMonteCarlo +res_ms = calibrate_model_multistart(petab_problem, IPNewton(), 10, "OptimizationRuns"; sampling_method=QuasiMonteCarlo.SobolSample()) +``` +The best result across all runs can still be retrieved using `get_ps(res_ms, petab_problem)`, with the results of the individual runs being stored in the `res_ms.runs` field. + +To load the result in a later session, we can call: +``` +res_ms = PEtabMultistartOptimisationResult("OptimizationRuns") +``` +If the OptimizationRuns folder contains the output from several runs, we can designate which to load using the `which_run` argument. Here we loads the second run to be saved in that folder: +``` +res_ms = PEtabMultistartOptimisationResult("OptimizationRuns"; which_run="2") +``` + + +## [Events](@id petab_events) +So far, we have assumed that all experiments, after initiation, run without interference. Experiments where conditions change, or where species are added/removed during the time course, can be represented through events (related to [callbacks](@ref advanced_simulations_callbacks)). In PEtab, an event is represented through the `PEtabEvent` structure. It takes three arguments: +1. The condition for triggering the event. This can either indicate a point in time, or a boolean condition. +2. A rule for updating the event's target +3. The event's target (either a species or parameter). + +Here we create an event which adds `0.5` units of `S` to the system at time `5.0`: +```petab1 +@unpack S = rn +event1 = PEtabEvent(5.0, S + 0.5, S) +``` +Here, the first argument is evaluated to a scalar value, in which case it is interpreted as a time point at which the event happens. If we instead want the event to happen whenever the concentration of `S` falls below `0.5` we set the first argument to a boolean condition indicating this: +```petab1 +event2 = PEtabEvent(S < 0.5, S + 0.5, S) +``` + +Whenever we have several events or not, we bundle them together in a single vector which is later passed to the `PEtabODEProblem` using the `events` argument: +```petab1 +events = [event, event2] +petab_model = PEtabModel(rn, observables, measurements, params; state_map=u0, events=events) +``` + +More details on how to use events, including how to create events with multiple targets, can be found in [PEtab.jl's documentation](https://sebapersson.github.io/PEtab.jl/stable/Julia_event/). + +## [Plot recipes](@id petab_plotting) +There exist various types of graphs that can be used to evaluate the parameter fitting process. These can be plotted using the `plot` command, where the input is either the result of a `calibrate_model` or a `calibrate_model_multistart` run. To be able to use this functionality, you have to ensure that PEtab.jl [records the optimisation process](@id petab_optimisation_path_recording) by providing the `save_trace=true` argument to the calibration functions. + +To, for a single start calibration run, plot, for each iteration of the optimization process, the best objective value achieved so far, run: + +```petab1 +plot(res) +``` +For a multi-start calibration run, the default output is instead a so-called waterfall plot: +```petab1 +plot(res_ms) +``` +![petab waterfall plot](../assets/petab_waterfall.svg) +(for this, and the next plot, we use a multi-start optimisation result from a different model, which yields less trivial optimisation runs than our catalysis one) + +In the waterfall plot, each dot shows the final objective value for a single run in the multi-start process. The runs are ordered by their objective values, and colours designate runs in the same local minimum. + +To instead use the best objective value plot for a multi-start run (with one curve for each run), the `plot_type` argument is used: +```petab1 +plot(res_ms; plot_type = :best_objective) +``` +![petab best objective plot](../assets/petab_best_objective.svg) + +There exist several types of plots for both types of calibration results. More details of the types of available plots, and how to customise them, can be found [here](@ref https://sebapersson.github.io/PEtab.jl/stable/optimisation_output_plotting/). + + +## [Citations](@id petab_citations) +If you use this functionality in your research, [in addition to Catalyst](@ref catalyst_citation), please cite the following papers to support the authors of the PEtab.jl package (currently there is no article associated with this package) and the PEtab standard: +``` +@misc{2023Petabljl, + author = {Ognissanti, Damiano and Arutjunjan, Rafael and Persson, Sebastian and Hasselgren, Viktor}, + title = {{2023Petabljl.jl}}, + howpublished = {\url{https://github.com/sebapersson/PEtab.jl/}}, + year = {2023} +} +``` +``` +@article{SchmiesterSch2021, + author = {Schmiester, Leonard AND Schälte, Yannik AND Bergmann, Frank T. AND Camba, Tacio AND Dudkin, Erika AND Egert, Janine AND Fröhlich, Fabian AND Fuhrmann, Lara AND Hauber, Adrian L. AND Kemmer, Svenja AND Lakrisenko, Polina AND Loos, Carolin AND Merkt, Simon AND Müller, Wolfgang AND Pathirana, Dilan AND Raimúndez, Elba AND Refisch, Lukas AND Rosenblatt, Marcus AND Stapor, Paul L. AND Städter, Philipp AND Wang, Dantong AND Wieland, Franz-Georg AND Banga, Julio R. AND Timmer, Jens AND Villaverde, Alejandro F. AND Sahle, Sven AND Kreutz, Clemens AND Hasenauer, Jan AND Weindl, Daniel}, + journal = {PLOS Computational Biology}, + title = {PEtab—Interoperable specification of parameter estimation problems in systems biology}, + year = {2021}, + month = {01}, + number = {1}, + pages = {1-10}, + volume = {17}, + doi = {10.1371/journal.pcbi.1008646}, + publisher = {Public Library of Science}, + url = {https://doi.org/10.1371/journal.pcbi.1008646}, +} +``` + +## References +[1]: [Schmiester, L et al. *PEtab—Interoperable specification of parameter estimation problems in systems biology*, PLOS Computational Biology (2021).](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1008646) \ No newline at end of file From 5fb6f57a5eccd8ab5446b74a284b2ae35408ec62 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 17 Oct 2023 14:24:25 -0400 Subject: [PATCH 02/32] Update docs/src/inverse_problems/petab_ode_param_fitting.md Co-authored-by: sebapersson <46872750+sebapersson@users.noreply.github.com> --- docs/src/inverse_problems/petab_ode_param_fitting.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index 44d7a7b88c..0215b33722 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -119,7 +119,7 @@ plot!(sol; idxs=4, label="Fitted solution") Here we use the `get_ps` function to retrieve a full parameter set using the optimal parameters (alternatively, the `ODEProblem` or fitted simulation can be retrieved directly using the `get_odeproblem` or `get_odesol` [functions](https://sebapersson.github.io/PEtab.jl/dev/API_choosen/#PEtab.get_odeproblem), respectively). The calibration result can also be found in `res.xmin`, however, note that PEtab automatically ([unless a linear scale is selected](petab_parameters_scales)) converts parameters to logarithmic scale, so typically `10 .^res.xmin` are the values of interest. If you investigate the result from this example you might note that the found parameter set (while it fits the data well) does not correspond to the true parameter set. This phenomenon is related to *identifiability*, and is very important for parameter fitting. ### Final notes -PEtab.jl also supports [multistart optimisation](@ref petab_multistart_optimisation), [automatic equilibration before simulations](https://sebapersson.github.io/PEtab.jl/stable/Brannmark/), and [events](@ref petab_events). Various [plot recipes](@ref petab_plotting) exist for investigating the optimisation process. Please read the [PETab.jl documentation](https://sebapersson.github.io/PEtab.jl/stable/) for a more complete description of the package's features. Below follows additional details of various options and features (generally, PEtab is able to find good default values for most options not specified). +PEtab.jl also supports [multistart optimisation](@ref petab_multistart_optimisation), [automatic pre-equilibration before simulations](https://sebapersson.github.io/PEtab.jl/stable/Brannmark/), and [events](@ref petab_events). Various [plot recipes](@ref petab_plotting) exist for investigating the optimisation process. Please read the [PETab.jl documentation](https://sebapersson.github.io/PEtab.jl/stable/) for a more complete description of the package's features. Below follows additional details of various options and features (generally, PEtab is able to find good default values for most options not specified). ## [Additional features: Observables](@id petab_observables) From aca01250b74e323c5fee51dbde143bd76513f671 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 17 Oct 2023 14:24:35 -0400 Subject: [PATCH 03/32] Update docs/src/inverse_problems/petab_ode_param_fitting.md Co-authored-by: sebapersson <46872750+sebapersson@users.noreply.github.com> --- docs/src/inverse_problems/petab_ode_param_fitting.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index 0215b33722..73de051ec3 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -259,7 +259,7 @@ nothing # hide Note that the `u0` we pass into `PEtabModel` through the `state_map` argument no longer contains the value of *S* (as it is provided by the conditions). ### [Varying parameters between different simulation conditions](@id petab_simulation_conditions) -Sometimes, the parameters that are used vary between the different conditions. Consider our catalysis example, if we had performed the experiment twice, using two different enzymes with different catalytic properties, this could have generated such conditions. The two enzymes could e.g. yield different rates (*kP₁* and *kP₂*) for the `SE --> P + E` reaction, but otherwise be identical. Here, the parameters *kP₁* and *kP₂* are unique to their respective conditions. PEtab.jl provides support for cases such as this, and [its documentation](https://sebapersson.github.io/PEtab.jl/dev/Beer/) provided instructions of how to handle them. +Sometimes, the parameters that are used vary between the different conditions. Consider our catalysis example, if we had performed the experiment twice, using two different enzymes with different catalytic properties, this could have generated such conditions. The two enzymes could e.g. yield different rates (*kP₁* and *kP₂*) for the `SE --> P + E` reaction, but otherwise be identical. Here, the parameters *kP₁* and *kP₂* are unique to their respective conditions. PEtab.jl provides support for cases such as this, and [its documentation](https://sebapersson.github.io/PEtab.jl/stable/Julia_condition_specific/) provided instructions of how to handle them. ## [Additional features: Initial conditions](@id petab_simulation_initial_conditions) From a0d6c27b65ade71730884afdeb2cb89e6700c847 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 17 Oct 2023 21:16:11 +0100 Subject: [PATCH 04/32] update --- .../petab_ode_param_fitting.md | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index 73de051ec3..507d88bb79 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -1,5 +1,5 @@ # [Parameter Fitting for ODEs using PEtab.jl](@id petab_parameter_fitting) -The [PEtab.jl package](https://github.com/sebapersson/PEtab.jl) implements the [PEtab format](https://petab.readthedocs.io/en/latest/) for fitting the parameters of deterministic CRN models to data ([please cite the corresponding papers if you use the PEtab approach in your research](@ref petab_citations)). PEtab.jl both implements methods for creating cost functions (determining how well parameter sets fit to data), and for minimizing these cost functions. The PEtab approach covers most cases of fitting deterministic (ODE) models to data. When possible, it is recommended to use it for parameter fitting (when not, an alternative, more general, approach will have to be used). This page describes how to combine PEtab.jl and Catalyst for parameter fitting, with the PEtab.jl package providing [a more extensive documentation](https://sebapersson.github.io/PEtab.jl/stable/) (this tutorial is partially an adaptation of this documentation). +The [PEtab.jl package](https://github.com/sebapersson/PEtab.jl) implements the [PEtab format](https://petab.readthedocs.io/en/latest/) for fitting the parameters of deterministic CRN models to data ([please cite the corresponding papers if you use the PEtab approach in your research](@ref petab_citations))[^1]. PEtab.jl both implements methods for creating cost functions (determining how well parameter sets fit to data), and for minimizing these cost functions. The PEtab approach covers most cases of fitting deterministic (ODE) models to data. When possible, it is recommended to use it for parameter fitting (when not, an alternative, more general, approach will have to be used). This page describes how to combine PEtab.jl and Catalyst for parameter fitting, with the PEtab.jl package providing [a more extensive documentation](https://sebapersson.github.io/PEtab.jl/stable/) (this tutorial is partially an adaptation of this documentation). ## Introductory example @@ -64,7 +64,7 @@ par_kP = PEtabParameter(:kP) params = [par_kB, par_kD, par_kP] nothing # hide ``` -For each parameter, it is also possible to set [a lower and/or upper bound](@ref petab_parameters_bounds), set whether to use [logarithmic or linear scale](@ref petab_parameters_scales), or add a [prior distribution of its value](@ref petab_parameters_priors). +For each parameter, it is also possible to set [a lower and/or upper bound](@ref petab_parameters_bounds) (by default, *(0.001,1000)* is used), set whether to use [logarithmic or linear scale](@ref petab_parameters_scales), or add a [prior distribution of its value](@ref petab_parameters_priors). ### Simulation conditions Sometimes, several different experiments are performed on a system (each potentially generating several measurements). An experiment could e.g. be the time development of a system from a specific initial condition. Since each experimental condition (during the optimisation procedure, for a guess of the unknown parameters) generates a distinct simulation, these are also called simulation conditions. In our example, all data comes from a single experiment, and the simulation condition input is not required. How to define and use different experimental conditions is described [here](@ref petab_simulation_conditions). @@ -102,6 +102,8 @@ We are now able to fit our model to the data. First, we create a `PEtabODEProble petab_problem = PEtabODEProblem(petab_model) nothing # hide ``` +Since no additional input is given, default optiosn are selected by PEtab.jl (and generally, its choices are good). + To fit a parameter set we use the `calibrate_model` function. In addition to our `PEtabODEProblem`, we must also provide an initial guess (which can be generated with the `generate_startguesses` function) and an optimisation algorithm (which needs to be imported specifically). PEtab.jl supports various [optimisation methods and options](@ref petab_optimisation_optimisers). ```petab1 using Optim @@ -116,7 +118,7 @@ oprob = ODEProblem(rn, u0, (0.0, 10.0), get_ps(res, petab_problem)) sol = solve(oprob, Tsit5(); saveat=0.1) plot!(sol; idxs=4, label="Fitted solution") ``` -Here we use the `get_ps` function to retrieve a full parameter set using the optimal parameters (alternatively, the `ODEProblem` or fitted simulation can be retrieved directly using the `get_odeproblem` or `get_odesol` [functions](https://sebapersson.github.io/PEtab.jl/dev/API_choosen/#PEtab.get_odeproblem), respectively). The calibration result can also be found in `res.xmin`, however, note that PEtab automatically ([unless a linear scale is selected](petab_parameters_scales)) converts parameters to logarithmic scale, so typically `10 .^res.xmin` are the values of interest. If you investigate the result from this example you might note that the found parameter set (while it fits the data well) does not correspond to the true parameter set. This phenomenon is related to *identifiability*, and is very important for parameter fitting. +Here we use the `get_ps` function to retrieve a full parameter set using the optimal parameters. Alternatively, the `ODEProblem` or fitted simulation can be retrieved directly using the `get_odeproblem` or `get_odesol` [functions](https://sebapersson.github.io/PEtab.jl/dev/API_choosen/#PEtab.get_odeproblem), respectively (and the initial condition using the `get_u0` function). The calibration result can also be found in `res.xmin`, however, note that PEtab automatically ([unless a linear scale is selected](petab_parameters_scales)) converts parameters to logarithmic scale, so typically `10 .^res.xmin` are the values of interest. If you investigate the result from this example you might note that the found parameter set (while it fits the data well) does not correspond to the true parameter set. This phenomenon is related to *identifiability*, and is very important for parameter fitting. ### Final notes PEtab.jl also supports [multistart optimisation](@ref petab_multistart_optimisation), [automatic pre-equilibration before simulations](https://sebapersson.github.io/PEtab.jl/stable/Brannmark/), and [events](@ref petab_events). Various [plot recipes](@ref petab_plotting) exist for investigating the optimisation process. Please read the [PETab.jl documentation](https://sebapersson.github.io/PEtab.jl/stable/) for a more complete description of the package's features. Below follows additional details of various options and features (generally, PEtab is able to find good default values for most options not specified). @@ -134,6 +136,8 @@ end obs_X = PEtabObservable(X1 + X2, 0.5) ``` +A common application for this is to define an [*offset* and a *scale* for each observable](https://sebapersson.github.io/PEtab.jl/stable/Define_in_julia/#Defining-the-Observable). + ### [Advanced observables noise formulas](@id petab_observables_noise_formula) In our basic example we assumed that the normally distributed noise had a standard deviation of `0.5`. However, this value may be a parameter (or indeed any algebraic expression). E.g, we could set ```petab1 @@ -144,6 +148,7 @@ and then let the parameter *σ* vary between different [simulation conditions](@ ```petab1 obs_P = PEtabObservable(P, 0.05P) ``` +It would also be possible to make *σ* a *parameter that is fitted as a part of the parameter fitting process*. PEtab.jl assumes that noise is normally distributed (with a standard deviation equal to the second argument provided to `PEtabObservable`). The only other (currently implemented) noise distribution is log-normally distributed noise, which is designated through the `transformation=:log` argument: ```petab1 @@ -175,7 +180,7 @@ to achieve the more conservative bound *[1e-2, 1e2]* for the parameter *kB*. ### [Parameter scales](@id petab_parameters_scales) -Internally, parameters that are fitted are converted to a logarithmic scale (generally, this is considered advantageous). To prevent this conversion, the `scale=:lin` argument can be used. Here we change the scale of *kB* to linear: +Internally, parameters that are fitted are converted to a logarithmic scale (generally, this is considered advantageous[^2]). To prevent this conversion, the `scale=:lin` argument can be used. Here we change the scale of *kB* to linear: ```petab1 par_kB = PEtabParameter(:kB; scale=:lin) @@ -446,7 +451,7 @@ plot(res_ms) ![petab waterfall plot](../assets/petab_waterfall.svg) (for this, and the next plot, we use a multi-start optimisation result from a different model, which yields less trivial optimisation runs than our catalysis one) -In the waterfall plot, each dot shows the final objective value for a single run in the multi-start process. The runs are ordered by their objective values, and colours designate runs in the same local minimum. +In the waterfall plot, each dot shows the final objective value for a single run in the multi-start process. The runs are ordered by their objective values, and colours designate runs in the same local minimum. A common use of waterfall plots is to check whether a sufficient number of runs (typically *>5*) has converged to the same best local minimum (in which case it is assumed to be the *global* minimum). To instead use the best objective value plot for a multi-start run (with one curve for each run), the `plot_type` argument is used: ```petab1 @@ -484,4 +489,5 @@ If you use this functionality in your research, [in addition to Catalyst](@ref c ``` ## References -[1]: [Schmiester, L et al. *PEtab—Interoperable specification of parameter estimation problems in systems biology*, PLOS Computational Biology (2021).](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1008646) \ No newline at end of file +[^1]: [Schmiester, L et al. *PEtab—Interoperable specification of parameter estimation problems in systems biology*, PLOS Computational Biology (2021).](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1008646) +[^2]: [Hass, H et al. *PBenchmark problems for dynamic modeling of intracellular processes*, Bioinformatics (2019).](https://academic.oup.com/bioinformatics/article/35/17/3073/5280731?login=false) \ No newline at end of file From 5d9a1ddb74d6560898fa281b8cd524eb8f91c011 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 24 Oct 2023 08:39:56 -0400 Subject: [PATCH 05/32] updates --- docs/src/inverse_problems/petab_ode_param_fitting.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index 507d88bb79..c4372badcb 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -197,6 +197,7 @@ to set a normally distributed prior (with mean `1.0` and standard deviation `0.2 ```petab1 par_kB = PEtabParameter(:kB; prior=Normal(1.0,0.2), prior_on_linear_scale=false) ``` +In this example, setting `prior_on_linear_scale=false` make sense as a (linear) normal distribution is non-zero for negative values (an alternative is to use a log-normal distribution, e.g. `prior=LogNormal(3.0, 3.0)`). ## [Simulation conditions](@id petab_simulation_conditions) Sometimes, we have data from several different experimental conditions. Here, when a potential parameter set is evaluated during the fitting process, each experimental condition corresponds to one simulation condition (which produces one simulation). To account for this, PEtab permits the user to define several different simulation conditions, with each condition being defined by specific values for some initial conditions and/or parameters. @@ -427,6 +428,7 @@ Here, the first argument is evaluated to a scalar value, in which case it is int ```petab1 event2 = PEtabEvent(S < 0.5, S + 0.5, S) ``` +Here, the event only triggers whenever the condition changes from `false` to `true`, and not while it remains `true` (or when changing from `true` to `false`). E.g. this event only triggers when `S` concentration passes from more than *5.0* to less that *5.0*. Whenever we have several events or not, we bundle them together in a single vector which is later passed to the `PEtabODEProblem` using the `events` argument: ```petab1 @@ -466,7 +468,7 @@ There exist several types of plots for both types of calibration results. More d If you use this functionality in your research, [in addition to Catalyst](@ref catalyst_citation), please cite the following papers to support the authors of the PEtab.jl package (currently there is no article associated with this package) and the PEtab standard: ``` @misc{2023Petabljl, - author = {Ognissanti, Damiano and Arutjunjan, Rafael and Persson, Sebastian and Hasselgren, Viktor}, + author = {Ognissanti, Damiano AND Arutjunjan, Rafael AND Persson, Sebastian AND Hasselgren, Viktor}, title = {{2023Petabljl.jl}}, howpublished = {\url{https://github.com/sebapersson/PEtab.jl/}}, year = {2023} From 861fd8aa0c4214d220032e61a0770f0fa1899584 Mon Sep 17 00:00:00 2001 From: Torkel Date: Sat, 28 Oct 2023 09:47:22 -0400 Subject: [PATCH 06/32] example networks -? Catalyst functionality --- docs/pages.jl | 6 +++--- .../example_networks/basic_CRN_examples.md | 0 .../example_networks/hodgkin_huxley_equation.md | 0 .../example_networks/smoluchowski_coagulation_equation.md | 0 4 files changed, 3 insertions(+), 3 deletions(-) rename docs/src/{ => catalyst_functionality}/example_networks/basic_CRN_examples.md (100%) rename docs/src/{ => catalyst_functionality}/example_networks/hodgkin_huxley_equation.md (100%) rename docs/src/{ => catalyst_functionality}/example_networks/smoluchowski_coagulation_equation.md (100%) diff --git a/docs/pages.jl b/docs/pages.jl index 2b67fd0ba5..46951d1e1b 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -7,9 +7,9 @@ pages = Any["Home" => "index.md", "catalyst_functionality/constraint_equations.md", "catalyst_functionality/parametric_stoichiometry.md", "catalyst_functionality/network_analysis.md", - "Model creation examples" => Any["example_networks/basic_CRN_examples.md", - "example_networks/hodgkin_huxley_equation.md", - "example_networks/smoluchowski_coagulation_equation.md"]], + "Model creation examples" => Any["catalyst_functionality/example_networks/basic_CRN_examples.md", + "catalyst_functionality/example_networks/hodgkin_huxley_equation.md", + "catalyst_functionality/example_networks/smoluchowski_coagulation_equation.md"]], "Catalyst Applications" => Any["catalyst_applications/simulation_structure_interfacing.md", "catalyst_applications/advanced_simulations.md", "catalyst_applications/homotopy_continuation.md", diff --git a/docs/src/example_networks/basic_CRN_examples.md b/docs/src/catalyst_functionality/example_networks/basic_CRN_examples.md similarity index 100% rename from docs/src/example_networks/basic_CRN_examples.md rename to docs/src/catalyst_functionality/example_networks/basic_CRN_examples.md diff --git a/docs/src/example_networks/hodgkin_huxley_equation.md b/docs/src/catalyst_functionality/example_networks/hodgkin_huxley_equation.md similarity index 100% rename from docs/src/example_networks/hodgkin_huxley_equation.md rename to docs/src/catalyst_functionality/example_networks/hodgkin_huxley_equation.md diff --git a/docs/src/example_networks/smoluchowski_coagulation_equation.md b/docs/src/catalyst_functionality/example_networks/smoluchowski_coagulation_equation.md similarity index 100% rename from docs/src/example_networks/smoluchowski_coagulation_equation.md rename to docs/src/catalyst_functionality/example_networks/smoluchowski_coagulation_equation.md From 1f7a1dd920bf529fa75772c71817927d6b2698a8 Mon Sep 17 00:00:00 2001 From: Torkel Date: Sat, 28 Oct 2023 10:30:00 -0400 Subject: [PATCH 07/32] update project file --- docs/Project.toml | 6 ++++++ docs/src/inverse_problems/petab_ode_param_fitting.md | 1 - 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/docs/Project.toml b/docs/Project.toml index 2b8a7745d7..52d47d092e 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,6 +1,7 @@ [deps] BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665" Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" @@ -8,9 +9,11 @@ HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" +Optim = "429524aa-4258-5aef-a3af-852621145aeb" Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +PEtab = "48d54b35-e43e-4a66-a5a1-dde6b987cf69" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" @@ -21,6 +24,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [compat] BifurcationKit = "0.3" Catalyst = "13" +DataFrames = "1" DifferentialEquations = "7.7" Distributions = "0.25" Documenter = "0.27" @@ -28,9 +32,11 @@ HomotopyContinuation = "2.6" Latexify = "0.15, 0.16" ModelingToolkit = "8.47" NonlinearSolve = "1.6.0, 2" +Optim = "1" Optimization = "3.19" OptimizationOptimisers = "0.1.1" OrdinaryDiffEq = "6" +PEtab = "2" Plots = "1.36" SciMLSensitivity = "7.19" Setfield = "1.1" diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index c4372badcb..58108bfd79 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -114,7 +114,6 @@ res = calibrate_model(petab_problem, p0, IPNewton()) We can now simulate our model for the fitted parameter set, and compare the result to the measurements and true solution. ```petab1 oprob = ODEProblem(rn, u0, (0.0, 10.0), get_ps(res, petab_problem)) - sol = solve(oprob, Tsit5(); saveat=0.1) plot!(sol; idxs=4, label="Fitted solution") ``` From 0f5cd0e6653b3ced6783afe0a716fe7b18ada78a Mon Sep 17 00:00:00 2001 From: Torkel Date: Sat, 28 Oct 2023 12:13:34 -0400 Subject: [PATCH 08/32] make gilespie name uniform in example network file plot titles. --- .../example_networks/basic_CRN_examples.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/catalyst_functionality/example_networks/basic_CRN_examples.md b/docs/src/catalyst_functionality/example_networks/basic_CRN_examples.md index 96af01386e..e151a9392d 100644 --- a/docs/src/catalyst_functionality/example_networks/basic_CRN_examples.md +++ b/docs/src/catalyst_functionality/example_networks/basic_CRN_examples.md @@ -34,7 +34,7 @@ jsol = solve(jprob, SSAStepper()) plot(plot(osol; title = "Reaction Rate Equation ODEs"), plot(ssol; title = "Chemical Langevin SDEs"), - plot(jsol; title = "Stochastic Chemical Kinetics Jump Processes"); + plot(jsol; title = "Gillespie Jump Simulation"); layout = (3, 1)) ``` @@ -63,7 +63,7 @@ jprob = JumpProblem(rs, dprob, Direct()) jsol = solve(jprob, SSAStepper()) plot(plot(osol; title = "Reaction Rate Equation ODEs"), - plot(jsol; title = "Stochastic Chemical Kinetics Jump Processes"); + plot(jsol; title = "Gillespie Jump Simulation"); layout = (2, 1)) ``` From d7801cf568226ddb09771b13ecf3f55b86b197d4 Mon Sep 17 00:00:00 2001 From: Torkel Date: Sat, 28 Oct 2023 12:13:42 -0400 Subject: [PATCH 09/32] reference fixes --- docs/src/catalyst_functionality/dsl_description.md | 2 +- .../smoluchowski_coagulation_equation.md | 2 +- .../inverse_problems/petab_ode_param_fitting.md | 14 +++++++------- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/docs/src/catalyst_functionality/dsl_description.md b/docs/src/catalyst_functionality/dsl_description.md index 7a52176f11..76e49076d9 100644 --- a/docs/src/catalyst_functionality/dsl_description.md +++ b/docs/src/catalyst_functionality/dsl_description.md @@ -412,7 +412,7 @@ sol = solve(oprob) plot(sol) ``` -## Setting initial conditions that depend on parameters +## [Setting initial conditions that depend on parameters](@id dsl_description_parametric_initial_conditions) It is possible to set the initial condition of one (or several) species so that they depend on some system parameter. This is done in a similar way as default initial conditions, but giving the parameter instead of a value. When doing this, we also need to ensure that the initial condition parameter is a variable of the system: ```@example tut2 rn = @reaction_network begin diff --git a/docs/src/catalyst_functionality/example_networks/smoluchowski_coagulation_equation.md b/docs/src/catalyst_functionality/example_networks/smoluchowski_coagulation_equation.md index adbacf4ca4..61c0eefd08 100644 --- a/docs/src/catalyst_functionality/example_networks/smoluchowski_coagulation_equation.md +++ b/docs/src/catalyst_functionality/example_networks/smoluchowski_coagulation_equation.md @@ -127,7 +127,7 @@ plot!(ϕ, sol[3,:]/Nₒ, line = (:dot,4,:purple), label="Analytical sol--X3", ``` For the **additive kernel** we find -![additive_kernel](../assets/additive_kernel.svg) +![additive_kernel](../../assets/additive_kernel.svg) --- ## References diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index 58108bfd79..19e7ce4a5a 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -97,7 +97,7 @@ nothing # hide ``` ### Fitting parameters -We are now able to fit our model to the data. First, we create a `PEtabODEProblem`. Here, we use `petab_model` as the only input, but it is also possible to set various [numeric solver and automatic differentiation option](@ref petab_simulation_conditions) (such as method or tolerance). +We are now able to fit our model to the data. First, we create a `PEtabODEProblem`. Here, we use `petab_model` as the only input, but it is also possible to set various [numeric solver and automatic differentiation options](@ref petab_simulation_options) (such as method or tolerance). ```petab1 petab_problem = PEtabODEProblem(petab_model) nothing # hide @@ -117,7 +117,7 @@ oprob = ODEProblem(rn, u0, (0.0, 10.0), get_ps(res, petab_problem)) sol = solve(oprob, Tsit5(); saveat=0.1) plot!(sol; idxs=4, label="Fitted solution") ``` -Here we use the `get_ps` function to retrieve a full parameter set using the optimal parameters. Alternatively, the `ODEProblem` or fitted simulation can be retrieved directly using the `get_odeproblem` or `get_odesol` [functions](https://sebapersson.github.io/PEtab.jl/dev/API_choosen/#PEtab.get_odeproblem), respectively (and the initial condition using the `get_u0` function). The calibration result can also be found in `res.xmin`, however, note that PEtab automatically ([unless a linear scale is selected](petab_parameters_scales)) converts parameters to logarithmic scale, so typically `10 .^res.xmin` are the values of interest. If you investigate the result from this example you might note that the found parameter set (while it fits the data well) does not correspond to the true parameter set. This phenomenon is related to *identifiability*, and is very important for parameter fitting. +Here we use the `get_ps` function to retrieve a full parameter set using the optimal parameters. Alternatively, the `ODEProblem` or fitted simulation can be retrieved directly using the `get_odeproblem` or `get_odesol` [functions](https://sebapersson.github.io/PEtab.jl/dev/API_choosen/#PEtab.get_odeproblem), respectively (and the initial condition using the `get_u0` function). The calibration result can also be found in `res.xmin`, however, note that PEtab automatically ([unless a linear scale is selected](@ref petab_parameters_scales)) converts parameters to logarithmic scale, so typically `10 .^res.xmin` are the values of interest. If you investigate the result from this example you might note that the found parameter set (while it fits the data well) does not correspond to the true parameter set. This phenomenon is related to *identifiability*, and is very important for parameter fitting. ### Final notes PEtab.jl also supports [multistart optimisation](@ref petab_multistart_optimisation), [automatic pre-equilibration before simulations](https://sebapersson.github.io/PEtab.jl/stable/Brannmark/), and [events](@ref petab_events). Various [plot recipes](@ref petab_plotting) exist for investigating the optimisation process. Please read the [PETab.jl documentation](https://sebapersson.github.io/PEtab.jl/stable/) for a more complete description of the package's features. Below follows additional details of various options and features (generally, PEtab is able to find good default values for most options not specified). @@ -263,13 +263,13 @@ nothing # hide ``` Note that the `u0` we pass into `PEtabModel` through the `state_map` argument no longer contains the value of *S* (as it is provided by the conditions). -### [Varying parameters between different simulation conditions](@id petab_simulation_conditions) +### Varying parameters between different simulation conditions Sometimes, the parameters that are used vary between the different conditions. Consider our catalysis example, if we had performed the experiment twice, using two different enzymes with different catalytic properties, this could have generated such conditions. The two enzymes could e.g. yield different rates (*kP₁* and *kP₂*) for the `SE --> P + E` reaction, but otherwise be identical. Here, the parameters *kP₁* and *kP₂* are unique to their respective conditions. PEtab.jl provides support for cases such as this, and [its documentation](https://sebapersson.github.io/PEtab.jl/stable/Julia_condition_specific/) provided instructions of how to handle them. ## [Additional features: Initial conditions](@id petab_simulation_initial_conditions) ### [Fitting initial conditions](@id petab_simulation_initial_conditions_fitted) -Sometimes, initial conditions are uncertain quantities which we wish to fit to the data. This is possible by [defining an initial condition as a parameter](dsl_description_parametric_initial_conditions): +Sometimes, initial conditions are uncertain quantities which we wish to fit to the data. This is possible by [defining an initial condition as a parameter](@ref dsl_description_parametric_initial_conditions): ```petab4 using Catalyst, PEtab # hide rn = @reaction_network begin @@ -298,7 +298,7 @@ nothing # hide and we can use our updated `rn`, `u0`, and `params` as input to our `PEtabModel`. ### [Uncertain initial conditions](@id petab_simulation_initial_conditions_uncertainty) -Often, while an initial condition has been reported for an experiment, its exact value is uncertain. This can be modelled by making the initial condition a [parameter that is fitted to the data](@id petab_simulation_initial_conditions_fitted) and attach a prior to it corresponding to our certainty about its value. +Often, while an initial condition has been reported for an experiment, its exact value is uncertain. This can be modelled by making the initial condition a [parameter that is fitted to the data](@ref petab_simulation_initial_conditions_fitted) and attach a prior to it corresponding to our certainty about its value. Let us consider our initial example, but where we want to add uncertainty to the initial conditions of `S` and `E`. We will add priors on these, assuming normal distributions with mean `1.0` and standard deviation `0.1`. For the synthetic measured data we will use the true values *S(0)=E(0)=1.0*. ```petab5 @@ -438,7 +438,7 @@ petab_model = PEtabModel(rn, observables, measurements, params; state_map=u0, ev More details on how to use events, including how to create events with multiple targets, can be found in [PEtab.jl's documentation](https://sebapersson.github.io/PEtab.jl/stable/Julia_event/). ## [Plot recipes](@id petab_plotting) -There exist various types of graphs that can be used to evaluate the parameter fitting process. These can be plotted using the `plot` command, where the input is either the result of a `calibrate_model` or a `calibrate_model_multistart` run. To be able to use this functionality, you have to ensure that PEtab.jl [records the optimisation process](@id petab_optimisation_path_recording) by providing the `save_trace=true` argument to the calibration functions. +There exist various types of graphs that can be used to evaluate the parameter fitting process. These can be plotted using the `plot` command, where the input is either the result of a `calibrate_model` or a `calibrate_model_multistart` run. To be able to use this functionality, you have to ensure that PEtab.jl [records the optimisation process](@ref petab_optimisation_path_recording) by providing the `save_trace=true` argument to the calibration functions. To, for a single start calibration run, plot, for each iteration of the optimization process, the best objective value achieved so far, run: @@ -460,7 +460,7 @@ plot(res_ms; plot_type = :best_objective) ``` ![petab best objective plot](../assets/petab_best_objective.svg) -There exist several types of plots for both types of calibration results. More details of the types of available plots, and how to customise them, can be found [here](@ref https://sebapersson.github.io/PEtab.jl/stable/optimisation_output_plotting/). +There exist several types of plots for both types of calibration results. More details of the types of available plots, and how to customise them, can be found [here](https://sebapersson.github.io/PEtab.jl/stable/optimisation_output_plotting/). ## [Citations](@id petab_citations) From d55bf904fe1d7a59dd6c863abe9b73ff3258ebd4 Mon Sep 17 00:00:00 2001 From: Torkel Date: Sun, 29 Oct 2023 09:58:05 -0400 Subject: [PATCH 10/32] update --- .../catalyst_applications/simulation_structure_interfacing.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/catalyst_applications/simulation_structure_interfacing.md b/docs/src/catalyst_applications/simulation_structure_interfacing.md index a5ec97571a..d588686466 100644 --- a/docs/src/catalyst_applications/simulation_structure_interfacing.md +++ b/docs/src/catalyst_applications/simulation_structure_interfacing.md @@ -99,7 +99,7 @@ sol[:k1] ``` Finally, we note that we cannot change the values of solution states or parameters (i.e. both `sol[:X1] = 0.0` and `sol[:k1] = 0.0` generate errors). -## Interfacing using symbolic representation +## [Interfacing using symbolic representation](@id simulation_structure_interfacing_symbolic_representation) Catalyst is built on an *intermediary representation* implemented by (ModelingToolkit.jl)[https://github.com/SciML/ModelingToolkit.jl]. ModelingToolkit is a modelling framework where one first declares a set of symbolic variables and parameters using e.g. ```@example ex2 From f04e73e9dd9e649dc27de5394fa3562b1ce8ea86 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Thu, 2 Nov 2023 15:49:58 -0400 Subject: [PATCH 11/32] add error message for parametric stoichiometry example --- docs/src/catalyst_functionality/parametric_stoichiometry.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/catalyst_functionality/parametric_stoichiometry.md b/docs/src/catalyst_functionality/parametric_stoichiometry.md index 87c5696215..e8ffaeff6a 100644 --- a/docs/src/catalyst_functionality/parametric_stoichiometry.md +++ b/docs/src/catalyst_functionality/parametric_stoichiometry.md @@ -77,7 +77,7 @@ sol = solve(oprob, Tsit5()) plot(sol) ``` *If we had used a vector to store parameters, `m` and `n` would be converted to -floating point giving an error when solving the system.* +floating point giving an error when solving the system.* **Note, currently a [bug]() in ModelingToolkit has broken this example by converting to floating point when using tuple parameters, see the alternative approach below for a workaround.** An alternative approach to avoid the issues of using mixed floating point and integer variables is to disable the rescaling of rate laws as described in From a9e0293b836b63fb44ee2bf4a4f8b59369165a67 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Thu, 2 Nov 2023 15:50:44 -0400 Subject: [PATCH 12/32] add link to bug --- docs/src/catalyst_functionality/parametric_stoichiometry.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/catalyst_functionality/parametric_stoichiometry.md b/docs/src/catalyst_functionality/parametric_stoichiometry.md index e8ffaeff6a..a067b9f977 100644 --- a/docs/src/catalyst_functionality/parametric_stoichiometry.md +++ b/docs/src/catalyst_functionality/parametric_stoichiometry.md @@ -77,7 +77,7 @@ sol = solve(oprob, Tsit5()) plot(sol) ``` *If we had used a vector to store parameters, `m` and `n` would be converted to -floating point giving an error when solving the system.* **Note, currently a [bug]() in ModelingToolkit has broken this example by converting to floating point when using tuple parameters, see the alternative approach below for a workaround.** +floating point giving an error when solving the system.* **Note, currently a [bug](https://github.com/SciML/ModelingToolkit.jl/issues/2296) in ModelingToolkit has broken this example by converting to floating point when using tuple parameters, see the alternative approach below for a workaround.** An alternative approach to avoid the issues of using mixed floating point and integer variables is to disable the rescaling of rate laws as described in From 72741c5ee35812046a33abf37836ba1a81f15d58 Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 2 Nov 2023 16:04:03 -0400 Subject: [PATCH 13/32] codebox fix --- .../petab_ode_param_fitting.md | 76 +++++++++---------- 1 file changed, 38 insertions(+), 38 deletions(-) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index 19e7ce4a5a..4a138c1c42 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -4,7 +4,7 @@ The [PEtab.jl package](https://github.com/sebapersson/PEtab.jl) implements the [ ## Introductory example Let us consider a simple catalysis network, where an enzyme (*E*) turns a substrate (*S*) into a product (*P*): -```petab1 +```@example petab1 using Catalyst, PEtab rn = @reaction_network begin @@ -14,7 +14,7 @@ rn = @reaction_network begin end ``` From some known initial condition, and a true parameter set (which we later want to recover from the data) we generate synthetic data (on which we will demonstrate the fitting process). -```petab1 +```@example petab1 # Define initial conditions and parameters. u0 = [:S => 1.0, :E => 1.0, :SE => 0.0, :P => 0.0] p_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5] @@ -43,7 +43,7 @@ The observables define the quantities that we may measure in our experiments. Ty In our example, we only have a single possible observable, the `P` species. We will assume that the noise is normally distributed with a standard deviation `0.5` (in our case this is not true, however, typically the noise distribution is unknown and a guess have to be made). We combine this information in a `PEtabObservable` struct (to access the `P` species we must use [`@unpack`](@ref simulation_structure_interfacing_symbolic_representation)). Finally, we store all our observables in a dictionary, giving each an id tag (which is later used in the measurements input). -```petab1 +```@example petab1 @unpack P = rn obs_P = PEtabObservable(P, 0.5) observables = Dict("obs_P" => obs_P) @@ -57,7 +57,7 @@ Each parameter of the system can either be 3. Be an unknown that we wish to fit to data. In our case, we wish to fit all three system parameters (*kB*, *kD*, and *kP*). For each, we create a single `PEtabParameter`, and then gather these into a single vector. -```petab1 +```@example petab1 par_kB = PEtabParameter(:kB) par_kD = PEtabParameter(:kD) par_kP = PEtabParameter(:kP) @@ -83,7 +83,7 @@ For cases where several simulation conditions are given, we also need to provide ([When [pre-equilibration](https://sebapersson.github.io/PEtab.jl/stable/Brannmark/) is used to initiate the system in a steady state, a fifth field is also required]) Measurements are provided in a `DataFrame` where each row corresponds to a measurement. The values are provided in the `obs_id`, `time`, `measurement`, and `simulation_id` columns. In our case we only need to fill in the first three: -```petab1 +```@example petab1 using DataFrames measurements = DataFrame(obs_id="obs_P", time=sol.t[10:10:end], measurement=data) ``` @@ -91,28 +91,28 @@ Since, in our example, all measurements are of the same observable, we can set ` ### Creating a PEtabModel Finally, we combine all inputs in a single `PEtabModel`. To it, we also pass the initial conditions of our simulations (using the `state_map` argument). It is also possible to have [initial conditions with uncertainty](@ref petab_simulation_initial_conditions_uncertainty), [that vary between different simulations](@ref petab_simulation_conditions), or [that we attempt to fit to the data](@ref petab_simulation_initial_conditions_fitted). -```petab1 +```@example petab1 petab_model = PEtabModel(rn, observables, measurements, params; state_map=u0) nothing # hide ``` ### Fitting parameters We are now able to fit our model to the data. First, we create a `PEtabODEProblem`. Here, we use `petab_model` as the only input, but it is also possible to set various [numeric solver and automatic differentiation options](@ref petab_simulation_options) (such as method or tolerance). -```petab1 +```@example petab1 petab_problem = PEtabODEProblem(petab_model) nothing # hide ``` Since no additional input is given, default optiosn are selected by PEtab.jl (and generally, its choices are good). To fit a parameter set we use the `calibrate_model` function. In addition to our `PEtabODEProblem`, we must also provide an initial guess (which can be generated with the `generate_startguesses` function) and an optimisation algorithm (which needs to be imported specifically). PEtab.jl supports various [optimisation methods and options](@ref petab_optimisation_optimisers). -```petab1 +```@example petab1 using Optim p0 = generate_startguesses(petab_problem, 1) res = calibrate_model(petab_problem, p0, IPNewton()) ``` We can now simulate our model for the fitted parameter set, and compare the result to the measurements and true solution. -```petab1 +```@example petab1 oprob = ODEProblem(rn, u0, (0.0, 10.0), get_ps(res, petab_problem)) sol = solve(oprob, Tsit5(); saveat=0.1) plot!(sol; idxs=4, label="Fitted solution") @@ -126,7 +126,7 @@ PEtab.jl also supports [multistart optimisation](@ref petab_multistart_optimisat ### [Defining non-trivial observables](@id petab_observables_observables) It is possible for observables to be any algebraic expression of species concentrations and parameters. E.g. in this example the total amount of `X` in the system is an observable: -```petab2 +```@example petab2 using Catalyst, PEtab # hide rn = @reaction_network begin (k1,k2), X1 <--> X2 @@ -139,18 +139,18 @@ A common application for this is to define an [*offset* and a *scale* for each o ### [Advanced observables noise formulas](@id petab_observables_noise_formula) In our basic example we assumed that the normally distributed noise had a standard deviation of `0.5`. However, this value may be a parameter (or indeed any algebraic expression). E.g, we could set -```petab1 +```@example petab1 @parameters σ obs_P = PEtabObservable(P, σ) ``` and then let the parameter *σ* vary between different [simulation conditions](@ref petab_simulation_conditions). Alternatively we could let the noise scale linearly with the species' concentration: -```petab1 +```@example petab1 obs_P = PEtabObservable(P, 0.05P) ``` It would also be possible to make *σ* a *parameter that is fitted as a part of the parameter fitting process*. PEtab.jl assumes that noise is normally distributed (with a standard deviation equal to the second argument provided to `PEtabObservable`). The only other (currently implemented) noise distribution is log-normally distributed noise, which is designated through the `transformation=:log` argument: -```petab1 +```@example petab1 obs_P = PEtabObservable(P, σ; transformation=:log) ``` @@ -158,21 +158,21 @@ obs_P = PEtabObservable(P, σ; transformation=:log) ### [Known parameters](@id petab_parameters_known) In our previous example, all parameters were unknowns that we wished to fit to the data. If any parameters have known values, it is possible to provide these to `PEtabModel` through the `parameter_map` argument. E.g if we had known that *kB = 1.0*, then we would only define *kD* and *kP* as parameters we wish to fit: -```petab1 +```@example petab1 par_kD = PEtabParameter(:kD) par_kP = PEtabParameter(:kP) params = [par_kD, par_kP] nothing # hide ``` We then provide `parameter_map=[:kB => 1.0]` when we assembly our model: -```petab1 +```@example petab1 petab_model = PEtabModel(rn, observables, measurements, params; state_map=u0, parameter_map=[:kB => 1.0]) nothing # hide ``` ### [Parameter bounds](@id petab_parameters_bounds) By default, when fitted, potential parameter values are assumed to be in the interval *[1e-3, 1e3]*. When declaring a `PEtabParameter` it is possible to change these values through the `lb` and `ub` arguments. E.g. we could use -```petab1 +```@example petab1 par_kB = PEtabParameter(:kB; lb=1e-2, ub=1e-2) ``` to achieve the more conservative bound *[1e-2, 1e2]* for the parameter *kB*. @@ -181,19 +181,19 @@ to achieve the more conservative bound *[1e-2, 1e2]* for the parameter *kB*. Internally, parameters that are fitted are converted to a logarithmic scale (generally, this is considered advantageous[^2]). To prevent this conversion, the `scale=:lin` argument can be used. Here we change the scale of *kB* to linear: -```petab1 +```@example petab1 par_kB = PEtabParameter(:kB; scale=:lin) ``` ### [Parameter priors](@id petab_parameters_priors) If we has prior knowledge about the distribution of a parameter, it is possible to incorporate this into the model. The prior can be any continuous, univariate, distribution from the [Distributions.jl package](https://github.com/JuliaStats/Distributions.jl). E.g we can use: -```petab1 +```@example petab1 using Distributions par_kB = PEtabParameter(:kB; prior=Normal(1.0,0.2)) ``` to set a normally distributed prior (with mean `1.0` and standard deviation `0.2`) on the value of *kB*. By default, the prior is assumed to be on the linear scale of the parameter (before any potential log transform). To specify that the prior is on the logarithmic scale, the `prior_on_linear_scale=false` argument can be provided: -```petab1 +```@example petab1 par_kB = PEtabParameter(:kB; prior=Normal(1.0,0.2), prior_on_linear_scale=false) ``` In this example, setting `prior_on_linear_scale=false` make sense as a (linear) normal distribution is non-zero for negative values (an alternative is to use a log-normal distribution, e.g. `prior=LogNormal(3.0, 3.0)`). @@ -202,20 +202,20 @@ In this example, setting `prior_on_linear_scale=false` make sense as a (linear) Sometimes, we have data from several different experimental conditions. Here, when a potential parameter set is evaluated during the fitting process, each experimental condition corresponds to one simulation condition (which produces one simulation). To account for this, PEtab permits the user to define several different simulation conditions, with each condition being defined by specific values for some initial conditions and/or parameters. If, for our previous catalysis example, we had measured the system for two different initial values of *S* (*S(0)=1.0* and *S(0)=0.5*), these would correspond to two different simulation conditions. For each condition we define a `Dict` mapping the species to their initial condition (here, *S* is the only species in each `Dict`): -```petab1 +```@example petab1 c1 = Dict(:S => 1.0) c2 = Dict(:S => 0.5) nothing # hide ``` Similarly as for observables, we then gather the conditions in another `Dict`, giving each an id tag: -```petab1 +```@example petab1 simulation_conditions = Dict("c1" => c1, "c2" => c2) nothing # hide ``` Again (like for observables), each measurement in the measurements `DataFrame` needs to be associated with a simulation condition id tag (describing which condition that measurements were taken from). Parameters, just like initial conditions, may vary between different conditions. If an initial condition (or parameter) occurs in one condition, it must occur in all of them. Here follows a complete version of our basic example, but with measurements both for *S(0)=1.0* and *S(0)=0.5*. -```petab3 +```@example petab3 using Catalyst, PEtab rn = @reaction_network begin @@ -270,7 +270,7 @@ Sometimes, the parameters that are used vary between the different conditions. C ### [Fitting initial conditions](@id petab_simulation_initial_conditions_fitted) Sometimes, initial conditions are uncertain quantities which we wish to fit to the data. This is possible by [defining an initial condition as a parameter](@ref dsl_description_parametric_initial_conditions): -```petab4 +```@example petab4 using Catalyst, PEtab # hide rn = @reaction_network begin @parameters E0 @@ -282,12 +282,12 @@ end nothing # hide ``` Here, the initial value of `E` is equal to the parameter `E0`. We modify our `u0` vector by removing `E` (which is no longer known): -```petab4 +```@example petab4 u0 = [:S => 1.0, :SE => 0.0, :P => 0.0] nothing # hide ``` Next, we add `E0` to the parameters we wish to fit: -```petab4 +```@example petab4 par_kB = PEtabParameter(:kB) par_kD = PEtabParameter(:kD) par_kP = PEtabParameter(:kP) @@ -301,7 +301,7 @@ and we can use our updated `rn`, `u0`, and `params` as input to our `PEtabModel` Often, while an initial condition has been reported for an experiment, its exact value is uncertain. This can be modelled by making the initial condition a [parameter that is fitted to the data](@ref petab_simulation_initial_conditions_fitted) and attach a prior to it corresponding to our certainty about its value. Let us consider our initial example, but where we want to add uncertainty to the initial conditions of `S` and `E`. We will add priors on these, assuming normal distributions with mean `1.0` and standard deviation `0.1`. For the synthetic measured data we will use the true values *S(0)=E(0)=1.0*. -```petab5 +```@example petab5 using Catalyst, PEtab @@ -345,7 +345,7 @@ Here, when we fit our data we will also gain values for `S0` and `E0`, however, While in our basic example, we do not provide any additional information to our `PEtabODEProblem`, this is an opportunity to specify how the model should be simulated, and what automatic differentiation techniques to use for the optimisation procedure (if none are provided, appropriate defaults are selected). Here is an example, taken from the [more detailed PEtab.jl documentation](https://sebapersson.github.io/PEtab.jl/dev/Boehm/#Creating-a-PEtabODEProblem) -```petab1 +```@example petab1 petab_problem = PEtabODEProblem(petab_model, ode_solver=ODESolver(Rodas5P(), abstol=1e-8, reltol=1e-8), gradient_method=:ForwardDiff, @@ -358,7 +358,7 @@ where we simulate our ODE model using the `Rodas5p` method (with absolute and re ### [Optimisation methods and options](@id petab_optimisation_optimisers) For our examples, we have used the `Optim.IPNewton` optimisation method. PEtab.jl supports [several additional optimisation methods](https://sebapersson.github.io/PEtab.jl/dev/Avaible_optimisers/#options_optimisers). Furthermore, `calibrate_model`'s `options` argument permits the customisation of the options for any used optimiser. E.g. to designate the maximum number of iterations of the `Optim.IPNewton` method we would use: -```petab1 +```@example petab1 res = calibrate_model(petab_problem, p0, IPNewton(); options=Optim.Options(iterations = 10000)) nothing # hide ``` @@ -372,7 +372,7 @@ To record all the parameter sets evaluated (and their objective values) during t ## Objective function extraction While PEtab.jl provides various tools for analysing the objective function generated by `PEtabODEProblem`, it is also possible to extract this function for customised analysis. Given a `PEtabODEProblem` -```petab1 +```@example petab1 petab_problem = PEtabODEProblem(petab_model) nothing # hide ``` @@ -395,7 +395,7 @@ And one additional optional argument: Because `calibrate_model_multistart` handles initial guess sampling, unlike for `calibrate_model`, no initial guess has to be provided. Here, we fit parameters through 10 independent optimisation runs, using QuasiMonteCarlo's `SobolSample` method, and save the result to the OptimizationRuns folder: -``` +```@example petab1 using Optim using QuasiMonteCarlo res_ms = calibrate_model_multistart(petab_problem, IPNewton(), 10, "OptimizationRuns"; sampling_method=QuasiMonteCarlo.SobolSample()) @@ -403,11 +403,11 @@ res_ms = calibrate_model_multistart(petab_problem, IPNewton(), 10, "Optimization The best result across all runs can still be retrieved using `get_ps(res_ms, petab_problem)`, with the results of the individual runs being stored in the `res_ms.runs` field. To load the result in a later session, we can call: -``` +```@example petab1 res_ms = PEtabMultistartOptimisationResult("OptimizationRuns") ``` If the OptimizationRuns folder contains the output from several runs, we can designate which to load using the `which_run` argument. Here we loads the second run to be saved in that folder: -``` +```@example petab1 res_ms = PEtabMultistartOptimisationResult("OptimizationRuns"; which_run="2") ``` @@ -419,18 +419,18 @@ So far, we have assumed that all experiments, after initiation, run without inte 3. The event's target (either a species or parameter). Here we create an event which adds `0.5` units of `S` to the system at time `5.0`: -```petab1 +```@example petab1 @unpack S = rn event1 = PEtabEvent(5.0, S + 0.5, S) ``` Here, the first argument is evaluated to a scalar value, in which case it is interpreted as a time point at which the event happens. If we instead want the event to happen whenever the concentration of `S` falls below `0.5` we set the first argument to a boolean condition indicating this: -```petab1 +```@example petab1 event2 = PEtabEvent(S < 0.5, S + 0.5, S) ``` Here, the event only triggers whenever the condition changes from `false` to `true`, and not while it remains `true` (or when changing from `true` to `false`). E.g. this event only triggers when `S` concentration passes from more than *5.0* to less that *5.0*. Whenever we have several events or not, we bundle them together in a single vector which is later passed to the `PEtabODEProblem` using the `events` argument: -```petab1 +```@example petab1 events = [event, event2] petab_model = PEtabModel(rn, observables, measurements, params; state_map=u0, events=events) ``` @@ -442,11 +442,11 @@ There exist various types of graphs that can be used to evaluate the parameter f To, for a single start calibration run, plot, for each iteration of the optimization process, the best objective value achieved so far, run: -```petab1 +```@example petab1 plot(res) ``` For a multi-start calibration run, the default output is instead a so-called waterfall plot: -```petab1 +```@example petab1 plot(res_ms) ``` ![petab waterfall plot](../assets/petab_waterfall.svg) @@ -455,7 +455,7 @@ plot(res_ms) In the waterfall plot, each dot shows the final objective value for a single run in the multi-start process. The runs are ordered by their objective values, and colours designate runs in the same local minimum. A common use of waterfall plots is to check whether a sufficient number of runs (typically *>5*) has converged to the same best local minimum (in which case it is assumed to be the *global* minimum). To instead use the best objective value plot for a multi-start run (with one curve for each run), the `plot_type` argument is used: -```petab1 +```@example petab1 plot(res_ms; plot_type = :best_objective) ``` ![petab best objective plot](../assets/petab_best_objective.svg) From f0d2351774c05504474790c6035d9904cdad2acd Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Thu, 2 Nov 2023 16:24:24 -0400 Subject: [PATCH 14/32] cap SciMLBase due to error --- docs/Project.toml | 2 ++ docs/src/catalyst_functionality/parametric_stoichiometry.md | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/Project.toml b/docs/Project.toml index 52d47d092e..500cd929ff 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -15,6 +15,7 @@ OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" PEtab = "48d54b35-e43e-4a66-a5a1-dde6b987cf69" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" @@ -38,6 +39,7 @@ OptimizationOptimisers = "0.1.1" OrdinaryDiffEq = "6" PEtab = "2" Plots = "1.36" +SciMLBase = "~2.5" SciMLSensitivity = "7.19" Setfield = "1.1" SpecialFunctions = "2.1" diff --git a/docs/src/catalyst_functionality/parametric_stoichiometry.md b/docs/src/catalyst_functionality/parametric_stoichiometry.md index a067b9f977..cc1fc66bcd 100644 --- a/docs/src/catalyst_functionality/parametric_stoichiometry.md +++ b/docs/src/catalyst_functionality/parametric_stoichiometry.md @@ -72,7 +72,7 @@ oprob = ODEProblem(osys, u₀, (0.0, 1.0), p) nothing # hide ``` We can now solve and plot the system -```@example s1 +```@julia sol = solve(oprob, Tsit5()) plot(sol) ``` From c548e8eb1e83eea2b83c01abdca5460de3cd5819 Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 2 Nov 2023 17:42:28 -0400 Subject: [PATCH 15/32] Update docs/src/inverse_problems/petab_ode_param_fitting.md Co-authored-by: Sam Isaacson --- docs/src/inverse_problems/petab_ode_param_fitting.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index 4a138c1c42..a0f26f192e 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -1,5 +1,5 @@ # [Parameter Fitting for ODEs using PEtab.jl](@id petab_parameter_fitting) -The [PEtab.jl package](https://github.com/sebapersson/PEtab.jl) implements the [PEtab format](https://petab.readthedocs.io/en/latest/) for fitting the parameters of deterministic CRN models to data ([please cite the corresponding papers if you use the PEtab approach in your research](@ref petab_citations))[^1]. PEtab.jl both implements methods for creating cost functions (determining how well parameter sets fit to data), and for minimizing these cost functions. The PEtab approach covers most cases of fitting deterministic (ODE) models to data. When possible, it is recommended to use it for parameter fitting (when not, an alternative, more general, approach will have to be used). This page describes how to combine PEtab.jl and Catalyst for parameter fitting, with the PEtab.jl package providing [a more extensive documentation](https://sebapersson.github.io/PEtab.jl/stable/) (this tutorial is partially an adaptation of this documentation). +The [PEtab.jl package](https://github.com/sebapersson/PEtab.jl) implements the [PEtab format](https://petab.readthedocs.io/en/latest/) for fitting the parameters of deterministic CRN models to data ([please cite the corresponding papers if you use the PEtab approach in your research](@ref petab_citations))[^1]. PEtab.jl both implements methods for creating cost functions (determining how well parameter sets fit to data), and for minimizing these cost functions. The PEtab approach covers most cases of fitting deterministic (ODE) models to data and is a good default choice when fitting reaction rate equation ODE models. This page describes how to combine PEtab.jl and Catalyst for parameter fitting, with the PEtab.jl package providing [a more extensive documentation](https://sebapersson.github.io/PEtab.jl/stable/) (this tutorial is partially an adaptation of this documentation). ## Introductory example From aea666d7296f05457aea06e5f6cbf92fbb7b0a82 Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 2 Nov 2023 17:42:39 -0400 Subject: [PATCH 16/32] Update docs/src/inverse_problems/petab_ode_param_fitting.md Co-authored-by: Sam Isaacson --- docs/src/inverse_problems/petab_ode_param_fitting.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index a0f26f192e..f0f9fc77a3 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -3,7 +3,7 @@ The [PEtab.jl package](https://github.com/sebapersson/PEtab.jl) implements the [ ## Introductory example -Let us consider a simple catalysis network, where an enzyme (*E*) turns a substrate (*S*) into a product (*P*): +Let us consider a simple catalysis network, where an enzyme ($E$) turns a substrate ($S$) into a product ($P$): ```@example petab1 using Catalyst, PEtab From 6a70665eb34fc4d9ee2cbcad3d5036f5b08857d0 Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 2 Nov 2023 17:43:44 -0400 Subject: [PATCH 17/32] Update docs/src/inverse_problems/petab_ode_param_fitting.md Co-authored-by: Sam Isaacson --- docs/src/inverse_problems/petab_ode_param_fitting.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index f0f9fc77a3..f85bc9b9ba 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -41,7 +41,7 @@ Generally, PEtab takes five different inputs to define an optimisation problem ( ### Observables The observables define the quantities that we may measure in our experiments. Typically, each corresponds to a single species, however, [more complicated observables are possible](@ref petab_observables_observables). For each observable, we also need a noise formula, defining the uncertainty in its measurements. By default, PEtab assumes normally distributed noise, with a mean equal to the true value and a standard deviation which we have to define. It is also possible to use [more advanced noise formulas](@ref petab_observables_noise_formula). -In our example, we only have a single possible observable, the `P` species. We will assume that the noise is normally distributed with a standard deviation `0.5` (in our case this is not true, however, typically the noise distribution is unknown and a guess have to be made). We combine this information in a `PEtabObservable` struct (to access the `P` species we must use [`@unpack`](@ref simulation_structure_interfacing_symbolic_representation)). Finally, we store all our observables in a dictionary, giving each an id tag (which is later used in the measurements input). +In our example, we only have a single possible observable, the `P` species. We will assume that the noise is normally distributed with a standard deviation `0.5` (in our case this is not true, however, as typically the noise distribution is unknown and a guess must be made). We combine this information in a `PEtabObservable` struct (to access the `P` species we must use [`@unpack`](@ref simulation_structure_interfacing_symbolic_representation)). Finally, we store all our observables in a dictionary, giving each an id tag (which is later used in the measurements input). ```@example petab1 @unpack P = rn From 6d642a311ba239fdd5b256a6f3f8aab49b4b6845 Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 2 Nov 2023 17:46:03 -0400 Subject: [PATCH 18/32] Update docs/src/inverse_problems/petab_ode_param_fitting.md Co-authored-by: Sam Isaacson --- docs/src/inverse_problems/petab_ode_param_fitting.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index f85bc9b9ba..3d38ad2720 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -450,6 +450,7 @@ For a multi-start calibration run, the default output is instead a so-called wat plot(res_ms) ``` ![petab waterfall plot](../assets/petab_waterfall.svg) + (for this, and the next plot, we use a multi-start optimisation result from a different model, which yields less trivial optimisation runs than our catalysis one) In the waterfall plot, each dot shows the final objective value for a single run in the multi-start process. The runs are ordered by their objective values, and colours designate runs in the same local minimum. A common use of waterfall plots is to check whether a sufficient number of runs (typically *>5*) has converged to the same best local minimum (in which case it is assumed to be the *global* minimum). From 3bc222835313aedb56e91bed2f9d7bd6dc3e76eb Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 2 Nov 2023 17:46:14 -0400 Subject: [PATCH 19/32] Update docs/src/inverse_problems/petab_ode_param_fitting.md Co-authored-by: Sam Isaacson --- docs/src/inverse_problems/petab_ode_param_fitting.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index 3d38ad2720..d342ab7973 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -406,7 +406,7 @@ To load the result in a later session, we can call: ```@example petab1 res_ms = PEtabMultistartOptimisationResult("OptimizationRuns") ``` -If the OptimizationRuns folder contains the output from several runs, we can designate which to load using the `which_run` argument. Here we loads the second run to be saved in that folder: +If the OptimizationRuns folder contains the output from several runs, we can designate which to load using the `which_run` argument. Here we load the second run to be saved in that folder: ```@example petab1 res_ms = PEtabMultistartOptimisationResult("OptimizationRuns"; which_run="2") ``` From a08503280b7d59ff5ac6690bf9d02a1328f667d3 Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 2 Nov 2023 17:47:00 -0400 Subject: [PATCH 20/32] Update docs/src/inverse_problems/petab_ode_param_fitting.md Co-authored-by: Sam Isaacson --- docs/src/inverse_problems/petab_ode_param_fitting.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index d342ab7973..8241b1f138 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -80,7 +80,7 @@ For cases where several simulation conditions are given, we also need to provide 4. The simulation condition which generates the measurement ([here](@ref petab_simulation_conditions) is an example where this is used). -([When [pre-equilibration](https://sebapersson.github.io/PEtab.jl/stable/Brannmark/) is used to initiate the system in a steady state, a fifth field is also required]) +*Note also, when [pre-equilibration](https://sebapersson.github.io/PEtab.jl/stable/Brannmark/) is used to initiate the system in a steady state, a fifth field is also required.* Measurements are provided in a `DataFrame` where each row corresponds to a measurement. The values are provided in the `obs_id`, `time`, `measurement`, and `simulation_id` columns. In our case we only need to fill in the first three: ```@example petab1 From ada38bb4c2b74e860f010f09ef14ccb882e32790 Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 2 Nov 2023 17:47:50 -0400 Subject: [PATCH 21/32] Update docs/src/inverse_problems/petab_ode_param_fitting.md Co-authored-by: Sam Isaacson --- docs/src/inverse_problems/petab_ode_param_fitting.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index 8241b1f138..f55fe37208 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -102,7 +102,7 @@ We are now able to fit our model to the data. First, we create a `PEtabODEProble petab_problem = PEtabODEProblem(petab_model) nothing # hide ``` -Since no additional input is given, default optiosn are selected by PEtab.jl (and generally, its choices are good). +Since no additional input is given, default options are selected by PEtab.jl (and generally, its choices are good). To fit a parameter set we use the `calibrate_model` function. In addition to our `PEtabODEProblem`, we must also provide an initial guess (which can be generated with the `generate_startguesses` function) and an optimisation algorithm (which needs to be imported specifically). PEtab.jl supports various [optimisation methods and options](@ref petab_optimisation_optimisers). ```@example petab1 From 5514cc9fbac2c47aeb3622d1f2980b8ce4a5d829 Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 2 Nov 2023 17:50:02 -0400 Subject: [PATCH 22/32] Update docs/src/inverse_problems/petab_ode_param_fitting.md Co-authored-by: Sam Isaacson --- docs/src/inverse_problems/petab_ode_param_fitting.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index f55fe37208..689d05ca65 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -82,7 +82,7 @@ For cases where several simulation conditions are given, we also need to provide *Note also, when [pre-equilibration](https://sebapersson.github.io/PEtab.jl/stable/Brannmark/) is used to initiate the system in a steady state, a fifth field is also required.* -Measurements are provided in a `DataFrame` where each row corresponds to a measurement. The values are provided in the `obs_id`, `time`, `measurement`, and `simulation_id` columns. In our case we only need to fill in the first three: +Each individual measurement is provided as a row of a `DataFrame`. The values are provided in the `obs_id`, `time`, `measurement`, and `simulation_id` columns. In our case we only need to fill in the first three: ```@example petab1 using DataFrames measurements = DataFrame(obs_id="obs_P", time=sol.t[10:10:end], measurement=data) From 29b148c5698f6281005e4f2a20835741db6365f2 Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 2 Nov 2023 17:55:41 -0400 Subject: [PATCH 23/32] Update docs/src/inverse_problems/petab_ode_param_fitting.md Co-authored-by: Sam Isaacson --- docs/src/inverse_problems/petab_ode_param_fitting.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index 689d05ca65..9469020771 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -186,7 +186,7 @@ par_kB = PEtabParameter(:kB; scale=:lin) ``` ### [Parameter priors](@id petab_parameters_priors) -If we has prior knowledge about the distribution of a parameter, it is possible to incorporate this into the model. The prior can be any continuous, univariate, distribution from the [Distributions.jl package](https://github.com/JuliaStats/Distributions.jl). E.g we can use: +If we have prior knowledge about the distribution of a parameter, it is possible to incorporate this into the model. The prior can be any continuous, univariate, distribution from the [Distributions.jl package](https://github.com/JuliaStats/Distributions.jl). E.g we can use: ```@example petab1 using Distributions From 4ccb24909fa078b6f5d053a848d635ebe1dd45f6 Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 2 Nov 2023 17:58:41 -0400 Subject: [PATCH 24/32] Update docs/src/inverse_problems/petab_ode_param_fitting.md Co-authored-by: Sam Isaacson --- docs/src/inverse_problems/petab_ode_param_fitting.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index 9469020771..bae68d3143 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -298,7 +298,7 @@ nothing # hide and we can use our updated `rn`, `u0`, and `params` as input to our `PEtabModel`. ### [Uncertain initial conditions](@id petab_simulation_initial_conditions_uncertainty) -Often, while an initial condition has been reported for an experiment, its exact value is uncertain. This can be modelled by making the initial condition a [parameter that is fitted to the data](@ref petab_simulation_initial_conditions_fitted) and attach a prior to it corresponding to our certainty about its value. +Often, while an initial condition has been reported for an experiment, its exact value is uncertain. This can be modelled by making the initial condition a [parameter that is fitted to the data](@ref petab_simulation_initial_conditions_fitted) and attaching a prior to it corresponding to our certainty about its value. Let us consider our initial example, but where we want to add uncertainty to the initial conditions of `S` and `E`. We will add priors on these, assuming normal distributions with mean `1.0` and standard deviation `0.1`. For the synthetic measured data we will use the true values *S(0)=E(0)=1.0*. ```@example petab5 From 947297af2531234aabd8565f74467c4c342ad92d Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 2 Nov 2023 18:03:50 -0400 Subject: [PATCH 25/32] Update docs/src/inverse_problems/petab_ode_param_fitting.md Co-authored-by: Sam Isaacson --- docs/src/inverse_problems/petab_ode_param_fitting.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index bae68d3143..d7a3e9d509 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -196,7 +196,7 @@ to set a normally distributed prior (with mean `1.0` and standard deviation `0.2 ```@example petab1 par_kB = PEtabParameter(:kB; prior=Normal(1.0,0.2), prior_on_linear_scale=false) ``` -In this example, setting `prior_on_linear_scale=false` make sense as a (linear) normal distribution is non-zero for negative values (an alternative is to use a log-normal distribution, e.g. `prior=LogNormal(3.0, 3.0)`). +In this example, setting `prior_on_linear_scale=false` makes sense as a (linear) normal distribution is non-zero for negative values (an alternative is to use a log-normal distribution, e.g. `prior=LogNormal(3.0, 3.0)`). ## [Simulation conditions](@id petab_simulation_conditions) Sometimes, we have data from several different experimental conditions. Here, when a potential parameter set is evaluated during the fitting process, each experimental condition corresponds to one simulation condition (which produces one simulation). To account for this, PEtab permits the user to define several different simulation conditions, with each condition being defined by specific values for some initial conditions and/or parameters. From de35757a9d2ac812fd5dc513632af1d86b49c1cc Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 2 Nov 2023 18:04:31 -0400 Subject: [PATCH 26/32] Update docs/src/inverse_problems/petab_ode_param_fitting.md Co-authored-by: Sam Isaacson --- docs/src/inverse_problems/petab_ode_param_fitting.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index d7a3e9d509..2d5e93de0d 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -212,7 +212,7 @@ Similarly as for observables, we then gather the conditions in another `Dict`, g simulation_conditions = Dict("c1" => c1, "c2" => c2) nothing # hide ``` -Again (like for observables), each measurement in the measurements `DataFrame` needs to be associated with a simulation condition id tag (describing which condition that measurements were taken from). Parameters, just like initial conditions, may vary between different conditions. If an initial condition (or parameter) occurs in one condition, it must occur in all of them. +Again (like for observables), each measurement in the measurements `DataFrame` needs to be associated with a simulation condition id tag (describing which condition those measurements were taken from). Parameters, just like initial conditions, may vary between different conditions. If an initial condition (or parameter) occurs in one condition, it must occur in all of them. Here follows a complete version of our basic example, but with measurements both for *S(0)=1.0* and *S(0)=0.5*. ```@example petab3 From d6ffad0763f82c988d68b1a75e333951cfe7c147 Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 2 Nov 2023 18:39:33 -0400 Subject: [PATCH 27/32] updates --- .../petab_ode_param_fitting.md | 63 +++++++++++-------- 1 file changed, 36 insertions(+), 27 deletions(-) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index 2d5e93de0d..a319067593 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -20,15 +20,17 @@ u0 = [:S => 1.0, :E => 1.0, :SE => 0.0, :P => 0.0] p_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5] # Generate synthetic data. -using OrdinaryDiffEq, Random -oprob = ODEProblem(rn, u0, (0.0, 10.0), p_true) -sol = solve(oprob, Tsit5(); saveat=0.1) -data = (0.8 .+ 0.4*rand(10)) .* sol[:P][10:10:end] +using OrdinaryDiffEq +oprob_true = ODEProblem(rn, u0, (0.0, 10.0), p_true) +true_sol = solve(oprob_true, Tsit5()) +data_sol = solve(oprob_true, Tsit5(); saveat=1.0) +data_ts = data_sol.t[2:end] +data_vals = (0.8 .+ 0.4*rand(10)) .* data_sol[:P][2:end] # Plots the true solutions and the (synthetic data) measurements. using Plots -plot(sol; idxs=:P, label="True solution") -plot!(sol.t[10:10:end], data; label="Measurements", seriestype=:scatter) +plot(true_sol; idxs=:P, label="True solution", lw=8) +plot!(data_ts, data_vals; label="Measurements", seriestype=:scatter, ms=6, color=:blue) ``` Generally, PEtab takes five different inputs to define an optimisation problem (the minimiser of which generates a fitted parameter set): @@ -85,7 +87,7 @@ For cases where several simulation conditions are given, we also need to provide Each individual measurement is provided as a row of a `DataFrame`. The values are provided in the `obs_id`, `time`, `measurement`, and `simulation_id` columns. In our case we only need to fill in the first three: ```@example petab1 using DataFrames -measurements = DataFrame(obs_id="obs_P", time=sol.t[10:10:end], measurement=data) +measurements = DataFrame(obs_id="obs_P", time=data_ts, measurement=data_vals) ``` Since, in our example, all measurements are of the same observable, we can set `obs_id="obs_P"`. If several different observables were measured, `obs_id` would be a vector of the same length as `time` and `measurement`. @@ -113,9 +115,9 @@ res = calibrate_model(petab_problem, p0, IPNewton()) We can now simulate our model for the fitted parameter set, and compare the result to the measurements and true solution. ```@example petab1 -oprob = ODEProblem(rn, u0, (0.0, 10.0), get_ps(res, petab_problem)) -sol = solve(oprob, Tsit5(); saveat=0.1) -plot!(sol; idxs=4, label="Fitted solution") +oprob_fitted = remake(oprob; p=get_ps(res, petab_problem)) +fitted_sol = solve(oprob_fitted, Tsit5()) +plot!(fitted_sol; idxs=:P, label="Fitted solution", linestyle=:dash, lw=6, color=:lightblue) ``` Here we use the `get_ps` function to retrieve a full parameter set using the optimal parameters. Alternatively, the `ODEProblem` or fitted simulation can be retrieved directly using the `get_odeproblem` or `get_odesol` [functions](https://sebapersson.github.io/PEtab.jl/dev/API_choosen/#PEtab.get_odeproblem), respectively (and the initial condition using the `get_u0` function). The calibration result can also be found in `res.xmin`, however, note that PEtab automatically ([unless a linear scale is selected](@ref petab_parameters_scales)) converts parameters to logarithmic scale, so typically `10 .^res.xmin` are the values of interest. If you investigate the result from this example you might note that the found parameter set (while it fits the data well) does not correspond to the true parameter set. This phenomenon is related to *identifiability*, and is very important for parameter fitting. @@ -228,16 +230,16 @@ u0 = [:E => 1.0, :SE => 0.0, :P => 0.0] p_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5] # Simulate data. -using OrdinaryDiffEq, Random -d1, t1 = let - oprob = ODEProblem(rn, [:S=>1.0; u0], (0.0, 10.0), p_true) - sol = solve(oprob, Tsit5(); saveat=0.1) - ((0.8 .+ 0.4*rand(10)) .* sol[:P][10:10:end], sol.t[10:10:end]) +using OrdinaryDiffEq +t1, d1 = let + oprob_true = ODEProblem(rn, [:S=>1.0; u0], (0.0, 10.0), p_true) + data_sol = solve(oprob_true, Tsit5(); saveat=1.0) + data_sol.t[2:end], (0.8 .+ 0.4*rand(10)) .* data_sol[:P][2:end] end -d2, t2 = let - oprob = ODEProblem(rn, [:S=>0.5; u0], (0.0, 10.0), p_true) - sol = solve(oprob, Tsit5(); saveat=0.1) - ((0.8 .+ 0.4*rand(10)) .* sol[:P][10:10:end], sol.t[10:10:end]) +t2, d2 = let + oprob_true = ODEProblem(rn, [:S=>0.5; u0], (0.0, 10.0), p_true) + data_sol = solve(oprob_true, Tsit5(); saveat=1.0) + data_sol.t[2:end], (0.8 .+ 0.4*rand(10)) .* data_sol[:P][2:end] end @unpack P = rn @@ -316,10 +318,12 @@ end u0 = [:SE => 0.0, :P => 0.0] p_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5, :S0=>1.0, :E0 => 1.0] -using OrdinaryDiffEq, Random -oprob = ODEProblem(rn, [:S=>1.0; :E => 1.0; u0], (0.0, 10.0), p_true) -sol = solve(oprob, Tsit5(); saveat=0.1) -data = (0.8 .+ 0.4*rand(10)) .* sol[:P][10:10:end] +using OrdinaryDiffEq +oprob_true = ODEProblem(rn, u0, (0.0, 10.0), p_true) +true_sol = solve(oprob_true, Tsit5()) +data_sol = solve(oprob_true, Tsit5(); saveat=1.0) +data_ts = data_sol.t[2:end] +data_vals = (0.8 .+ 0.4*rand(10)) .* data_sol[:P][2:end] @unpack P = rn obs_P = PEtabObservable(P, 0.5) @@ -333,7 +337,7 @@ par_E0 = PEtabParameter(:E0) params = [par_kB, par_kD, par_kP, par_S0, par_E0] using DataFrames -measurements = DataFrame(obs_id="obs_P", time=sol.t[10:10:end], measurement=data) +measurements = DataFrame(obs_id="obs_P", time=data_ts, measurement=data_vals) petab_model = PEtabModel(rn, observables, measurements, params; state_map=u0) nothing # hide @@ -386,7 +390,7 @@ We can find the: To avoid the optimisation process returning a local minimum, it is often advised to run it multiple times, using different initial guesses. PEtab.jl supports this through the `calibrate_model_multistart` function. This is identical to the `calibrate_model` function, but takes two additional arguments: 1. `n_multistart`: The number of runs to perform. -2. `dir_save`: A location to which the output is saved. If `dir_save=nothing`, no output is saved. +2. `dir_save`: A location to which the output is automatically saved. If `dir_save=nothing`, no output is saved. And one additional optional argument: @@ -406,10 +410,11 @@ To load the result in a later session, we can call: ```@example petab1 res_ms = PEtabMultistartOptimisationResult("OptimizationRuns") ``` -If the OptimizationRuns folder contains the output from several runs, we can designate which to load using the `which_run` argument. Here we load the second run to be saved in that folder: +where `"OptimizationRuns"` is the name of the save directory (specified in `calibrate_model_multistart`). If the OptimizationRuns folder contains the output from several runs, we can designate which to load using the `which_run` argument. Here we load the second run to be saved in that folder: ```@example petab1 res_ms = PEtabMultistartOptimisationResult("OptimizationRuns"; which_run="2") ``` +By default, `which_run` loads the first run saved to that directory. ## [Events](@id petab_events) @@ -441,13 +446,16 @@ More details on how to use events, including how to create events with multiple There exist various types of graphs that can be used to evaluate the parameter fitting process. These can be plotted using the `plot` command, where the input is either the result of a `calibrate_model` or a `calibrate_model_multistart` run. To be able to use this functionality, you have to ensure that PEtab.jl [records the optimisation process](@ref petab_optimisation_path_recording) by providing the `save_trace=true` argument to the calibration functions. To, for a single start calibration run, plot, for each iteration of the optimization process, the best objective value achieved so far, run: - ```@example petab1 plot(res) +nothing # hide ``` +![petab single best objective plot](../assets/petab_best_objective_single_run.svg) + For a multi-start calibration run, the default output is instead a so-called waterfall plot: ```@example petab1 plot(res_ms) +nothing # hide ``` ![petab waterfall plot](../assets/petab_waterfall.svg) @@ -458,6 +466,7 @@ In the waterfall plot, each dot shows the final objective value for a single run To instead use the best objective value plot for a multi-start run (with one curve for each run), the `plot_type` argument is used: ```@example petab1 plot(res_ms; plot_type = :best_objective) +nothing # hide ``` ![petab best objective plot](../assets/petab_best_objective.svg) From 3986a300e4a9bc69f19caaf1f41319ee8ebd13b7 Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 2 Nov 2023 18:54:38 -0400 Subject: [PATCH 28/32] update --- docs/src/inverse_problems/petab_ode_param_fitting.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index a319067593..84f76d2698 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -31,7 +31,9 @@ data_vals = (0.8 .+ 0.4*rand(10)) .* data_sol[:P][2:end] using Plots plot(true_sol; idxs=:P, label="True solution", lw=8) plot!(data_ts, data_vals; label="Measurements", seriestype=:scatter, ms=6, color=:blue) +nothing # hide ``` +![petab data solution](../assets/petab_data.svg) Generally, PEtab takes five different inputs to define an optimisation problem (the minimiser of which generates a fitted parameter set): 1. **Model**: The model which we want to fit to the data (a `ReactionSystem`). @@ -118,7 +120,10 @@ We can now simulate our model for the fitted parameter set, and compare the resu oprob_fitted = remake(oprob; p=get_ps(res, petab_problem)) fitted_sol = solve(oprob_fitted, Tsit5()) plot!(fitted_sol; idxs=:P, label="Fitted solution", linestyle=:dash, lw=6, color=:lightblue) +nothing # hide ``` +![petab fitted solution](../assets/petab_fitted_sol.svg) + Here we use the `get_ps` function to retrieve a full parameter set using the optimal parameters. Alternatively, the `ODEProblem` or fitted simulation can be retrieved directly using the `get_odeproblem` or `get_odesol` [functions](https://sebapersson.github.io/PEtab.jl/dev/API_choosen/#PEtab.get_odeproblem), respectively (and the initial condition using the `get_u0` function). The calibration result can also be found in `res.xmin`, however, note that PEtab automatically ([unless a linear scale is selected](@ref petab_parameters_scales)) converts parameters to logarithmic scale, so typically `10 .^res.xmin` are the values of interest. If you investigate the result from this example you might note that the found parameter set (while it fits the data well) does not correspond to the true parameter set. This phenomenon is related to *identifiability*, and is very important for parameter fitting. ### Final notes From b7b88ece893147533389a770f2d8ee7afd47d056 Mon Sep 17 00:00:00 2001 From: Torkel Date: Fri, 3 Nov 2023 08:54:33 -0400 Subject: [PATCH 29/32] updates --- .../constraint_equations.md | 2 +- .../example_networks/basic_CRN_examples.md | 6 +- .../petab_ode_param_fitting.md | 70 +++++++++++++------ 3 files changed, 54 insertions(+), 24 deletions(-) diff --git a/docs/src/catalyst_functionality/constraint_equations.md b/docs/src/catalyst_functionality/constraint_equations.md index 99ab2ae3a0..29dcbb2fc1 100644 --- a/docs/src/catalyst_functionality/constraint_equations.md +++ b/docs/src/catalyst_functionality/constraint_equations.md @@ -85,7 +85,7 @@ sol = solve(oprob, Tsit5()) plot(sol) ``` -## Adding events +## [Adding events](@id constraint_equations_events) Our current model is unrealistic in assuming the cell will grow exponentially forever. Let's modify it such that the cell divides in half every time its volume reaches a size of `2`. We also assume we lose half of the protein upon diff --git a/docs/src/catalyst_functionality/example_networks/basic_CRN_examples.md b/docs/src/catalyst_functionality/example_networks/basic_CRN_examples.md index e151a9392d..abadae12b6 100644 --- a/docs/src/catalyst_functionality/example_networks/basic_CRN_examples.md +++ b/docs/src/catalyst_functionality/example_networks/basic_CRN_examples.md @@ -34,7 +34,7 @@ jsol = solve(jprob, SSAStepper()) plot(plot(osol; title = "Reaction Rate Equation ODEs"), plot(ssol; title = "Chemical Langevin SDEs"), - plot(jsol; title = "Gillespie Jump Simulation"); + plot(jsol; title = "Stochastic Chemical Kinetics Jump Processes"); layout = (3, 1)) ``` @@ -63,7 +63,7 @@ jprob = JumpProblem(rs, dprob, Direct()) jsol = solve(jprob, SSAStepper()) plot(plot(osol; title = "Reaction Rate Equation ODEs"), - plot(jsol; title = "Gillespie Jump Simulation"); + plot(jsol; title = "Stochastic Chemical Kinetics Jump Processes"); layout = (2, 1)) ``` @@ -89,7 +89,7 @@ jprob = JumpProblem(rs, dprob, Direct()) jsol = solve(jprob, SSAStepper()) plot(plot(osol; title = "Reaction Rate Equation ODEs"), - plot(jsol; title = "Gillespie Jump Simulation"); + plot(jsol; title = "Stochastic Chemical Kinetics Jump Processes"); layout = (2, 1)) ``` diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index 84f76d2698..961fa49dcf 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -60,7 +60,7 @@ Each parameter of the system can either be 2. Depend on experimental/simulation conditions (described [here](@ref petab_simulation_conditions)). 3. Be an unknown that we wish to fit to data. -In our case, we wish to fit all three system parameters (*kB*, *kD*, and *kP*). For each, we create a single `PEtabParameter`, and then gather these into a single vector. +In our case, we wish to fit all three system parameters ($kB$, $kD$, and $kP$). For each, we create a single `PEtabParameter`, and then gather these into a single vector. ```@example petab1 par_kB = PEtabParameter(:kB) par_kD = PEtabParameter(:kD) @@ -68,7 +68,7 @@ par_kP = PEtabParameter(:kP) params = [par_kB, par_kD, par_kP] nothing # hide ``` -For each parameter, it is also possible to set [a lower and/or upper bound](@ref petab_parameters_bounds) (by default, *(0.001,1000)* is used), set whether to use [logarithmic or linear scale](@ref petab_parameters_scales), or add a [prior distribution of its value](@ref petab_parameters_priors). +For each parameter, it is also possible to set [a lower and/or upper bound](@ref petab_parameters_bounds) (by default, $(0.001,1000)$ is used), set whether to use [logarithmic or linear scale](@ref petab_parameters_scales), or add a [prior distribution of its value](@ref petab_parameters_priors). ### Simulation conditions Sometimes, several different experiments are performed on a system (each potentially generating several measurements). An experiment could e.g. be the time development of a system from a specific initial condition. Since each experimental condition (during the optimisation procedure, for a guess of the unknown parameters) generates a distinct simulation, these are also called simulation conditions. In our example, all data comes from a single experiment, and the simulation condition input is not required. How to define and use different experimental conditions is described [here](@ref petab_simulation_conditions). @@ -91,7 +91,8 @@ Each individual measurement is provided as a row of a `DataFrame`. The values ar using DataFrames measurements = DataFrame(obs_id="obs_P", time=data_ts, measurement=data_vals) ``` -Since, in our example, all measurements are of the same observable, we can set `obs_id="obs_P"`. If several different observables were measured, `obs_id` would be a vector of the same length as `time` and `measurement`. +Since, in our example, all measurements are of the same observable, we can set `obs_id="obs_P"`. However, it is also possible [include measurements from several different observables](petab_simulation_measurements_several_observables). + ### Creating a PEtabModel Finally, we combine all inputs in a single `PEtabModel`. To it, we also pass the initial conditions of our simulations (using the `state_map` argument). It is also possible to have [initial conditions with uncertainty](@ref petab_simulation_initial_conditions_uncertainty), [that vary between different simulations](@ref petab_simulation_conditions), or [that we attempt to fit to the data](@ref petab_simulation_initial_conditions_fitted). @@ -150,11 +151,11 @@ In our basic example we assumed that the normally distributed noise had a standa @parameters σ obs_P = PEtabObservable(P, σ) ``` -and then let the parameter *σ* vary between different [simulation conditions](@ref petab_simulation_conditions). Alternatively we could let the noise scale linearly with the species' concentration: +and then let the parameter $σ$ vary between different [simulation conditions](@ref petab_simulation_conditions). Alternatively we could let the noise scale linearly with the species' concentration: ```@example petab1 obs_P = PEtabObservable(P, 0.05P) ``` -It would also be possible to make *σ* a *parameter that is fitted as a part of the parameter fitting process*. +It would also be possible to make $σ$ a *parameter that is fitted as a part of the parameter fitting process*. PEtab.jl assumes that noise is normally distributed (with a standard deviation equal to the second argument provided to `PEtabObservable`). The only other (currently implemented) noise distribution is log-normally distributed noise, which is designated through the `transformation=:log` argument: ```@example petab1 @@ -164,7 +165,7 @@ obs_P = PEtabObservable(P, σ; transformation=:log) ## [Additional features: Parameters](@id petab_parameters) ### [Known parameters](@id petab_parameters_known) -In our previous example, all parameters were unknowns that we wished to fit to the data. If any parameters have known values, it is possible to provide these to `PEtabModel` through the `parameter_map` argument. E.g if we had known that *kB = 1.0*, then we would only define *kD* and *kP* as parameters we wish to fit: +In our previous example, all parameters were unknowns that we wished to fit to the data. If any parameters have known values, it is possible to provide these to `PEtabModel` through the `parameter_map` argument. E.g if we had known that $kB = 1.0$, then we would only define $kD$ and $kP$ as parameters we wish to fit: ```@example petab1 par_kD = PEtabParameter(:kD) par_kP = PEtabParameter(:kP) @@ -178,15 +179,15 @@ nothing # hide ``` ### [Parameter bounds](@id petab_parameters_bounds) -By default, when fitted, potential parameter values are assumed to be in the interval *[1e-3, 1e3]*. When declaring a `PEtabParameter` it is possible to change these values through the `lb` and `ub` arguments. E.g. we could use +By default, when fitted, potential parameter values are assumed to be in the interval $(1e-3, 1e3)$. When declaring a `PEtabParameter` it is possible to change these values through the `lb` and `ub` arguments. E.g. we could use ```@example petab1 par_kB = PEtabParameter(:kB; lb=1e-2, ub=1e-2) ``` -to achieve the more conservative bound *[1e-2, 1e2]* for the parameter *kB*. +to achieve the more conservative bound $(1e-2, 1e2)$ for the parameter $kB$. ### [Parameter scales](@id petab_parameters_scales) -Internally, parameters that are fitted are converted to a logarithmic scale (generally, this is considered advantageous[^2]). To prevent this conversion, the `scale=:lin` argument can be used. Here we change the scale of *kB* to linear: +Internally, parameters that are fitted are converted to a logarithmic scale (generally, this is considered advantageous[^2]). To prevent this conversion, the `scale=:lin` argument can be used. Here we change the scale of $kB$ to linear: ```@example petab1 par_kB = PEtabParameter(:kB; scale=:lin) @@ -199,7 +200,7 @@ If we have prior knowledge about the distribution of a parameter, it is possible using Distributions par_kB = PEtabParameter(:kB; prior=Normal(1.0,0.2)) ``` -to set a normally distributed prior (with mean `1.0` and standard deviation `0.2`) on the value of *kB*. By default, the prior is assumed to be on the linear scale of the parameter (before any potential log transform). To specify that the prior is on the logarithmic scale, the `prior_on_linear_scale=false` argument can be provided: +to set a normally distributed prior (with mean `1.0` and standard deviation `0.2`) on the value of $kB$. By default, the prior is assumed to be on the linear scale of the parameter (before any potential log transform). To specify that the prior is on the logarithmic scale, the `prior_on_linear_scale=false` argument can be provided: ```@example petab1 par_kB = PEtabParameter(:kB; prior=Normal(1.0,0.2), prior_on_linear_scale=false) ``` @@ -208,7 +209,7 @@ In this example, setting `prior_on_linear_scale=false` makes sense as a (linear) ## [Simulation conditions](@id petab_simulation_conditions) Sometimes, we have data from several different experimental conditions. Here, when a potential parameter set is evaluated during the fitting process, each experimental condition corresponds to one simulation condition (which produces one simulation). To account for this, PEtab permits the user to define several different simulation conditions, with each condition being defined by specific values for some initial conditions and/or parameters. -If, for our previous catalysis example, we had measured the system for two different initial values of *S* (*S(0)=1.0* and *S(0)=0.5*), these would correspond to two different simulation conditions. For each condition we define a `Dict` mapping the species to their initial condition (here, *S* is the only species in each `Dict`): +If, for our previous catalysis example, we had measured the system for two different initial values of $S$ ($S(0)=1.0$ and $S(0)=\tfrac{1}{2}$), these would correspond to two different simulation conditions. For each condition we define a `Dict` mapping the species to their initial condition (here, $S$ is the only species in each `Dict`): ```@example petab1 c1 = Dict(:S => 1.0) c2 = Dict(:S => 0.5) @@ -221,7 +222,7 @@ nothing # hide ``` Again (like for observables), each measurement in the measurements `DataFrame` needs to be associated with a simulation condition id tag (describing which condition those measurements were taken from). Parameters, just like initial conditions, may vary between different conditions. If an initial condition (or parameter) occurs in one condition, it must occur in all of them. -Here follows a complete version of our basic example, but with measurements both for *S(0)=1.0* and *S(0)=0.5*. +Here follows a complete version of our basic example, but with measurements both for $S(0)=1.0$ and $S(0)=\tfrac{1}{2}$. ```@example petab3 using Catalyst, PEtab @@ -268,10 +269,36 @@ measurements = vcat(m1,m2) petab_model = PEtabModel(rn, simulation_conditions, observables, measurements, params; state_map=u0) nothing # hide ``` -Note that the `u0` we pass into `PEtabModel` through the `state_map` argument no longer contains the value of *S* (as it is provided by the conditions). +Note that the `u0` we pass into `PEtabModel` through the `state_map` argument no longer contains the value of $S$ (as it is provided by the conditions). + +## [Additional features: Measurements](@id petab_simulation_measurements) + +### [Measurements of several observables](@id petab_simulation_measurements_several_observables) +In our previous example, all our measurements were from a single observable, `obs_P`. If we also had collected measurements of both $S$ and $P$: +```@example petab1 +data_ts = data_sol.t[2:end] +data_vals_S = (0.8 .+ 0.4*rand(10)) .* data_sol[:S][2:end] +data_vals_P = (0.8 .+ 0.4*rand(10)) .* data_sol[:P][2:end] +nothing # hide +``` +and then corresponding observables: +```@example petab1 +@unpack S, P = rn +obs_S = PEtabObservable(S, 0.5) +obs_P = PEtabObservable(P, 0.5) +observables = Dict("obs_S" => obs_P, "obs_P" => obs_P) +nothing # hide +``` +we are able to include all these measurements in the same `measurements` `DataFrame`: +```@example petab1 +m1 = DataFrame(obs_id="obs_P", time=data_ts, measurement=data_vals_S) +m2 = DataFrame(obs_id="obs_S", time=data_ts, measurement=data_vals_P) +measurements = vcat(m1,m2) +``` +which then can be used as input to `PEtabModel`. ### Varying parameters between different simulation conditions -Sometimes, the parameters that are used vary between the different conditions. Consider our catalysis example, if we had performed the experiment twice, using two different enzymes with different catalytic properties, this could have generated such conditions. The two enzymes could e.g. yield different rates (*kP₁* and *kP₂*) for the `SE --> P + E` reaction, but otherwise be identical. Here, the parameters *kP₁* and *kP₂* are unique to their respective conditions. PEtab.jl provides support for cases such as this, and [its documentation](https://sebapersson.github.io/PEtab.jl/stable/Julia_condition_specific/) provided instructions of how to handle them. +Sometimes, the parameters that are used vary between the different conditions. Consider our catalysis example, if we had performed the experiment twice, using two different enzymes with different catalytic properties, this could have generated such conditions. The two enzymes could e.g. yield different rates ($kP_1$ and $kP_2$) for the `SE --> P + E` reaction, but otherwise be identical. Here, the parameters $kP_1$ and $kP_2$ are unique to their respective conditions. PEtab.jl provides support for cases such as this, and [its documentation](https://sebapersson.github.io/PEtab.jl/stable/Julia_condition_specific/) provided instructions of how to handle them. ## [Additional features: Initial conditions](@id petab_simulation_initial_conditions) @@ -307,7 +334,7 @@ and we can use our updated `rn`, `u0`, and `params` as input to our `PEtabModel` ### [Uncertain initial conditions](@id petab_simulation_initial_conditions_uncertainty) Often, while an initial condition has been reported for an experiment, its exact value is uncertain. This can be modelled by making the initial condition a [parameter that is fitted to the data](@ref petab_simulation_initial_conditions_fitted) and attaching a prior to it corresponding to our certainty about its value. -Let us consider our initial example, but where we want to add uncertainty to the initial conditions of `S` and `E`. We will add priors on these, assuming normal distributions with mean `1.0` and standard deviation `0.1`. For the synthetic measured data we will use the true values *S(0)=E(0)=1.0*. +Let us consider our initial example, but where we want to add uncertainty to the initial conditions of `S` and `E`. We will add priors on these, assuming normal distributions with mean `1.0` and standard deviation `0.1`. For the synthetic measured data we will use the true values $S(0) = E(0) = 1.0$. ```@example petab5 using Catalyst, PEtab @@ -386,9 +413,9 @@ petab_problem = PEtabODEProblem(petab_model) nothing # hide ``` We can find the: -1. Objective function in the `petab_problem.compute_cost`. It takes a single argument (`p`) and returns the objective value. -2. Gradient in the `petab_problem.compute_gradient!` field. It takes two arguments (`g` and `p`) with the updated gradient values being written to `g`. -3. Hessian in the `petab_problem.compute_hessian` field. It takes two arguments (`H` and `p`) with the updated hessian values being written to `H`. +1. Objective function as the `petab_problem.compute_cost`. It takes a single argument (`p`) and returns the objective value. +2. Gradient as the `petab_problem.compute_gradient!` field. It takes two arguments (`g` and `p`) with the updated gradient values being written to `g`. +3. Hessian as the `petab_problem.compute_hessian` field. It takes two arguments (`H` and `p`) with the updated hessian values being written to `H`. ## [Multi-start optimisation](@id petab_multistart_optimisation) @@ -437,7 +464,7 @@ Here, the first argument is evaluated to a scalar value, in which case it is int ```@example petab1 event2 = PEtabEvent(S < 0.5, S + 0.5, S) ``` -Here, the event only triggers whenever the condition changes from `false` to `true`, and not while it remains `true` (or when changing from `true` to `false`). E.g. this event only triggers when `S` concentration passes from more than *5.0* to less that *5.0*. +Here, the event only triggers whenever the condition changes from `false` to `true`, and not while it remains `true` (or when changing from `true` to `false`). E.g. this event only triggers when `S` concentration passes from more than $5.0$ to less that $5.0$. Whenever we have several events or not, we bundle them together in a single vector which is later passed to the `PEtabODEProblem` using the `events` argument: ```@example petab1 @@ -447,6 +474,9 @@ petab_model = PEtabModel(rn, observables, measurements, params; state_map=u0, ev More details on how to use events, including how to create events with multiple targets, can be found in [PEtab.jl's documentation](https://sebapersson.github.io/PEtab.jl/stable/Julia_event/). +!!! note + PEtab currently does not support events [created as a part of Catalyst's `ReactionStructure`s](@ref constraint_equations_events) and [implemented through `callbacks`](@ref advanced_simulations_callbacks). Instead, events have to use this interface. + ## [Plot recipes](@id petab_plotting) There exist various types of graphs that can be used to evaluate the parameter fitting process. These can be plotted using the `plot` command, where the input is either the result of a `calibrate_model` or a `calibrate_model_multistart` run. To be able to use this functionality, you have to ensure that PEtab.jl [records the optimisation process](@ref petab_optimisation_path_recording) by providing the `save_trace=true` argument to the calibration functions. @@ -466,7 +496,7 @@ nothing # hide (for this, and the next plot, we use a multi-start optimisation result from a different model, which yields less trivial optimisation runs than our catalysis one) -In the waterfall plot, each dot shows the final objective value for a single run in the multi-start process. The runs are ordered by their objective values, and colours designate runs in the same local minimum. A common use of waterfall plots is to check whether a sufficient number of runs (typically *>5*) has converged to the same best local minimum (in which case it is assumed to be the *global* minimum). +In the waterfall plot, each dot shows the final objective value for a single run in the multi-start process. The runs are ordered by their objective values, and colours designate runs in the same local minimum. A common use of waterfall plots is to check whether a sufficient number of runs (typically $>5$) has converged to the same best local minimum (in which case it is assumed to be the *global* minimum). To instead use the best objective value plot for a multi-start run (with one curve for each run), the `plot_type` argument is used: ```@example petab1 From 576ae139f65c476226b7aafd01aa6ef49c6fdf8f Mon Sep 17 00:00:00 2001 From: Torkel Date: Fri, 3 Nov 2023 09:54:46 -0400 Subject: [PATCH 30/32] Update docs/src/inverse_problems/petab_ode_param_fitting.md Co-authored-by: Sam Isaacson --- docs/src/inverse_problems/petab_ode_param_fitting.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index 961fa49dcf..7cea798faa 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -475,7 +475,7 @@ petab_model = PEtabModel(rn, observables, measurements, params; state_map=u0, ev More details on how to use events, including how to create events with multiple targets, can be found in [PEtab.jl's documentation](https://sebapersson.github.io/PEtab.jl/stable/Julia_event/). !!! note - PEtab currently does not support events [created as a part of Catalyst's `ReactionStructure`s](@ref constraint_equations_events) and [implemented through `callbacks`](@ref advanced_simulations_callbacks). Instead, events have to use this interface. + PEtab currently ignores events [created as a part of a Catalyst `ReactionSystem` model](@ref constraint_equations_events), and does not support SciML-style events [implemented through `callbacks` to `solve`](@ref advanced_simulations_callbacks). Instead, events have to use the preceding interface. ## [Plot recipes](@id petab_plotting) There exist various types of graphs that can be used to evaluate the parameter fitting process. These can be plotted using the `plot` command, where the input is either the result of a `calibrate_model` or a `calibrate_model_multistart` run. To be able to use this functionality, you have to ensure that PEtab.jl [records the optimisation process](@ref petab_optimisation_path_recording) by providing the `save_trace=true` argument to the calibration functions. From 1193c5a575cddc9fc4a9eab78f0a90e9018ebd55 Mon Sep 17 00:00:00 2001 From: Torkel Date: Sat, 4 Nov 2023 17:44:20 -0400 Subject: [PATCH 31/32] rewording --- docs/src/inverse_problems/petab_ode_param_fitting.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index 7cea798faa..e142277e4a 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -125,7 +125,7 @@ nothing # hide ``` ![petab fitted solution](../assets/petab_fitted_sol.svg) -Here we use the `get_ps` function to retrieve a full parameter set using the optimal parameters. Alternatively, the `ODEProblem` or fitted simulation can be retrieved directly using the `get_odeproblem` or `get_odesol` [functions](https://sebapersson.github.io/PEtab.jl/dev/API_choosen/#PEtab.get_odeproblem), respectively (and the initial condition using the `get_u0` function). The calibration result can also be found in `res.xmin`, however, note that PEtab automatically ([unless a linear scale is selected](@ref petab_parameters_scales)) converts parameters to logarithmic scale, so typically `10 .^res.xmin` are the values of interest. If you investigate the result from this example you might note that the found parameter set (while it fits the data well) does not correspond to the true parameter set. This phenomenon is related to *identifiability*, and is very important for parameter fitting. +Here we use the `get_ps` function to retrieve a full parameter set using the optimal parameters. Alternatively, the `ODEProblem` or fitted simulation can be retrieved directly using the `get_odeproblem` or `get_odesol` [functions](https://sebapersson.github.io/PEtab.jl/dev/API_choosen/#PEtab.get_odeproblem), respectively (and the initial condition using the `get_u0` function). The calibration result can also be found in `res.xmin`, however, note that PEtab automatically ([unless a linear scale is selected](@ref petab_parameters_scales)) converts parameters to logarithmic scale, so typically `10 .^res.xmin` are the values of interest. If you investigate the result from this example you might note, even if PEtab.jl have found the global optimum (which fits the data well), this does not actually correspond to the true parameter set. This phenomenon is related to *identifiability*, and is very important for parameter fitting. ### Final notes PEtab.jl also supports [multistart optimisation](@ref petab_multistart_optimisation), [automatic pre-equilibration before simulations](https://sebapersson.github.io/PEtab.jl/stable/Brannmark/), and [events](@ref petab_events). Various [plot recipes](@ref petab_plotting) exist for investigating the optimisation process. Please read the [PETab.jl documentation](https://sebapersson.github.io/PEtab.jl/stable/) for a more complete description of the package's features. Below follows additional details of various options and features (generally, PEtab is able to find good default values for most options not specified). From c3550a966c911a33a1aa43f9d498e263b7b11128 Mon Sep 17 00:00:00 2001 From: Torkel Date: Sat, 4 Nov 2023 17:45:51 -0400 Subject: [PATCH 32/32] re-rewording --- docs/src/inverse_problems/petab_ode_param_fitting.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index e142277e4a..360e8f5025 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -125,7 +125,7 @@ nothing # hide ``` ![petab fitted solution](../assets/petab_fitted_sol.svg) -Here we use the `get_ps` function to retrieve a full parameter set using the optimal parameters. Alternatively, the `ODEProblem` or fitted simulation can be retrieved directly using the `get_odeproblem` or `get_odesol` [functions](https://sebapersson.github.io/PEtab.jl/dev/API_choosen/#PEtab.get_odeproblem), respectively (and the initial condition using the `get_u0` function). The calibration result can also be found in `res.xmin`, however, note that PEtab automatically ([unless a linear scale is selected](@ref petab_parameters_scales)) converts parameters to logarithmic scale, so typically `10 .^res.xmin` are the values of interest. If you investigate the result from this example you might note, even if PEtab.jl have found the global optimum (which fits the data well), this does not actually correspond to the true parameter set. This phenomenon is related to *identifiability*, and is very important for parameter fitting. +Here we use the `get_ps` function to retrieve a full parameter set using the optimal parameters. Alternatively, the `ODEProblem` or fitted simulation can be retrieved directly using the `get_odeproblem` or `get_odesol` [functions](https://sebapersson.github.io/PEtab.jl/dev/API_choosen/#PEtab.get_odeproblem), respectively (and the initial condition using the `get_u0` function). The calibration result can also be found in `res.xmin`, however, note that PEtab automatically ([unless a linear scale is selected](@ref petab_parameters_scales)) converts parameters to logarithmic scale, so typically `10 .^res.xmin` are the values of interest. If you investigate the result from this example you might note, that even if PEtab.jl has found the global optimum (which fits the data well), this does not actually correspond to the true parameter set. This phenomenon is related to the concept of *identifiability*, which is very important for parameter fitting. ### Final notes PEtab.jl also supports [multistart optimisation](@ref petab_multistart_optimisation), [automatic pre-equilibration before simulations](https://sebapersson.github.io/PEtab.jl/stable/Brannmark/), and [events](@ref petab_events). Various [plot recipes](@ref petab_plotting) exist for investigating the optimisation process. Please read the [PETab.jl documentation](https://sebapersson.github.io/PEtab.jl/stable/) for a more complete description of the package's features. Below follows additional details of various options and features (generally, PEtab is able to find good default values for most options not specified).