Catalyst integration cont.

docs/make.jl

@@ -21,14 +21,15 @@ makedocs(;
         "Importing a PEtab problem" => "Boehm.md",
         "Defining a PEtab problem in Julia" => Any["Getting started" => "Define_in_julia.md",
                                                    "Pre-equilibration (steady-state simulations)" => "Julia_steady_state.md",
-                                                   "Noise and observable parameters" => "Julia_obs_noise.md"],
+                                                   "Noise and observable parameters" => "Julia_obs_noise.md", 
+                                                   "Condition specific system/model parameters" => "Julia_condition_specific.md", 
+                                                   "Events (callbacks, dosages etc...)" => "Julia_event.md"],
         "Selecting options for a PEtab-problem" => Any["Models with pre-equilibration (steady-state simulation)" => "Brannmark.md",
                            "Medium sized models and adjoint sensitivity analysis" => "Bachmann.md",
                            "Models with many conditions specific parameters" => "Beer.md"],
         "Parameter estimation" => Any["Parameter estimation" => "Parameter_estimation.md",
                                       "Available optimisers" => "Avaible_optimisers.md",
-                                      "Model selection (PEtab select)" => "Model_selection.md",
-                                      "Condition specific system/model parameters" => "Julia_condition_specific.md"],
+                                      "Model selection (PEtab select)" => "Model_selection.md"],
         "Supported gradient and hessian methods" => "Gradient_hessian_support.md",
         "Choosing the best options for a PEtab problem" => "Best_options.md",
         "API" => "API_choosen.md"

docs/src/API_choosen.md

@@ -5,6 +5,7 @@ PEtabModel
 PEtabODEProblem
 PEtabObservable
 PEtabParameter
+PEtabEvent
 ODESolver
 SteadyStateSolver
 remake_PEtab_problem

docs/src/Define_in_julia.md

@@ -200,6 +200,8 @@ This example has covered the fundamental aspects of setting up a parameter estim
 
 - **Condition Specific System/Model Parameters**: Sometimes a subset of model parameters, like protein synthesis rates, vary between simulation conditions, while other parameters remain constant across all conditions. To handle conditions specific parameters, see [this](@ref define_conditions) tutorial.
 
+- **Events**: Sometimes a model incorporates events like substrate addition at specific time points, and/or parameter changes when a state/species reaches certain values. To manage these events/callbacks, see [this](@ref define_events) tutorial.
+
 For guidance on choosing the best options for your specific PEtab problem, we recommend the [Choosing the Best Options for a PEtab Problem](@ref best_options) section and refer to the [Supported Gradient and Hessian Methods](@ref gradient_support) section for more information on available gradient and hessian methods.
 
 ## Runnable Example

docs/src/Julia_condition_specific.md

@@ -1,6 +1,6 @@
 # [Condition Specific System/Model Parameters](@id define_conditions)
 
-In some cases, certain model such like a substrate enzyme binding rate may vary between experimental conditions while other parameters remain constant. These condition-specific parameters can be defined via the simulation conditions. To demonstrate how, let us consider the same enzyme kinetics model as used in the [Creating a PEtab Parameter Estimation Problem in Julia](@ref define_in_julia) tutorial.
+In some cases, certain model parameters such like a substrate enzyme binding rate may vary between experimental conditions while other parameters remain constant. These condition-specific parameters can be defined via the simulation conditions. To demonstrate how, let us consider the same enzyme kinetics model as used in the [Creating a PEtab Parameter Estimation Problem in Julia](@ref define_in_julia) tutorial.
 
 ```julia
 using Catalyst

docs/src/Julia_event.md

@@ -0,0 +1,200 @@
+# [Events (callbacks, dosages etc...)](@id define_events)
+
+To account for experimental interventions, such as the addition of a substrate, changes in experimental conditions (e.g., temperature), or automatic dosages when a species exceeds a certain concentration, events (often called callbacks, dosages etc...) can be used. In a `PEtabModel`, events can be encoded via `PEtabEvent`:
+
+```julia
+PEtabEvent(condition, affect, target)
+```
+
+where
+
+- `condition` is the event trigger and can be i) a constant value or a model parameter that are assumed to be the time the event is activated, or ii) a Boolean expression triggered when, for example, a state exceeds a certain value.
+
+- `affect` is the effect of the event, which can either be a value or an algebraic expression involving model parameters and/or states.
+
+- `target` specifies the target on which the effect acts. It must be either a model state or parameter.
+
+This section provides examples of how to use `PEtabEvent` to encode different types of events, and uses a modified version of the example in the [Creating a PEtab Parameter Estimation Problem in Julia](@ref define_in_julia) tutorial:
+
+```julia
+using Catalyst
+using DataFrames
+using Distributions
+using PEtab
+
+
+system = @reaction_network begin
+    @parameters se0
+    @species SE(t) = se0  # se0 = initial value for S
+    c1, S + E --> SE
+    c2, SE --> S + E
+    c3, SE --> P + E
+end
+
+# Define state and parameter maps
+state_map =  [:E => 1.0, :P => 0.0]
+parameter_map = [:c1 => 1.0]
+
+# Define observables
+@unpack P, E, SE = system
+obs_P = PEtabObservable(P, 1.0, transformation=:lin)
+obs_Sum = PEtabObservable(E + SE, 3.0, transformation=:log)
+observables = Dict("obs_P" => obs_P,
+                   "obs_Sum" => obs_Sum)
+
+# Define parameters for estimation
+_c3 = PEtabParameter(:c3, scale=:log10)
+_se0 = PEtabParameter(:se0, prior=LogNormal(1.0, 0.5), prior_on_linear_scale=true)
+_c2 = PEtabParameter(:c2)
+petab_parameters = [_c2, _c3, _se0]
+
+# Define simulation conditions
+condition_c0 = Dict(:S => 5.0)
+condition_c1 = Dict(:S => 2.0)
+simulation_conditions = Dict("cond0" => condition_c0,
+                             "cond1" => condition_c1)
+
+# Define measurement data
+measurements = DataFrame(
+    simulation_id=["cond0", "cond0", "cond1", "cond1"],
+    obs_id=["obs_P", "obs_Sum", "obs_P", "obs_Sum"],
+    time=[1.0, 10.0, 1.0, 20.0],
+    measurement=[0.7, 0.1, 1.0, 1.5]
+)
+```
+
+In this section we cover two types of events: those triggered at specific time-points and those triggered by a state (e.g., when a state exceeds a certain concentration).
+
+!!! note
+    Even though events can be directly encoded in a Catalyst or ModellingToolkit model, we recommend `PEtabEvent` for optimal performance (e.g. use `DiscreteCallback` when possible), and for ensuring correct evaluation of the objective function and its derivatives.
+
+## Time-Triggered Events
+
+Time-triggered events are activated at specific time-points. The condition can either a constant value (e.g. 1.0) or a model parameter (e.g. `c2`). For example, to trigger an event at `t = 2` where the value 2 is added to the state `S` write:
+
+```julia
+@unpack S = system
+event = PEtabEvent(2.0, S + 2, S)
+```
+
+When building the `PEtabModel` the event is then provided via the `events` keyword:
+
+```julia
+petab_model = PEtabModel(
+    system, simulation_conditions, observables, measurements,
+    petab_parameters, state_map=state_map, parameter_map=parameter_map,
+    events=event, verbose=true
+)
+petab_problem = PEtabODEProblem(petab_model)
+```
+
+The trigger time can also be a model parameter. For instance, to trigger the event when `t == c2` and to set `c1` to 2.0 write:
+
+```julia
+event = PEtabEvent(:c2, 2.0, :c1)
+```
+
+!!! note
+    If the condition and target are single parameters or states, they can be specied as `Num` (from unpack) or a `Symbol`.If the event involves multiple parameters or states, you must provide them as either a `Num` (as shown below) or a `String`.
+
+### Modifying Event Parameters for Different Simulation Conditions
+
+The trigger trime (`condition`) and/or `affect` can be made specific to different simulation conditions by introducing control parameters (here `c_time` and `c_value`) and setting their values accordingly in the simulation conditions:
+
+```julia
+system = @reaction_network begin
+    @parameters se0 c_time c_value
+    @species SE(t) = se0  # se0 represents the initial value of S
+    c1, S + E --> SE
+    c2, SE --> S + E
+    c3, SE --> P + E
+end
+
+# Define state and parameter maps
+state_map =  [:E => 1.0, :P => 0.0]
+# c_value defaults to 1 prior to the event trigger
+parameter_map = [:c1 => 1.0]
+
+condition_c0 = Dict(:S => 5.0, :c_time => 1.0, :c_value => 2.0)
+condition_c1 = Dict(:S => 2.0, :c_time => 4.0, :c_value => 3.0)
+simulation_conditions = Dict("cond0" => condition_c0,
+                             "cond1" => condition_c1)
+```
+
+In this setup, when the event is defined as:
+
+```julia
+event = PEtabEvent(:c_time, :c_value, :c1)
+```
+
+the `c_time` parameter controls when the event is triggered, so for condition `c0`, the event is triggered at `t=1.0`, while for condition `c1`, it is triggered at `t=4.0`. Additionally, for conditions `cond0` and `cond1`, the parameter `c1` takes on the corresponding `c_value` values, which is `2.0` and `3.0`, respectively.
+
+## State Triggered Events
+
+State-triggered events are activated when a species/state fulfills a certain condition. For example, suppose we have a dosage machine that trigger when the substrate `S` drops below the threshold value of `0.1`, and at this time the machine adds up the substrate so it reaches the value 1.0:
+
+```julia
+@unpack S = system
+event = PEtabEvent(S == 0.2, 1.0, S)
+```
+
+Here, `S == 0.1` means that the event is triggered when `S` reaches the value of 0.1.
+
+With state-triggered events, the direction of the condition can matter. For instance, with `S == 0.2`, the event is triggered whether `S` approaches `0.2` from above or below. If we want it to activate only when `S` enters from above, write:
+
+```julia
+@unpack S = system
+event = PEtabEvent(S < 0.2, 1.0, S)
+```
+
+Here, `S < 0.2` means that the event will trigger only when the expression goes from `false` to `true`, in this case when `S` approaches `0.2` from above. To trigger the event only when `S` approaches `0.2` from below, write `S > 0.2`.
+
+### Modifying Event Parameters for Different Simulation Conditions
+
+For models with multiple simulation conditions, it can be relevant to vary the trigger value and potentially change `affect` value for different simulations. This can be done via by introducing control parameters (here `s_trigger` and `s_vakye`) and setting their values accordingly in the simulation conditions:
+
+```julia
+system = @reaction_network begin
+    @parameters se0 s_trigger s_value
+    @species SE(t) = se0  # se0 = initial value for S
+    c1 * c_controll, S + E --> SE
+    c2, SE --> S + E
+    c3, SE --> P + E
+end
+
+# Define state and parameter maps
+state_map =  [:E => 1.0, :P => 0.0]
+parameter_map = [:c1 => 1.0]
+
+condition_c0 = Dict(:S => 5.0, :s_trigger => 0.2, :s_value => 1.0)
+condition_c1 = Dict(:S => 2.0, :s_trigger => 0.3, :s_value => 2.0)
+simulation_conditions = Dict("cond0" => condition_c0,
+                             "cond1" => condition_c1)
+
+@unpack S, s_trigger, s_value = system
+event = PEtabEvent(S == s_trigger, s_value, S)
+```
+
+Now, for conditions `c0` and `c1`, the event is triggered when `S == 0.2` and `S == 0.3`, respectively. Additionally, the value of `S` changes to 1.0 and 2.0 for conditions `c0` and `c1`, respectively.
+
+## Multiple Events
+
+A model can have multiple events, which can be easily defined as a `PEtabModel` accepts a `Vector` of `PEtabEvent` as input. For example, suppose we have an event triggered when the substrate `S` satisfies `S < 0.2`, where `S` changes its value to `1.0`. Additionally, we have another event triggered when `t == 1.0`, where the parameter `c1` changes its value to `2.0`. This can be encoded as follows:
+
+```julia
+@unpack S, c1 = system
+event1 = PEtabEvent(S < 0.2, 1.0, S)
+event2 = PEtabEvent(1.0, 2.0, :c1)
+events = [event1, event2]
+```
+
+These events can then be provided when building the `PEtabModel` with the `events` keyword:
+
+```julia
+petab_model = PEtabModel(
+    system, simulation_conditions, observables, measurements,
+    petab_parameters, state_map=state_map, parameter_map=parameter_map,
+    events=events, verbose=true
+)
+petab_problem = PEtabODEProblem(petab_model)
+```

ext/CatalystExtension/Create_PEtab_model.jl

@@ -5,24 +5,26 @@ function PEtab.PEtabModel(system::ReactionSystem,
                           petab_parameters::Vector{PEtab.PEtabParameter};
                           state_map::Union{Nothing, Vector{Pair{T1, Float64}}}=nothing,
                           parameter_map::Union{Nothing, Vector{Pair{T2, Float64}}}=nothing,
-                          verbose::Bool=false)::PEtab.PEtabModel where {T1<:Union{Symbol, Any}, T2<:Union{Symbol, Any}, T<:Dict}
+                          events::Union{T3, Vector{T3}, Nothing}=nothing,
+                          verbose::Bool=false)::PEtab.PEtabModel where {T1<:Union{Symbol, Any}, T2<:Union{Symbol, Any}, T<:Dict, T3<:PEtabEvent}
 
     model_name = "ReactionSystemModel"                          
     return PEtab._PEtabModel(system, model_name, simulation_conditions, observables, measurements, 
-                             petab_parameters, state_map, parameter_map, verbose)
+                             petab_parameters, state_map, parameter_map, events, verbose)
 end
 function PEtab.PEtabModel(system::ReactionSystem,
                           observables::Dict{String, PEtab.PEtabObservable},
                           measurements::DataFrame,
                           petab_parameters::Vector{PEtab.PEtabParameter};
                           state_map::Union{Nothing, Vector{Pair{T1, Float64}}}=nothing,
                           parameter_map::Union{Nothing, Vector{Pair{T2, Float64}}}=nothing,
-                          verbose::Bool=false)::PEtab.PEtabModel where {T1<:Union{Symbol, Any}, T2<:Union{Symbol, Any}}
+                          events::Union{T3, Vector{T3}, Nothing}=nothing,
+                          verbose::Bool=false)::PEtab.PEtabModel where {T1<:Union{Symbol, Any}, T2<:Union{Symbol, Any}, T3<:PEtabEvent}
 
     simulation_conditions = Dict("__c0__" => Dict())                        
     model_name = "ReactionSystemModel"                          
     return PEtab._PEtabModel(system, model_name, simulation_conditions, observables, measurements, 
-                             petab_parameters, state_map, parameter_map, verbose)
+                             petab_parameters, state_map, parameter_map, events, verbose)
 end
 
 

ext/FidesExtension/Setup_fides.jl

@@ -59,6 +59,7 @@ function PEtab.calibrate_model(petab_problem::PEtabODEProblem,
                                    fmin, 
                                    _p0,
                                    xmin, 
+                                   petab_problem.θ_names,
                                    converged, 
                                    runtime)
 end

ext/IpoptExtension/Setup_ipopt.jl

@@ -60,6 +60,7 @@ function PEtab.calibrate_model(petab_problem::PEtabODEProblem,
                                    fmin, 
                                    _p0,
                                    xmin, 
+                                   petab_problem.θ_names,
                                    converged, 
                                    runtime)
 end

ext/OptimExtension/Setup_optim.jl

@@ -110,6 +110,7 @@ function PEtab.calibrate_model(petab_problem::PEtabODEProblem,
                                    fmin, 
                                    _p0,
                                    xmin, 
+                                   petab_problem.θ_names,
                                    converged, 
                                    runtime)
 end

src/Model_callibration/Common.jl

@@ -306,6 +306,7 @@ function _calibrate_model_multistart(petab_problem::PEtabODEProblem,
     xmin = res_best.xmin
     sampling_method_str = string(sampling_method)[1:findfirst(x -> x == '(', string(sampling_method))][1:end-1]
     results = PEtabMultistartOptimisationResult(xmin,
+                                                petab_problem.θ_names,
                                                 fmin,
                                                 n_multistarts,
                                                 res_best.alg,

src/PEtab.jl

@@ -98,7 +98,7 @@ include(joinpath("Show.jl"))
     end
 end
 
-export PEtabModel, PEtabODEProblem, ODESolver, SteadyStateSolver, PEtabModel, PEtabODEProblem, remake_PEtab_problem, Fides, solve_SBML, PEtabOptimisationResult, IpoptOptions, IpoptOptimiser, PEtabParameter, PEtabObservable, PEtabMultistartOptimisationResult, generate_startguesses, get_fitted_ps, get_fitted_u0
+export PEtabModel, PEtabODEProblem, ODESolver, SteadyStateSolver, PEtabModel, PEtabODEProblem, remake_PEtab_problem, Fides, solve_SBML, PEtabOptimisationResult, IpoptOptions, IpoptOptimiser, PEtabParameter, PEtabObservable, PEtabMultistartOptimisationResult, generate_startguesses, get_fitted_ps, get_fitted_u0, PEtabEvent
 
 # These are given as extensions, but their docstrings are availble in the 
 # general documentation