From ef5b9fa74d913d45d27b73dadd7ca9439fba84c3 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 30 Apr 2024 11:31:21 -0400 Subject: [PATCH 1/3] init --- .../events_and_callbacks.md | 405 ++++++++++++++++++ 1 file changed, 405 insertions(+) create mode 100644 docs/src/catalyst_functionality/events_and_callbacks.md diff --git a/docs/src/catalyst_functionality/events_and_callbacks.md b/docs/src/catalyst_functionality/events_and_callbacks.md new file mode 100644 index 0000000000..09fd8b651d --- /dev/null +++ b/docs/src/catalyst_functionality/events_and_callbacks.md @@ -0,0 +1,405 @@ +# [Simulation events](@id events) + +Sometimes one wishes to model events, describing things that can happen to a system during its simulation. Examples can include: + - A chemical system where an amount of some species is added at a time point after the simulation's initiation. + - A simulation of a circadian rhythm, where light is turned on/off every 12 hours. + - A cell divides when some size variable reaches a certain threshold, halving the amount of each species in the system. + +Generally, there are two approaches for creating events: +1. Including them as a part of the model when its is created. +2. Creating a separate *callback*, which is then supplied to the `solve` function. + +Generally, the first approach is more convenient, and should be used if possible. However, the second approach is more potent, and can be used to create more general types of events. Below, we will describe both approaches. + +## [Creating models with incorporated events](@id events_model) +Sometimes one wishes to model events, describing things that can happen to it during a simulation. + - A chemical system where an amount of some species is added at a time point after the simulation's initiation. + - A simulation of a circadian rhythm, where light is turned on/off every 12 hours. + - A cell divides when some size variable reaches a certain threshold, halving the amount of each species in the system. + +Events are divided into *continuous* and *discrete* events, and these can be added directly to a system using the `continuous_events` and `discrete_events` options. Events can also be modelled through *callbacks*. These are different in that they are supplied in the simulation step (rather than on system creation), and generally provide more flexibility in how they may affect the system. Callbacks are described on a separate [page](@ref advanced_simulations_callbacks). + +The notation described below for creating continuous and discrete events is the same which is used in [ModelingToolkit to create events](https://docs.sciml.ai/ModelingToolkit/stable/basics/Events/), and which is used for [events for programmatic model creation](@ref ref). + +### [Continuous vs discrete events](@id events_model_continuous_vs_discrete) +Both continuous and discrete events combine some condition (for triggering the event) with some affect (describing their effects on the system). They differ in the following ways: +- They use slightly different notation. +- Discrete events' conditions are checked at *the end of* each simulation time step. For continuous events, the simulation instead finds the *exact time point* when the event is triggered at. +- Continuous events cannot be supplied to jump simulations. + +### [Continuous events](@id events_model_continuous) +Let us consider a simple system where species `X` degraded at a constant rate `d`. Next, we wish to add an event which adds `2.0` units of `X` whenever `X` reaches a critical threshold `1.0`. This can be done in the following manner: +```@example dsl_advanced_events +using Catalyst # hide +rn = @reaction_network begin + @continuous_events begin + X ~ 1.0 => [X ~ X + 2.0] + end + d, X --> 0 +end +nothing # hide +``` +Here, the `@continuous_events` option is followed by a `begin ... end` block. Next, each line corresponds to a separate event. Each event is created in the following manner: +- It combines a *condition* (denoting when the event will happen) with one (or several) *affects* (denoting what the event does to the system when it happens). +- The condition (left) and affect (right) are separated by a `=>`. +- The condition takes the form of an [equation](). Here, the event is triggered whenever the equation's two sides (separated by a `~`) are equal. +- The affect(s) are enclosed within `[]`. If there are multiple affects, these are separated by `,` (the example above contains a single affect). +- Each affect is a single equation that describes how a parameter, species, or [variable](@ref dsl_advanced_options_variables) is updated when the event is triggered. +- Each affect's equation's left-hand side must contain only the parameter/species/variable whose value should be updated. +- Each affect's equation's right-hand side is an expression describing its updated value. + +We can simulate the model we declared, just like any other model: +```@example dsl_advanced_events +using OrdinaryDiffEq, Plots +u0 = [:X => 2.0] +tspan = (0.0, 10.0) +ps = [:d => 1.0] + +oprob = ODEProblem(rn, u0, tspan, ps) +sol = solve(ODEProblem, Tsit5()) +plot(sol) +``` +Inspecting the solution, we can confirm that whenever `X` reaches a value of `1.0`, `2.0` units of `X` is added to the system. + +In our example, we can also denote the critical quantities using parameters: +```@example dsl_advanced_events +rn = @reaction_network begin + @parameters X_thres X_add + @continuous_events begin + X ~ X_thres => [X ~ X + X_add] + end + d, X --> 0 +end +nothing # hide +``` +Here, since `X_thres` and `X_add` do not appear in any reactions, Catalyst cannot determine whether they are parameters or species. hence, they must be [explicitly designated as parameters by using the `@parameters` option](@ref dsl_advanced_options_declaring_species_and_parameters). Next, these can be designated as any value, and supplied to the `ODEProblem`: +```@example dsl_advanced_events +u0 = [:X => 2.0] +tspan = (0.0, 10.0) +ps = [:d => 1.0, :X_thres => 0.5, :X_add => 3.0] + +oprob = ODEProblem(rn, u0, tspan, ps) +sol = solve(ODEProblem, Tsit5()) +plot(sol) +``` + +As previously noted, each continuous event can have multiple affects. The following system has two components (`X` and `Y`, one being produced and one being degraded). When their concentrations are equal, a continuous events reduce the concentration of `X` while increasing the concentration of `Y`: +```@example dsl_advanced_events +rn = @reaction_network begin + @continuous_events begin + X ~ Y => [X ~ X - 1.0, Y ~ Y + 1.0] + end + p, 0 --> X + d, Y --> 0 +end + +u0 = [:X => 1.0, :Y => 3.0] +tspan = (0.0, 10.0) +ps = [:p => 1.0, :d => 1.0] + +oprob = ODEProblem(rn, u0, tspan, ps) +sol = solve(ODEProblem, Tsit5()) +plot(sol) +``` + +!!!warn + A single event (continuous or discrete) can update the value of either (one or several) species (and variables), or of (one or several) parameters. It is not possible for an event to update the values of both species/variables and parameters. + +In the above examples we have modelled a system with a single event. In these cases, the `begin end` block is not required, and the event can be provided on the same line as the `@continuous_events` option: +```@example dsl_advanced_events +rn = @reaction_network begin + @continuous_events X ~ Y => [X ~ X - 1.0, Y ~ Y + 1.0] + p, 0 --> X + d, Y --> 0 +end +nothing # hide +``` + +### [Discrete events](@id events_model_discrete) +Just like [continuous events](dsl_advanced_options_events_continuous), discrete events combine a condition with one or more affect statements. Here, discrete events' affects are created identically to those for continuous events. Discrete events' conditions are different. There exist 3 different types of discrete events, each with a different type of condition. All three types are created using the `@discrete_events` option, and a single system can contain a mix of all types. The three types are: +- Preset-time discrete events. +- Periodic discrete events. +- Conditional discrete events. + +#### [Preset-time discrete events](@id events_model_discrete_presettime) +*Present-time events* are events that happen at specific time points. Here, the condition is a vector with all the time points at which an event is triggered. E.g. here we create a production/degradation loop, where `2.0` units of `X` is added at time points `3.0` and `7.0` +```@example dsl_advanced_events +rn = @reaction_network begin + @discrete_events begin + [3.0, 7.0] => [X ~ X + 2.0] + end + (p,d), 0 <--> X +end + +u0 = [:X => 0.1, :Y => 3.0] +tspan = (0.0, 10.0) +ps = [:p => 1.0, :d => 0.5] + +oprob = ODEProblem(rn, u0, tspan, ps) +sol = solve(ODEProblem, Tsit5()) +plot(sol) +``` + +The preset time points can also be parameters (in which case, they have to be [designated as such using the `@parameters` option](@ref dsl_advanced_options_declaring_species_and_parameters)): +```@example dsl_advanced_events +rn = @reaction_network begin + @parameters t1 t2 + @discrete_events begin + [t1, t2] => [X ~ X + 2.0] + end + (p,d), 0 <--> X +end +nothing +``` + +#### [Periodic discrete events](@id events_model_discrete_periodic) +When a discrete event's condition is a vector, a preset-time event is created. If it instead is a single value, a *periodic event* is created. These occur repeatedly throughout a simulation, with its period set by the affect term. E.g. here we create a system where `0.5` units of `X` is added every `1.0` time units. +```@example dsl_advanced_events +rn = @reaction_network begin + @discrete_events begin + 1.0 => [X ~ X + 0.5] + end + d, X --> 0 +end + +u0 = [:X => 1.0] +tspan = (0.0, 10.0) +ps = [:d => 1.0] + +oprob = ODEProblem(rn, u0, tspan, ps) +sol = solve(ODEProblem, Tsit5()) +plot(sol) +``` +Like for preset-time events, periodic events' affects may contain parameters. + +#### [Conditional discrete events](@id events_model_discrete_conditional) +Finally, discrete events' condition may be a boolean expression (consisting of parameters, species, variables, and the time variable). Let's say that we want to create an event which, if the concentration of `X` is below a threshold `1.0`, adds `1.0` units of `X` to the system, then we can use the following discrete event: +```@example dsl_advanced_events +rn = @reaction_network begin + @discrete_events begin + X < 1.0 => [X ~ X + 2.0] + end + d, X --> 0 +end +``` +If we simulate the system using the same conditions as for our [similar, continuous, example](@ref dsl_advanced_options_events_continuous) we get the same result: +```@example dsl_advanced_events +u0 = [:X => 2.0] +tspan = (0.0, 10.0) +ps = [:d => 1.0] + +oprob = ODEProblem(rn, u0, tspan, ps) +sol = solve(ODEProblem, Tsit5()) +plot(sol) +``` +So, how is modelling this event as a discrete or continuous event different? There are four differences: +1) For continuous events, the simulation method finds the exact time point when the condition triggers. Discrete events are triggered at the first time step when the condition holds. +2) This discrete event will be triggered whenever `X < 1.0` holds, not just when the concentration of `X` passes the threshold. E.g. it will be triggered if the initial concentration of `X` is less than `1.0`. +3) Only the discrete event event can be used with jump simulations. +4) The discrete event can be used to create more advanced conditions. + +E.g. using (4), we can modify our system so that the event is only triggered when time is less than `5.0` (after which `X` decays towards `0`): +```@example dsl_advanced_events +rn = @reaction_network begin + @discrete_events begin + (X < 1.0) & (t < 5.0) => [X ~ X + 2.0] + end + d, X --> 0 +end + +u0 = [:X => 2.0] +tspan = (0.0, 10.0) +ps = [:d => 1.0] + +oprob = ODEProblem(rn, u0, tspan, ps) +sol = solve(ODEProblem, Tsit5()) +plot(sol) +``` +!!!note + When we form composite boolean conditions for conditional discrete events, we use `&` to denote the AND operator (not `&&`, as this is currently not supported). + +!!!warn + Generally, discrete events including equality (`==`) will not be triggered. The reason is that the condition is only checked at the end of every time step. Hence, unless special precautions are taken to ensure that the simulator stops when the condition holds, it will unlikely be triggered. + + +## [Creating event through callbacks](@id events_model_callbacks) + +### [Discrete callbacks](@id events_model_callbacks_discrete) + +### [Present time callbacks](@id events_model_callbacks_presenttime) + +### [Continuous callbacks](@id events_model_callbacks_continuous) + +### [Callbacks during jump simulations](@id events_model_callbacks_jump_callbacks) + + +Sometimes one wishes to add discrete events during simulations. Examples could include: + - A chemical system where an amount of some species is added at a time point + after the simulation's initiation. + - A simulation of a circadian rhythm, where light is turned on/off every 12 hours. + - A cell divides when some size variable reaches a certain threshold, randomly + allocating all species to two daughter cells. + +In simple cases events such as these can be modelled symbolically, as described +in the [Constraint Equations and Events](@ref constraint_equations) tutorial. A +more flexible, but low-level, interface is also available via the callback +functionality of DifferentialEquations.jl. A callback is a function that is +passed to the `solve()` command, combing an `affect!` function (defining how the +callback changes the system) with a `condition` function (a condition for +triggering a callback). For a thorough introduction, please read [the section +about callbacks in the DifferentialEquations.jl +documentation](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/). + +There exist three types of callbacks, `PresetTimeCallback`s `DiscreteCallback`s, +and `ContinuousCallback`s. Here, we will limit ourselves to introducing the +`PresetTimeCallback`. For our example, we are going to use a simple network +where a single component, `X`, degrades linearly. +```@example ex2 +using Catalyst +degradation_model = @reaction_network begin + d, X --> 0 +end +nothing # hide +``` +we can simulate the model without using a callback: +```@example ex2 +using DifferentialEquations, Plots +u0 = [:X => 10.0] +tspan = (0.0, 10.0) +p = [:d => 1.0] + +oprob = ODEProblem(degradation_model, u0, tspan, p) +sol = solve(oprob) +plot(sol) +``` +We now wish to modify our simulation so that at the times `t = 3.0` and `t = +7.0` we add `5` units of `X` to the system. For this we create a +`PresetTimeCallback`: +```@example ex2 +condition = [3.0, 7.0] +function affect!(integrator) + integrator[:X] += 5.0 +end +ps_cb = PresetTimeCallback(condition, affect!) +nothing # hide +``` +Here, `condition` is simply a vector with all the time points during which we +want the callback to trigger. The `affect!` function determines what happens to +the simulation when the callback is triggered. It takes a single object, an +`integrator` and makes some modification to it (please read more about +integrators [here](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/)). +Here, we access the system's unknown `X` through the `integrator`, and add +`5.0` to its amount. We can now simulate our system using the +callback: +```@example ex2 +sol = solve(oprob; callback = ps_cb) +plot(sol) +``` + +Next, we can also use a callback to change the parameters of a system. The following code +plots the concentration of a two-state system, as we change the equilibrium +constant between the two states: +```@example ex2 +rn = @reaction_network begin + (k,1), X1 <--> X2 +end +u0 = [:X1 => 10.0, :X2 => 0.0] +tspan = (0.0, 20.0) +p = [:k => 1.0] +oprob = ODEProblem(rn, u0, tspan, p) + +condition = [5.0] +affect!(integrator) = integrator[:k] = 5.0 +ps_cb = PresetTimeCallback(condition, affect!) + +sol = solve(oprob; callback = ps_cb) +plot(sol) +``` +The result looks as expected. However, what happens if we attempt to run the +simulation again? +```@example ex2 +sol = solve(oprob; callback = ps_cb) +plot(sol) +``` +The plot looks different, even though we simulate the same problem. Furthermore, +the callback does not seem to have any effect on the system. If we check our +`ODEProblem` +```@example ex2 +oprob.p +``` +we note that `k = 5.0`, rather than `k = 1.0` as we initially specified. This is +because the callback modifies our `ODEProblem` during the simulation, and this +modification remains during the second simulation. An improved workflow to avoid +this issue is: +```@example ex2 +rn = @reaction_network begin + (k,1), X1 <--> X2 +end +u0 = [:X1 => 10.0,:X2 => 0.0] +tspan = (0.0, 20.0) +p = [:k => 1.0] +oprob = ODEProblem(rn, u0, tspan, p) + +condition = [5.0] +affect!(integrator) = integrator[:k] = 5.0 +ps_cb = PresetTimeCallback(condition, affect!) + +sol = solve(deepcopy(oprob); callback = ps_cb) +plot(sol) +``` +where we parse a copy of our `ODEProblem` to the solver (using `deepcopy`). We can now run +```@example ex2 +sol = solve(deepcopy(oprob); callback = ps_cb) +plot(sol) +``` +and get the expected result. + +It is possible to give several callbacks to the `solve()` command. To do so, one +has to bundle them together in a `CallbackSet`, here follows one example: +```@example ex2 +rn = @reaction_network begin + (k,1), X1 <--> X2 +end +u0 = [:X1 => 10.0,:X2 => 0.0] +tspan = (0.0, 20.0) +p = [:k => 1.0] +oprob = ODEProblem(rn, u0, tspan, p) + +ps_cb_1 = PresetTimeCallback([3.0, 7.0], integ -> integ[:X1] += 5.0) +ps_cb_2 = PresetTimeCallback([5.0], integ -> integ[:k] = 5.0) + +sol = solve(deepcopy(oprob); callback=CallbackSet(ps_cb_1, ps_cb_2)) +plot(sol) +``` + +The difference between the `PresetTimeCallback`s and the `DiscreteCallback`s and +`ContiniousCallback`s is that the latter two allow the condition to be a +function, permitting the user to give more general conditions for the callback +to be triggered. An example could be a callback that triggers whenever a species +surpasses some threshold value. + +### [Callbacks during SSA simulations](@id advanced_simulations_ssa_callbacks) +An assumption of (most) SSA simulations is that the state of the system is unchanged between reaction events. However, callbacks that affect the system's state can violate this assumption. To prevent erroneous simulations, users must inform a SSA solver when the state has been updated in a callback. This allows the solver to reinitialize any internal state information that may have changed. This can be done through the `reset_aggregated_jumps!` function, see the following example: + +```@example ex2 +rn = @reaction_network begin + (k,1), X1 <--> X2 +end +u0 = [:X1 => 10.0,:X2 => 0.0] +tspan = (0.0, 20.0) +p = [:k => 1.0] +dprob = DiscreteProblem(rn, u0, tspan, p) +jprob = JumpProblem(rn, dprob, Direct()) + +condition = [5.0] +function affect!(integrator) + integrator[:X1] += 5.0 + integrator[:k] += 2.0 + reset_aggregated_jumps!(integrator) + nothing +end +cb = PresetTimeCallback(condition, affect!) + +sol = solve(deepcopy(jprob), SSAStepper(); callback=cb) +plot(sol) +``` \ No newline at end of file From 8e898334e0a69d06dae4c2ea14591c2c199d935c Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 30 Apr 2024 15:52:08 -0400 Subject: [PATCH 2/3] save progress --- .../events_and_callbacks.md | 91 +++++++++---------- 1 file changed, 41 insertions(+), 50 deletions(-) diff --git a/docs/src/catalyst_functionality/events_and_callbacks.md b/docs/src/catalyst_functionality/events_and_callbacks.md index 09fd8b651d..119b7d6ba0 100644 --- a/docs/src/catalyst_functionality/events_and_callbacks.md +++ b/docs/src/catalyst_functionality/events_and_callbacks.md @@ -3,52 +3,48 @@ Sometimes one wishes to model events, describing things that can happen to a system during its simulation. Examples can include: - A chemical system where an amount of some species is added at a time point after the simulation's initiation. - A simulation of a circadian rhythm, where light is turned on/off every 12 hours. - - A cell divides when some size variable reaches a certain threshold, halving the amount of each species in the system. + - A cell divides when some size variable reaches a certain threshold, halving its size. Generally, there are two approaches for creating events: -1. Including them as a part of the model when its is created. +1. Including them as a part of the model when it is created. 2. Creating a separate *callback*, which is then supplied to the `solve` function. Generally, the first approach is more convenient, and should be used if possible. However, the second approach is more potent, and can be used to create more general types of events. Below, we will describe both approaches. -## [Creating models with incorporated events](@id events_model) -Sometimes one wishes to model events, describing things that can happen to it during a simulation. - - A chemical system where an amount of some species is added at a time point after the simulation's initiation. - - A simulation of a circadian rhythm, where light is turned on/off every 12 hours. - - A cell divides when some size variable reaches a certain threshold, halving the amount of each species in the system. +!!! warn + Independently of how they are created, there are some additional considerations when creating events. These are described [at the end of this tutorial](@ref ref). Please briefly check these cases, and if any apply to you, read the corresponding section thoroughly. -Events are divided into *continuous* and *discrete* events, and these can be added directly to a system using the `continuous_events` and `discrete_events` options. Events can also be modelled through *callbacks*. These are different in that they are supplied in the simulation step (rather than on system creation), and generally provide more flexibility in how they may affect the system. Callbacks are described on a separate [page](@ref advanced_simulations_callbacks). - -The notation described below for creating continuous and discrete events is the same which is used in [ModelingToolkit to create events](https://docs.sciml.ai/ModelingToolkit/stable/basics/Events/), and which is used for [events for programmatic model creation](@ref ref). +## [Creating models with incorporated events](@id events_model) +Catalyst's DSL provides two options, `@continuous_events` and `@discrete_events`, which permits users to add events to the models when they are constructed. Each corresponds to the creation of a specific type of events, and both are described both. The syntax for creating evens in the DSL is identical to that used for [programmatic] model creation (described in detail [here]) and which is used by the [ModelingToolkit.jl package](https://docs.sciml.ai/ModelingToolkit/stable/basics/Events/) -### [Continuous vs discrete events](@id events_model_continuous_vs_discrete) -Both continuous and discrete events combine some condition (for triggering the event) with some affect (describing their effects on the system). They differ in the following ways: +Both continuous and discrete events combine a *condition* (for triggering the event) with an *affect* (describing its effect on the system). They differ in the following ways: - They use slightly different notation. - Discrete events' conditions are checked at *the end of* each simulation time step. For continuous events, the simulation instead finds the *exact time point* when the event is triggered at. - Continuous events cannot be supplied to jump simulations. ### [Continuous events](@id events_model_continuous) -Let us consider a simple system where species `X` degraded at a constant rate `d`. Next, we wish to add an event which adds `2.0` units of `X` whenever `X` reaches a critical threshold `1.0`. This can be done in the following manner: +Let us consider a simple system where species `X` degrades at a constant rate `d`. We wish to create an event that adds `2.0` units of `X` whenever `X` reaches a critical threshold `1.0`. This can be done in the following manner: ```@example dsl_advanced_events using Catalyst # hide rn = @reaction_network begin @continuous_events begin - X ~ 1.0 => [X ~ X + 2.0] + [X ~ 1.0] => [X ~ X + 2.0] end d, X --> 0 end nothing # hide ``` -Here, the `@continuous_events` option is followed by a `begin ... end` block. Next, each line corresponds to a separate event. Each event is created in the following manner: -- It combines a *condition* (denoting when the event will happen) with one (or several) *affects* (denoting what the event does to the system when it happens). -- The condition (left) and affect (right) are separated by a `=>`. -- The condition takes the form of an [equation](). Here, the event is triggered whenever the equation's two sides (separated by a `~`) are equal. -- The affect(s) are enclosed within `[]`. If there are multiple affects, these are separated by `,` (the example above contains a single affect). -- Each affect is a single equation that describes how a parameter, species, or [variable](@ref dsl_advanced_options_variables) is updated when the event is triggered. -- Each affect's equation's left-hand side must contain only the parameter/species/variable whose value should be updated. +Here, the `@continuous_events` option is followed by a `begin ... end` block. Next, each line corresponds to a separate event. Each continuous event is created in the following manner: +- It combines a set of *conditions* (describing when the event is triggered) with a set of *affects* (describing the event's effect on the system once triggered). +- The conditions and affects are vectors (containing each containing any number of condition or affect equations). +- The conditions (left) and affects (right) are separated by `=>`. +- Each condition is an [equation](@ref ref). Here, the event is triggered whenever any of its conditions equations' hold (i.e. the equation's left and right ahnd sides are identical). +- Each affect is a single [equation](@ref ref) that describes how a parameter or species (or [variable](@ref ref)) is updated when the event is triggered. +- Each affect's equation's left-hand side must contain only the parameter/species(/variable) whose value should be updated. - Each affect's equation's right-hand side is an expression describing its updated value. +- A single event's affects can updates species (and variables) only, or parameters only (but not parameters and species/variables). -We can simulate the model we declared, just like any other model: +The declared model can be simulated using standard syntax: ```@example dsl_advanced_events using OrdinaryDiffEq, Plots u0 = [:X => 2.0] @@ -66,13 +62,13 @@ In our example, we can also denote the critical quantities using parameters: rn = @reaction_network begin @parameters X_thres X_add @continuous_events begin - X ~ X_thres => [X ~ X + X_add] + [X ~ X_thres] => [X ~ X + X_add] end d, X --> 0 end nothing # hide ``` -Here, since `X_thres` and `X_add` do not appear in any reactions, Catalyst cannot determine whether they are parameters or species. hence, they must be [explicitly designated as parameters by using the `@parameters` option](@ref dsl_advanced_options_declaring_species_and_parameters). Next, these can be designated as any value, and supplied to the `ODEProblem`: +Here, since `X_thres` and `X_add` do not appear in any reactions, Catalyst cannot determine whether they are parameters or species. Hence, they must be [explicitly designated as parameters by using the `@parameters` option](@ref dsl_advanced_options_declaring_species_and_parameters). These values can now be designated when the `ODEProblem` is created: ```@example dsl_advanced_events u0 = [:X => 2.0] tspan = (0.0, 10.0) @@ -87,7 +83,7 @@ As previously noted, each continuous event can have multiple affects. The follow ```@example dsl_advanced_events rn = @reaction_network begin @continuous_events begin - X ~ Y => [X ~ X - 1.0, Y ~ Y + 1.0] + [X ~ Y] => [X ~ X - 1.0, Y ~ Y + 1.0] end p, 0 --> X d, Y --> 0 @@ -102,9 +98,6 @@ sol = solve(ODEProblem, Tsit5()) plot(sol) ``` -!!!warn - A single event (continuous or discrete) can update the value of either (one or several) species (and variables), or of (one or several) parameters. It is not possible for an event to update the values of both species/variables and parameters. - In the above examples we have modelled a system with a single event. In these cases, the `begin end` block is not required, and the event can be provided on the same line as the `@continuous_events` option: ```@example dsl_advanced_events rn = @reaction_network begin @@ -116,13 +109,13 @@ nothing # hide ``` ### [Discrete events](@id events_model_discrete) -Just like [continuous events](dsl_advanced_options_events_continuous), discrete events combine a condition with one or more affect statements. Here, discrete events' affects are created identically to those for continuous events. Discrete events' conditions are different. There exist 3 different types of discrete events, each with a different type of condition. All three types are created using the `@discrete_events` option, and a single system can contain a mix of all types. The three types are: +Discrete events are similar to continuous events in that they combine a condition for being triggered with an affect once triggered. The affect(s) of discrete events are declared, and works, identical to [those of continuous events](@ref events_model_continuous). However, the condition is different. There exist 3 different types of discrete events, each with a different type of condition. All three types are created using the `@discrete_events` option, and a single system can contain a mix of all types. The three types are: - Preset-time discrete events. - Periodic discrete events. - Conditional discrete events. #### [Preset-time discrete events](@id events_model_discrete_presettime) -*Present-time events* are events that happen at specific time points. Here, the condition is a vector with all the time points at which an event is triggered. E.g. here we create a production/degradation loop, where `2.0` units of `X` is added at time points `3.0` and `7.0` +*Present-time events* are events that happen at specific time points. Their conditions are vectors with all the time points at which the event is triggered. E.g. here we create a production/degradation loop, where `2.0` units of `X` is added at time points `3.0` and `7.0` ```@example dsl_advanced_events rn = @reaction_network begin @discrete_events begin @@ -140,17 +133,7 @@ sol = solve(ODEProblem, Tsit5()) plot(sol) ``` -The preset time points can also be parameters (in which case, they have to be [designated as such using the `@parameters` option](@ref dsl_advanced_options_declaring_species_and_parameters)): -```@example dsl_advanced_events -rn = @reaction_network begin - @parameters t1 t2 - @discrete_events begin - [t1, t2] => [X ~ X + 2.0] - end - (p,d), 0 <--> X -end -nothing -``` +Currently, the time points of preset cannot be parameters. If you want to create preset time event with a callback, either use a [continuous callback with condition `[t ~ preset_t]`](@ref events_model_continuous) (note, not possible for jump simulations), or [create a conditional discrete event and supply a `tstops` to the `solve` command](@ref events_additional_considerations_tstops) #### [Periodic discrete events](@id events_model_discrete_periodic) When a discrete event's condition is a vector, a preset-time event is created. If it instead is a single value, a *periodic event* is created. These occur repeatedly throughout a simulation, with its period set by the affect term. E.g. here we create a system where `0.5` units of `X` is added every `1.0` time units. @@ -170,10 +153,10 @@ oprob = ODEProblem(rn, u0, tspan, ps) sol = solve(ODEProblem, Tsit5()) plot(sol) ``` -Like for preset-time events, periodic events' affects may contain parameters. +Like for preset-time events, periodic events' conditions cannot contain parameters. #### [Conditional discrete events](@id events_model_discrete_conditional) -Finally, discrete events' condition may be a boolean expression (consisting of parameters, species, variables, and the time variable). Let's say that we want to create an event which, if the concentration of `X` is below a threshold `1.0`, adds `1.0` units of `X` to the system, then we can use the following discrete event: +Finally, discrete events' condition may be a boolean expression (consisting of any combinations of parameters, species, variables, numbers, and the time variable). Let's say that we want to create an event which, if the concentration of `X` is below a threshold `1.0`, adds `1.0` units of `X` to the system, then we can use the following discrete event: ```@example dsl_advanced_events rn = @reaction_network begin @discrete_events begin @@ -182,7 +165,7 @@ rn = @reaction_network begin d, X --> 0 end ``` -If we simulate the system using the same conditions as for our [similar, continuous, example](@ref dsl_advanced_options_events_continuous) we get the same result: +If we simulate the system using the same conditions as for our [similar, continuous, example](@ref dsl_advanced_options_events_continuous), the result is similar: ```@example dsl_advanced_events u0 = [:X => 2.0] tspan = (0.0, 10.0) @@ -193,10 +176,10 @@ sol = solve(ODEProblem, Tsit5()) plot(sol) ``` So, how is modelling this event as a discrete or continuous event different? There are four differences: -1) For continuous events, the simulation method finds the exact time point when the condition triggers. Discrete events are triggered at the first time step when the condition holds. +1) For continuous events, the simulation method finds the exact time point when the condition triggers. Discrete events are triggered after a time step after which the condition holds. 2) This discrete event will be triggered whenever `X < 1.0` holds, not just when the concentration of `X` passes the threshold. E.g. it will be triggered if the initial concentration of `X` is less than `1.0`. -3) Only the discrete event event can be used with jump simulations. -4) The discrete event can be used to create more advanced conditions. +3) Only discrete event events can be used with jump simulations. +4) Discrete event can be used to create more advanced conditions. E.g. using (4), we can modify our system so that the event is only triggered when time is less than `5.0` (after which `X` decays towards `0`): ```@example dsl_advanced_events @@ -215,10 +198,10 @@ oprob = ODEProblem(rn, u0, tspan, ps) sol = solve(ODEProblem, Tsit5()) plot(sol) ``` -!!!note - When we form composite boolean conditions for conditional discrete events, we use `&` to denote the AND operator (not `&&`, as this is currently not supported). +!!! note + When we form composite boolean conditions for conditional discrete events, we use `&` to denote the AND operator (not `&&`, as this is currently not supported). Similarly, `|` is sued instead of `||`. -!!!warn +!!! warn Generally, discrete events including equality (`==`) will not be triggered. The reason is that the condition is only checked at the end of every time step. Hence, unless special precautions are taken to ensure that the simulator stops when the condition holds, it will unlikely be triggered. @@ -378,6 +361,14 @@ function, permitting the user to give more general conditions for the callback to be triggered. An example could be a callback that triggers whenever a species surpasses some threshold value. +## [Additional consideration when using events (or callbacks)](@id events_additional_considerations) + +### [Events affecting changing model parameter values](@id events_additional_considerations_parameters) + +### [Events affecting changing model parameter values during jump simulations](@id events_additional_considerations_jump_parameters) + +### [Events occurring at specific time points](@id events_additional_considerations_tstops) + ### [Callbacks during SSA simulations](@id advanced_simulations_ssa_callbacks) An assumption of (most) SSA simulations is that the state of the system is unchanged between reaction events. However, callbacks that affect the system's state can violate this assumption. To prevent erroneous simulations, users must inform a SSA solver when the state has been updated in a callback. This allows the solver to reinitialize any internal state information that may have changed. This can be done through the `reset_aggregated_jumps!` function, see the following example: From 05c0aa5ef188d5ff7c9e979d531574f972446ff6 Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 9 May 2024 13:20:46 -0400 Subject: [PATCH 3/3] move to new doc structure --- docs/pages.jl | 2 +- .../05_events_and_callbacks.md} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename docs/src/{catalyst_functionality/events_and_callbacks.md => 02_model_creation/05_events_and_callbacks.md} (100%) diff --git a/docs/pages.jl b/docs/pages.jl index ee1f4ac113..35d8ad8485 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -10,7 +10,7 @@ pages = Any[ # DSL - Advanced "02_model_creation/03_programmatic_CRN_construction.md", "02_model_creation/04_constraint_equations.md", - # Events. + "02_model_creation/05_events_and_callbacks.md", "02_model_creation/06_compositional_modeling.md", # Distributed parameters, rates, and initial conditions. # Loading and writing models to files. diff --git a/docs/src/catalyst_functionality/events_and_callbacks.md b/docs/src/02_model_creation/05_events_and_callbacks.md similarity index 100% rename from docs/src/catalyst_functionality/events_and_callbacks.md rename to docs/src/02_model_creation/05_events_and_callbacks.md