From de2ea4d579ecefd875ff6ba78f255ddb7e33d677 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 10 Jun 2024 19:09:30 -0400 Subject: [PATCH 1/5] init --- docs/pages.jl | 3 + .../examples/periodic_events_simulation.md | 95 +++++++++++++++++++ 2 files changed, 98 insertions(+) create mode 100644 docs/src/model_simulation/examples/periodic_events_simulation.md diff --git a/docs/pages.jl b/docs/pages.jl index a8772bd62b..1bc2ad6b05 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -28,6 +28,9 @@ pages = Any[ "model_simulation/simulation_structure_interfacing.md", "model_simulation/ensemble_simulations.md", "model_simulation/ode_simulation_performance.md" + "Model simulation examples" => Any[ + "model_simulation/examples/periodic_events_simulation.md" + ] ], "Steady state analysis" => Any[ "steady_state_functionality/homotopy_continuation.md", diff --git a/docs/src/model_simulation/examples/periodic_events_simulation.md b/docs/src/model_simulation/examples/periodic_events_simulation.md new file mode 100644 index 0000000000..8536b2eee4 --- /dev/null +++ b/docs/src/model_simulation/examples/periodic_events_simulation.md @@ -0,0 +1,95 @@ +# [Modelling a periodic event during ODE and jump simulations](@id periodic_event_simulation_example) +This tutorial will describe how to simulate systems with periodic events in ODE (straightforward) and jump (slightly less, but still fairly, straightforward) simulations (SDEs use identical syntax to ODEs). We will consider a model with a [circadian rhythm](https://en.wikipedia.org/wiki/Circadian_rhythm), where a parameter represents the level of light. While outdoor light varies smoothly, in experimental settings a lamp is often simply turned on/off every 12 hours. Here we will model this toggling of the light using a periodic event that is triggered every 12 hours. + +## [Modelling a circadian periodic event in an ODE simulation](@id periodic_event_simulation_example_ode) +We will consider a simple circadian model, consisting of a single protein ($X$), which is phosphorylated ($X \to Xᴾ$) in the presence of light ($l$). Here, the light parameter can either be $0$ (night) or $1$ (day). We can model this using a simple periodic event which switches the value of $l$ every 12 hours (here, `%` is the [remainder operator](https://docs.julialang.org/en/v1/manual/mathematical-operations/#Arithmetic-Operators)). +```@example periodic_event_example +using Catalyst +circadian_model = @reaction_network begin + @discrete_events begin + 12 => [l ~ (l + 1)%2] + end + (kₚ*l,kᵢ), X <--> Xᴾ +end +``` +We can now simulate this model, observing how a 24-hour cycle is reached +```@example periodic_event_example +using OrdinaryDiffEq, Plots +u0 = [:X => 150.0, :Xᴾ => 50.0] +ps = [:kₚ => 0.1, :kᵢ => 0.1, :l => 1.0] +oprob = ODEProblem(circadian_model, u0, (0.0, 100.0), ps) +sol = solve(oprob) +plot(sol) +``` + +## [Modelling a circadian periodic event in a jump simulation](@id periodic_event_simulation_example_ode) +While defining periodic events is easy for ODE and SDE simulations, due to how events are internally implemented, these cannot currently be used for jump simulations. Instead, we create a new model (and its corresponding `JumpProblem`), without the event: +```@example periodic_event_example +circadian_model = @reaction_network begin + (kₚ*l,kᵢ), X <--> Xᴾ +end + +using JumpProcesses +u0 = [:X => 150, :Xᴾ => 50] +ps = [:kₚ => 0.1, :kᵢ => 0.1, :l => 1.0] +dprob = DiscreteProblem(circadian_model, u0, (0.0, 100.0), ps) +jprob = JumpProblem(circadian_model, dprob, Direct()) +nothing # hide +``` +Next, we implement our event through a `DiscreteCallback`. It triggers every 12 time units, and when it does, it flips the value of `l`. Since we will interface with $l$'s value multiple times, we [create specific `getl` and `setl` functions](@ref simulation_structure_interfacing_functions) to do this (in practice, due to this model's small size, the performance impact is negliable). +```@example periodic_event_example +getl = ModelingToolkit.getp(jprob, :l) +setl = ModelingToolkit.setp(jprob, :l) +function condition(u, t, integrator) + (t % 12) == 0.0 +end +function affect!(integrator) + setl(integrator, (getl(integrator) + 1) % 2) + reset_aggregated_jumps!(integrator) +end +callback = DiscreteCallback(condition, affect!) +nothing # hide +``` +!!! danger + Whenever the integrator of a jump simulator is modified, the `reset_aggregated_jumps!` function must be called (with the integrator as its input) afterwards. Omitting to do so will result in erroneous simulations. + +Next, if we simulate our model, we note that the events do not seem to be triggered +```@example periodic_event_example +sol = solve(jprob, SSAStepper(); callback) +plot(sol) +``` +The reason is that discrete callbacks' conditions are only checked at the end of each simulation time step, and the probability that these will coincide with the callback's trigger times ($t = 12, 24, 36, ...$) is infinitely small. Hence, we must directly instruct our simulation to stop at these time points. We can do this using the `tstops` argument: +```@example periodic_event_example +tstops = 12.0:12.0:dprob.tspan[2] +sol = solve(jprob, SSAStepper(); callback, tstops) +plot(sol) +``` + +## [Plotting the light level](@id periodic_event_simulation_plotting_light) +Sometimes when simulating models with periodic parameters, one would like to plot the parameter's value across the simulation. To do this, we must turn our parameter $l$ to a *variable* (so that its value is recorded throughout the simulation): +```@example periodic_event_example +circadian_model = @reaction_network begin + @variables l(t) + @discrete_events begin + 12 => [l ~ (l + 1)%2] + end + (kₚ*l,kᵢ), X <--> Xᴾ +end +nothing # hide +``` +Next, we simulate our model like before (but providing $l$'s value as an initial condition): +```@example periodic_event_example +u0 = [:X => 150.0, :Xᴾ => 50.0, :l => 1.0] +ps = [:kₚ => 0.1, :kᵢ => 0.1] +oprob = ODEProblem(circadian_model, u0, (0.0, 100.0), ps) +sol = solve(oprob) +nothing # hide +``` +If we directly plot $l$'s value, it will be too small (compared to $X$ and $Xᴾ$ to be discernible). We instead [`@unpack` our variables](@ref dsl_advanced_options_symbolics_and_DSL_unpack), and then plot a re-scaled version: +```@example periodic_event_example +@unpack X, Xᴾ, l = circadian_model +plot(sol; idxs = [X, Xᴾ, 150*l], labels = ["X" "Xᴾ" "Light amplitude"]) +``` + +!!! note + If you wish to reproduce this in a jump simulation, remember to make appropriate modifications (like using `setu` instead of `setp`). \ No newline at end of file From a89c2c0afa6490a44c9836458294382a96b96de7 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 10 Jun 2024 20:13:25 -0400 Subject: [PATCH 2/5] pages update --- docs/pages.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/pages.jl b/docs/pages.jl index 1bc2ad6b05..9e5e1cdfa7 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -27,7 +27,7 @@ pages = Any[ "model_simulation/simulation_plotting.md", "model_simulation/simulation_structure_interfacing.md", "model_simulation/ensemble_simulations.md", - "model_simulation/ode_simulation_performance.md" + "model_simulation/ode_simulation_performance.md", "Model simulation examples" => Any[ "model_simulation/examples/periodic_events_simulation.md" ] From 532ebc56bce36cc856f56bc4c471f9aeb60a824c Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Wed, 17 Jul 2024 09:21:50 -0400 Subject: [PATCH 3/5] Update docs/src/model_simulation/examples/periodic_events_simulation.md Co-authored-by: Sam Isaacson --- .../src/model_simulation/examples/periodic_events_simulation.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/model_simulation/examples/periodic_events_simulation.md b/docs/src/model_simulation/examples/periodic_events_simulation.md index 8536b2eee4..753d0085fc 100644 --- a/docs/src/model_simulation/examples/periodic_events_simulation.md +++ b/docs/src/model_simulation/examples/periodic_events_simulation.md @@ -51,7 +51,7 @@ callback = DiscreteCallback(condition, affect!) nothing # hide ``` !!! danger - Whenever the integrator of a jump simulator is modified, the `reset_aggregated_jumps!` function must be called (with the integrator as its input) afterwards. Omitting to do so will result in erroneous simulations. + Whenever the integrator of a jump simulator is modified, the `reset_aggregated_jumps!` function must be called (with the integrator as its input) afterwards. Omitting to do so will result in simulations (incorrectly) using the old values of parameters/states internally instead of the updated values. Next, if we simulate our model, we note that the events do not seem to be triggered ```@example periodic_event_example From 759f44e05d0eec93bca7782976d2b1c38b91e534 Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Wed, 17 Jul 2024 09:21:55 -0400 Subject: [PATCH 4/5] Update docs/src/model_simulation/examples/periodic_events_simulation.md Co-authored-by: Sam Isaacson --- .../src/model_simulation/examples/periodic_events_simulation.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/model_simulation/examples/periodic_events_simulation.md b/docs/src/model_simulation/examples/periodic_events_simulation.md index 753d0085fc..8b81d0ee51 100644 --- a/docs/src/model_simulation/examples/periodic_events_simulation.md +++ b/docs/src/model_simulation/examples/periodic_events_simulation.md @@ -36,7 +36,7 @@ dprob = DiscreteProblem(circadian_model, u0, (0.0, 100.0), ps) jprob = JumpProblem(circadian_model, dprob, Direct()) nothing # hide ``` -Next, we implement our event through a `DiscreteCallback`. It triggers every 12 time units, and when it does, it flips the value of `l`. Since we will interface with $l$'s value multiple times, we [create specific `getl` and `setl` functions](@ref simulation_structure_interfacing_functions) to do this (in practice, due to this model's small size, the performance impact is negliable). +Next, we implement our event through a `DiscreteCallback`. It triggers every 12 time units, and when it does, it flips the value of `l`. Since we will interface with $l$'s value multiple times, we [create specific `getl` and `setl` functions](@ref simulation_structure_interfacing_functions) to do this (in practice, due to this model's small size, the performance increase is minimal). ```@example periodic_event_example getl = ModelingToolkit.getp(jprob, :l) setl = ModelingToolkit.setp(jprob, :l) From 59b854ea9e4d3a284260f8999bb686526597a0ec Mon Sep 17 00:00:00 2001 From: Torkel Date: Sun, 21 Jul 2024 15:31:56 -0400 Subject: [PATCH 5/5] up --- .../examples/periodic_events_simulation.md | 32 ++++++------------- 1 file changed, 10 insertions(+), 22 deletions(-) diff --git a/docs/src/model_simulation/examples/periodic_events_simulation.md b/docs/src/model_simulation/examples/periodic_events_simulation.md index 8b81d0ee51..ab171a68d3 100644 --- a/docs/src/model_simulation/examples/periodic_events_simulation.md +++ b/docs/src/model_simulation/examples/periodic_events_simulation.md @@ -23,10 +23,13 @@ plot(sol) ``` ## [Modelling a circadian periodic event in a jump simulation](@id periodic_event_simulation_example_ode) -While defining periodic events is easy for ODE and SDE simulations, due to how events are internally implemented, these cannot currently be used for jump simulations. Instead, we create a new model (and its corresponding `JumpProblem`), without the event: +While defining periodic events is easy for ODE and SDE simulations, due to how events are internally implemented, these cannot currently be used for jump simulations. Instead, there is a workaround which includes firs creating a [conditional discrete event](@ref constraint_equations_events) which is designed to trigger every 13 time units: ```@example periodic_event_example circadian_model = @reaction_network begin - (kₚ*l,kᵢ), X <--> Xᴾ + @discrete_events begin + (t % 12 == 0) => [l ~ (l + 1)%2] + end + (kₚ*l,kᵢ), X <--> Xᴾ end using JumpProcesses @@ -36,37 +39,22 @@ dprob = DiscreteProblem(circadian_model, u0, (0.0, 100.0), ps) jprob = JumpProblem(circadian_model, dprob, Direct()) nothing # hide ``` -Next, we implement our event through a `DiscreteCallback`. It triggers every 12 time units, and when it does, it flips the value of `l`. Since we will interface with $l$'s value multiple times, we [create specific `getl` and `setl` functions](@ref simulation_structure_interfacing_functions) to do this (in practice, due to this model's small size, the performance increase is minimal). -```@example periodic_event_example -getl = ModelingToolkit.getp(jprob, :l) -setl = ModelingToolkit.setp(jprob, :l) -function condition(u, t, integrator) - (t % 12) == 0.0 -end -function affect!(integrator) - setl(integrator, (getl(integrator) + 1) % 2) - reset_aggregated_jumps!(integrator) -end -callback = DiscreteCallback(condition, affect!) -nothing # hide -``` -!!! danger - Whenever the integrator of a jump simulator is modified, the `reset_aggregated_jumps!` function must be called (with the integrator as its input) afterwards. Omitting to do so will result in simulations (incorrectly) using the old values of parameters/states internally instead of the updated values. - Next, if we simulate our model, we note that the events do not seem to be triggered ```@example periodic_event_example -sol = solve(jprob, SSAStepper(); callback) +sol = solve(jprob, SSAStepper()) plot(sol) +plot(sol, density = 10000, fmt = :png) # hide ``` The reason is that discrete callbacks' conditions are only checked at the end of each simulation time step, and the probability that these will coincide with the callback's trigger times ($t = 12, 24, 36, ...$) is infinitely small. Hence, we must directly instruct our simulation to stop at these time points. We can do this using the `tstops` argument: ```@example periodic_event_example tstops = 12.0:12.0:dprob.tspan[2] -sol = solve(jprob, SSAStepper(); callback, tstops) +sol = solve(jprob, SSAStepper(); tstops) plot(sol) +plot(sol, density = 10000, fmt = :png) # hide ``` ## [Plotting the light level](@id periodic_event_simulation_plotting_light) -Sometimes when simulating models with periodic parameters, one would like to plot the parameter's value across the simulation. To do this, we must turn our parameter $l$ to a *variable* (so that its value is recorded throughout the simulation): +Sometimes when simulating models with periodic parameters, one would like to plot the parameter's value across the simulation. For this, there are two potential strategies. One includes creating a [*saving callback*](https://docs.sciml.ai/DiffEqCallbacks/stable/output_saving/#DiffEqCallbacks.SavingCallback). The other one, which we will demonstrate here, includes turning the parameter $l$ to a *variable* (so that its value is recorded throughout the simulation): ```@example periodic_event_example circadian_model = @reaction_network begin @variables l(t)