From dce78c4aa897aac1ff86a562d304e5c867ffe6f9 Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Tue, 1 Apr 2025 22:35:06 +0100 Subject: [PATCH 1/8] init --- docs/pages.jl | 1 + ...ctivation_time_distribution_measurement.md | 53 +++++++++++++++++++ 2 files changed, 54 insertions(+) create mode 100644 docs/src/model_simulation/examples/activation_time_distribution_measurement.md diff --git a/docs/pages.jl b/docs/pages.jl index a5aaeae785..06b5abc75c 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -34,6 +34,7 @@ pages = Any[ "model_simulation/sde_simulation_performance.md", "Examples" => Any[ "model_simulation/examples/periodic_events_simulation.md", + "model_simulation/examples/activation_time_distribution_measurement.md", "model_simulation/examples/interactive_brusselator_simulation.md" ] ], diff --git a/docs/src/model_simulation/examples/activation_time_distribution_measurement.md b/docs/src/model_simulation/examples/activation_time_distribution_measurement.md new file mode 100644 index 0000000000..3e89449b7c --- /dev/null +++ b/docs/src/model_simulation/examples/activation_time_distribution_measurement.md @@ -0,0 +1,53 @@ +# [Measuring the Distribution of System Activation Times](@id activation_time_distribution_measurement) +In this example we will consider a model which, while initially inactive, activates in response to an input. The model is *stochastic*, causing the activation times to be *random*. By combining events, callbacks, and stochastic ensemble simulations, we will measure the probability distribution of the activation times (so called[*first passage times](https://en.wikipedia.org/wiki/First-hitting-time_model)). + +Our model will be a version of the [simple self-activation loop](@ref basic_CRN_library_self_activation) (the ensemble simulations of we have [considered previously](@ref ensemble_simulations_monte_carlo)). Here, we will consider the activation threshold parameter ($K$) to be activated by an input (at time $t = 0$). Before the input, $K$ is very large (essentially keeping the system inactive). After the input, it is reduced to a lower value (which permits the activation of the system). We will model this using two additional parameters ($Kᵢ$ and $Kₐ$, describing the pre and post-activation values of $K$, respectively). Initially, $K$ will [default to](@ref dsl_advanced_options_default_vals) $Kᵢ$. Next, at the activation time ($t = 0$), a discrete event will change $K$'s value to $Kᵢ$. +```@example activation_time_distribution_measurement +sa_model = @reaction_network begin + @parameters Kᵢ Kₐ K = Kᵢ + @discrete_events [0.0] => [K ~ Kₐ] + v0 + hill(X,v,K,n), 0 --> X + deg, X --> 0 +end +``` +Next, to perform stochastic simulations of the system we will create an `SDEProblem`. Here, we will need to assign parameter values to $Kᵢ$ and $Kₐ$, but not to $K$ (as its value is controlled by its default and the event). Also note that we start the simulation at a time $t < 0$. This ensures that by the input time ($t = 0$), the system have (more or less) reached its (inactive state) steady state distribution. It also means that the activation time can be measured exactly as the system time (as this is the time from the input at $t = 0$). +```@example activation_time_distribution_measurement +u0 = [:X => 10.0] +tspan = (-200.0, 2000.0) +ps = [:v0 => 0.1, :v => 2.5, :Kᵢ => 1000.0, :Kₐ => 40.0, :n => 3.0, :deg => 0.01] +sprob = SDEProblem(sa_model, u0, tspan, ps) +nothing # hide +``` +We can now create a simple `EnsembleProblem` and perform an ensemble simulation (as described [here](@ref ensemble_simulations)). Please note that the system contain an event which modifies its parameters, hence we must add the `safetycopy = true` argument to `EnsembleProblem` (else, subsequent simulations would start with $K = Kₐ$). +```@example activation_time_distribution_measurement +using Plots, StochasticDiffEq +eprob = EnsembleProblem(sprob; safetycopy = true) +esol = solve(eprob, ImplicitEM(); trajectories = 10) +plot(esol) +``` +Here we see how, after the input time, the system (randomly) switches from the inactive state to the active one (several examples of this, bistability-based, activation have been studied in the literature, both in models and experiments). + +Next, we wish to measure the distribution of these activation times. First we will create a [callback](@ref https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/) which terminates the simulation once it have reached a threshold. This both ensures that we do not have to expend unnecessary computer time on the simulation after its activation, and also enables us to measure the activation time as the final time point of the simulation. Here we will use a [discrete callback](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#SciMLBase.DiscreteCallback) (for an ODE simulation, a [continuous callback](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#ContinuousCallback) would instead have been appropriate, however, these combine less well with stochastic models). From looking at the previous simulations, $X = 100$ seems like a good threshold between the inactive and active state, so we will use this as our activation threshold. We use the `terminate!` function to terminate the simulation once the activation threshold has been reached. +```@example activation_time_distribution_measurement +condition(u, t, integrator) = integrator[:X] > 100.0 +affect!(integrator) = terminate!(integrator) +callback = DiscreteCallback(condition, affect!) +nothing # hide +``` +Next, we will perform our ensemble simulation. By default, for each simulation, these save the full trajectory. Here, however, we are only interested in teh activation time. To do this, we can utilise an [*output function*](https://docs.sciml.ai/DiffEqDocs/dev/features/ensemble/#Building-a-Problem). This will be called at the end of each single simulation and determine what to save (it can also be used to potentially rerun individual simulations, however, we will not use this feature here). Here we create an output function which saves only the simulation final time point. We also make it throw a warning if the simulation reached the final time point (which indicates a simulation that never activated, the occurrences of such simulation would cause us to underestimate the activation times). Finally, just like previously, we must set `safetycopy = true`. +```@example activation_time_distribution_measurement +function output_func(sol, i) + (sol.t[end] == tspan[2]) && @warn "A simulation did not activate during the given time span." + return (sol.t[end], false) +end +eprob = EnsembleProblem(sprob; output_func, safetycopy = true) +esol = solve(eprob, STrapezoid(); trajectories = 250, callback) +nothing # hide +``` +Finally, we can plot the distribution of activation times directly. For this, we will use the StatsPlots.jl package's `density` function (essentially a smoothed histogram). The input to `density` is the activation times (which our output function will save to `esol.u`). + +```@example activation_time_distribution_measurement +using StatsPlots +density(esol.u; label = "Activation times") +``` +Here we that the activation times take some form of long tail distribution (for non-trivial models, it is generally not possible to identify the activation times as any known statistical distribution). \ No newline at end of file From 7d91fdc4ccea8b3a8d5f64b696ea657945ef0710 Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Tue, 1 Apr 2025 22:50:39 +0100 Subject: [PATCH 2/8] second pass --- ...ctivation_time_distribution_measurement.md | 24 +++++++++---------- 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/docs/src/model_simulation/examples/activation_time_distribution_measurement.md b/docs/src/model_simulation/examples/activation_time_distribution_measurement.md index 3e89449b7c..1066e5e5e2 100644 --- a/docs/src/model_simulation/examples/activation_time_distribution_measurement.md +++ b/docs/src/model_simulation/examples/activation_time_distribution_measurement.md @@ -1,16 +1,17 @@ # [Measuring the Distribution of System Activation Times](@id activation_time_distribution_measurement) -In this example we will consider a model which, while initially inactive, activates in response to an input. The model is *stochastic*, causing the activation times to be *random*. By combining events, callbacks, and stochastic ensemble simulations, we will measure the probability distribution of the activation times (so called[*first passage times](https://en.wikipedia.org/wiki/First-hitting-time_model)). +In this example we will consider a model which, while initially inactive, activates in response to an input. The model is *stochastic*, causing the activation times to be *random*. By combining events, callbacks, and stochastic ensemble simulations, we will measure the probability distribution of the activation times (so-called [*first passage times()](https://en.wikipedia.org/wiki/First-hitting-time_model)). -Our model will be a version of the [simple self-activation loop](@ref basic_CRN_library_self_activation) (the ensemble simulations of we have [considered previously](@ref ensemble_simulations_monte_carlo)). Here, we will consider the activation threshold parameter ($K$) to be activated by an input (at time $t = 0$). Before the input, $K$ is very large (essentially keeping the system inactive). After the input, it is reduced to a lower value (which permits the activation of the system). We will model this using two additional parameters ($Kᵢ$ and $Kₐ$, describing the pre and post-activation values of $K$, respectively). Initially, $K$ will [default to](@ref dsl_advanced_options_default_vals) $Kᵢ$. Next, at the activation time ($t = 0$), a discrete event will change $K$'s value to $Kᵢ$. +Our model will be a version of the [simple self-activation loop](@ref basic_CRN_library_self_activation) (the ensemble simulations of which we have [considered previously](@ref ensemble_simulations_monte_carlo)). Here, we will consider the activation threshold parameter ($K$) to be activated by an input (at an input time $t = 0$). Before the input, $K$ is very large (essentially keeping the system inactive). After the input, it is reduced to a lower value (which permits the system to activate). We will model this using two additional parameters ($Kᵢ$ and $Kₐ$, describing the pre and post-activation values of $K$, respectively). Initially, $K$ will [default to](@ref dsl_advanced_options_default_vals) $Kᵢ$. Next, at the input time ($t = 0$), an event will change $K$'s value to $Kᵢ$. ```@example activation_time_distribution_measurement +using Catalyst sa_model = @reaction_network begin - @parameters Kᵢ Kₐ K = Kᵢ + @parameters Kᵢ Kₐ K=Kᵢ @discrete_events [0.0] => [K ~ Kₐ] v0 + hill(X,v,K,n), 0 --> X deg, X --> 0 end ``` -Next, to perform stochastic simulations of the system we will create an `SDEProblem`. Here, we will need to assign parameter values to $Kᵢ$ and $Kₐ$, but not to $K$ (as its value is controlled by its default and the event). Also note that we start the simulation at a time $t < 0$. This ensures that by the input time ($t = 0$), the system have (more or less) reached its (inactive state) steady state distribution. It also means that the activation time can be measured exactly as the system time (as this is the time from the input at $t = 0$). +Next, to perform stochastic simulations of the system we will create an `SDEProblem`. Here, we will need to assign parameter values to $Kᵢ$ and $Kₐ$, but not to $K$ (as its value is controlled by its default and the event). Also note that we start the simulation at a time $t = -200 < 0$. This ensures that by the input time ($t = 0$), the system has (more or less) reached its (inactive) steady state distribution. It also means that the activation time can be measured exactly as the simulation time at the time of activation (as this will be the time from the input at $t = 0$). ```@example activation_time_distribution_measurement u0 = [:X => 10.0] tspan = (-200.0, 2000.0) @@ -18,7 +19,7 @@ ps = [:v0 => 0.1, :v => 2.5, :Kᵢ => 1000.0, :Kₐ => 40.0, :n => 3.0, :deg => sprob = SDEProblem(sa_model, u0, tspan, ps) nothing # hide ``` -We can now create a simple `EnsembleProblem` and perform an ensemble simulation (as described [here](@ref ensemble_simulations)). Please note that the system contain an event which modifies its parameters, hence we must add the `safetycopy = true` argument to `EnsembleProblem` (else, subsequent simulations would start with $K = Kₐ$). +We can now create a simple `EnsembleProblem` and perform an ensemble simulation (as described [here](@ref ensemble_simulations)). Please note that the system s an event which modifies its parameters, hence we must add the `safetycopy = true` argument to `EnsembleProblem` (else, subsequent simulations would start with $K = Kₐ$). ```@example activation_time_distribution_measurement using Plots, StochasticDiffEq eprob = EnsembleProblem(sprob; safetycopy = true) @@ -27,27 +28,26 @@ plot(esol) ``` Here we see how, after the input time, the system (randomly) switches from the inactive state to the active one (several examples of this, bistability-based, activation have been studied in the literature, both in models and experiments). -Next, we wish to measure the distribution of these activation times. First we will create a [callback](@ref https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/) which terminates the simulation once it have reached a threshold. This both ensures that we do not have to expend unnecessary computer time on the simulation after its activation, and also enables us to measure the activation time as the final time point of the simulation. Here we will use a [discrete callback](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#SciMLBase.DiscreteCallback) (for an ODE simulation, a [continuous callback](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#ContinuousCallback) would instead have been appropriate, however, these combine less well with stochastic models). From looking at the previous simulations, $X = 100$ seems like a good threshold between the inactive and active state, so we will use this as our activation threshold. We use the `terminate!` function to terminate the simulation once the activation threshold has been reached. +Next, we wish to measure the distribution of these activation times. First we will create a [callback](@ref https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/) which terminates the simulation once it has reached a threshold. This both ensures that we do not have to expend unnecessary computer time on the simulation after its activation, and also enables us to measure the activation time as the final time point of the simulation. Here we will use a [discrete callback](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#SciMLBase.DiscreteCallback) (for an ODE simulation, a [continuous callback](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#ContinuousCallback) would instead have been appropriate, however, these combine less well with stochastic models). By looking at the previous simulations, we determine $X = 100$ as a suitable activation threshold. We use the `terminate!` function to terminate the simulation once this value has been reached. ```@example activation_time_distribution_measurement condition(u, t, integrator) = integrator[:X] > 100.0 affect!(integrator) = terminate!(integrator) callback = DiscreteCallback(condition, affect!) nothing # hide ``` -Next, we will perform our ensemble simulation. By default, for each simulation, these save the full trajectory. Here, however, we are only interested in teh activation time. To do this, we can utilise an [*output function*](https://docs.sciml.ai/DiffEqDocs/dev/features/ensemble/#Building-a-Problem). This will be called at the end of each single simulation and determine what to save (it can also be used to potentially rerun individual simulations, however, we will not use this feature here). Here we create an output function which saves only the simulation final time point. We also make it throw a warning if the simulation reached the final time point (which indicates a simulation that never activated, the occurrences of such simulation would cause us to underestimate the activation times). Finally, just like previously, we must set `safetycopy = true`. +Next, we will perform our ensemble simulation. By default, for each simulation, these save the full trajectory. Here, however, we are only interested in the activation time. This enables us to utilise an [*output function*](https://docs.sciml.ai/DiffEqDocs/dev/features/ensemble/#Building-a-Problem). This will be called at the end of every single simulation and determine what to save (it can also be used to potentially rerun individual simulations, however, we will not use this feature here). Here we create an output function which saves only the simulation's final time point. We also make it throw a warning if the simulation reaches the end of the simulation time frame ($t = 2000$) (this indicates a simulation that never activated, the occurrences of such simulation would cause us to underestimate the activation times). Finally, just like previously, we must set `safetycopy = true`. ```@example activation_time_distribution_measurement function output_func(sol, i) (sol.t[end] == tspan[2]) && @warn "A simulation did not activate during the given time span." return (sol.t[end], false) end eprob = EnsembleProblem(sprob; output_func, safetycopy = true) -esol = solve(eprob, STrapezoid(); trajectories = 250, callback) +esol = solve(eprob, ImplicitEM(); trajectories = 250, callback) nothing # hide ``` -Finally, we can plot the distribution of activation times directly. For this, we will use the StatsPlots.jl package's `density` function (essentially a smoothed histogram). The input to `density` is the activation times (which our output function will save to `esol.u`). - +Finally, we can plot the distribution of activation times. For this, we will use the [StatsPlots.jl](https://docs.juliaplots.org/latest/generated/statsplots/) package's `density` function (essentially a smoothed histogram). The input to `density` is the activation times (which our output function has saved to `esol.u`). ```@example activation_time_distribution_measurement using StatsPlots -density(esol.u; label = "Activation times") +density(esol.u; label = "Activation time distribution", xlabel = "Activation time") ``` -Here we that the activation times take some form of long tail distribution (for non-trivial models, it is generally not possible to identify the activation times as any known statistical distribution). \ No newline at end of file +Here we that the activation times take some form of long tail distribution (for non-trivial models like this one, it is generally not possible to identify the activation times as any known statistical distribution). \ No newline at end of file From 7a23dee5b9100de82c0e084409fe49daff744db9 Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Tue, 1 Apr 2025 23:12:17 +0100 Subject: [PATCH 3/8] third pass --- .../activation_time_distribution_measurement.md | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/docs/src/model_simulation/examples/activation_time_distribution_measurement.md b/docs/src/model_simulation/examples/activation_time_distribution_measurement.md index 1066e5e5e2..54c6317896 100644 --- a/docs/src/model_simulation/examples/activation_time_distribution_measurement.md +++ b/docs/src/model_simulation/examples/activation_time_distribution_measurement.md @@ -1,5 +1,5 @@ # [Measuring the Distribution of System Activation Times](@id activation_time_distribution_measurement) -In this example we will consider a model which, while initially inactive, activates in response to an input. The model is *stochastic*, causing the activation times to be *random*. By combining events, callbacks, and stochastic ensemble simulations, we will measure the probability distribution of the activation times (so-called [*first passage times()](https://en.wikipedia.org/wiki/First-hitting-time_model)). +In this example we will consider a model which, while initially inactive, activates in response to an input. The model is *stochastic*, causing the activation times to be *random*. By combining events, callbacks, and stochastic ensemble simulations, we will measure the probability distribution of the activation times (so-called [*first passage times*](https://en.wikipedia.org/wiki/First-hitting-time_model)). Our model will be a version of the [simple self-activation loop](@ref basic_CRN_library_self_activation) (the ensemble simulations of which we have [considered previously](@ref ensemble_simulations_monte_carlo)). Here, we will consider the activation threshold parameter ($K$) to be activated by an input (at an input time $t = 0$). Before the input, $K$ is very large (essentially keeping the system inactive). After the input, it is reduced to a lower value (which permits the system to activate). We will model this using two additional parameters ($Kᵢ$ and $Kₐ$, describing the pre and post-activation values of $K$, respectively). Initially, $K$ will [default to](@ref dsl_advanced_options_default_vals) $Kᵢ$. Next, at the input time ($t = 0$), an event will change $K$'s value to $Kᵢ$. ```@example activation_time_distribution_measurement @@ -26,7 +26,7 @@ eprob = EnsembleProblem(sprob; safetycopy = true) esol = solve(eprob, ImplicitEM(); trajectories = 10) plot(esol) ``` -Here we see how, after the input time, the system (randomly) switches from the inactive state to the active one (several examples of this, bistability-based, activation have been studied in the literature, both in models and experiments). +Here we see how, after the input time, the system (randomly) switches from the inactive state to the active one (several examples of this, bistability-based, activation have been studied in the literature, both in models and experiments[^1][^2]). Next, we wish to measure the distribution of these activation times. First we will create a [callback](@ref https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/) which terminates the simulation once it has reached a threshold. This both ensures that we do not have to expend unnecessary computer time on the simulation after its activation, and also enables us to measure the activation time as the final time point of the simulation. Here we will use a [discrete callback](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#SciMLBase.DiscreteCallback) (for an ODE simulation, a [continuous callback](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#ContinuousCallback) would instead have been appropriate, however, these combine less well with stochastic models). By looking at the previous simulations, we determine $X = 100$ as a suitable activation threshold. We use the `terminate!` function to terminate the simulation once this value has been reached. ```@example activation_time_distribution_measurement @@ -50,4 +50,9 @@ Finally, we can plot the distribution of activation times. For this, we will use using StatsPlots density(esol.u; label = "Activation time distribution", xlabel = "Activation time") ``` -Here we that the activation times take some form of long tail distribution (for non-trivial models like this one, it is generally not possible to identify the activation times as any known statistical distribution). \ No newline at end of file +Here we that the activation times take some form of long tail distribution (for non-trivial models like this one, it is generally not possible to identify the activation times as any known statistical distribution). + +--- +## References +[^1]: [David Frigola, Laura Casanellas, José M. Sancho, Marta Ibañes, *Asymmetric Stochastic Switching Driven by Intrinsic Molecular Noise*, PLoS One (2012).](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0031407) +[^2]: [Christian P Schwall et al., *Tunable phenotypic variability through an autoregulatory alternative sigma factor circuit*, Molecular Systems Biology (2021).](https://www.embopress.org/doi/full/10.15252/msb.20209832) \ No newline at end of file From 6ef8c46f0684ca62cb3e971b4350dfd4dde50b95 Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Tue, 1 Apr 2025 23:27:06 +0100 Subject: [PATCH 4/8] add StatsPlots dependency to docs --- docs/Project.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/Project.toml b/docs/Project.toml index e153f3a71c..837d874798 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -40,6 +40,7 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f" StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544" @@ -86,6 +87,7 @@ SciMLBase = "2.46" SciMLSensitivity = "7.60" SpecialFunctions = "2.4" StaticArrays = "1.9" +StatsPlots = "0.15.7" SteadyStateDiffEq = "2.2" StochasticDiffEq = "6.65" StructuralIdentifiability = "0.5.11" From 971d73a42ed0884722f335786dcfd3642ba44bc6 Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Wed, 2 Apr 2025 07:28:52 +0100 Subject: [PATCH 5/8] cross reference fix --- .../examples/activation_time_distribution_measurement.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/model_simulation/examples/activation_time_distribution_measurement.md b/docs/src/model_simulation/examples/activation_time_distribution_measurement.md index 54c6317896..55cc1ebc14 100644 --- a/docs/src/model_simulation/examples/activation_time_distribution_measurement.md +++ b/docs/src/model_simulation/examples/activation_time_distribution_measurement.md @@ -28,7 +28,7 @@ plot(esol) ``` Here we see how, after the input time, the system (randomly) switches from the inactive state to the active one (several examples of this, bistability-based, activation have been studied in the literature, both in models and experiments[^1][^2]). -Next, we wish to measure the distribution of these activation times. First we will create a [callback](@ref https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/) which terminates the simulation once it has reached a threshold. This both ensures that we do not have to expend unnecessary computer time on the simulation after its activation, and also enables us to measure the activation time as the final time point of the simulation. Here we will use a [discrete callback](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#SciMLBase.DiscreteCallback) (for an ODE simulation, a [continuous callback](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#ContinuousCallback) would instead have been appropriate, however, these combine less well with stochastic models). By looking at the previous simulations, we determine $X = 100$ as a suitable activation threshold. We use the `terminate!` function to terminate the simulation once this value has been reached. +Next, we wish to measure the distribution of these activation times. First we will create a [callback](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/) which terminates the simulation once it has reached a threshold. This both ensures that we do not have to expend unnecessary computer time on the simulation after its activation, and also enables us to measure the activation time as the final time point of the simulation. Here we will use a [discrete callback](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#SciMLBase.DiscreteCallback). By looking at the previous simulations, we determine $X = 100$ as a suitable activation threshold. We use the `terminate!` function to terminate the simulation once this value has been reached. ```@example activation_time_distribution_measurement condition(u, t, integrator) = integrator[:X] > 100.0 affect!(integrator) = terminate!(integrator) From 5abb04841e79332155e18b14c53e342ecc484e65 Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Wed, 9 Apr 2025 10:52:28 +0200 Subject: [PATCH 6/8] Update docs/src/model_simulation/examples/activation_time_distribution_measurement.md Co-authored-by: Sam Isaacson --- .../examples/activation_time_distribution_measurement.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/model_simulation/examples/activation_time_distribution_measurement.md b/docs/src/model_simulation/examples/activation_time_distribution_measurement.md index 55cc1ebc14..39d5dc4324 100644 --- a/docs/src/model_simulation/examples/activation_time_distribution_measurement.md +++ b/docs/src/model_simulation/examples/activation_time_distribution_measurement.md @@ -19,7 +19,7 @@ ps = [:v0 => 0.1, :v => 2.5, :Kᵢ => 1000.0, :Kₐ => 40.0, :n => 3.0, :deg => sprob = SDEProblem(sa_model, u0, tspan, ps) nothing # hide ``` -We can now create a simple `EnsembleProblem` and perform an ensemble simulation (as described [here](@ref ensemble_simulations)). Please note that the system s an event which modifies its parameters, hence we must add the `safetycopy = true` argument to `EnsembleProblem` (else, subsequent simulations would start with $K = Kₐ$). +We can now create a simple `EnsembleProblem` and perform an ensemble simulation (as described [here](@ref ensemble_simulations)). Please note that the system has an event which modifies its parameters, hence we must add the `safetycopy = true` argument to `EnsembleProblem` (else, subsequent simulations would start with $K = Kₐ$). ```@example activation_time_distribution_measurement using Plots, StochasticDiffEq eprob = EnsembleProblem(sprob; safetycopy = true) From 7932ff427912ab54416b49fd518d2b34ccc03913 Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Wed, 9 Apr 2025 11:46:42 +0100 Subject: [PATCH 7/8] remove statsplots --- docs/Project.toml | 2 -- .../examples/activation_time_distribution_measurement.md | 5 ++--- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 837d874798..e153f3a71c 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -40,7 +40,6 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" -StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f" StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544" @@ -87,7 +86,6 @@ SciMLBase = "2.46" SciMLSensitivity = "7.60" SpecialFunctions = "2.4" StaticArrays = "1.9" -StatsPlots = "0.15.7" SteadyStateDiffEq = "2.2" StochasticDiffEq = "6.65" StructuralIdentifiability = "0.5.11" diff --git a/docs/src/model_simulation/examples/activation_time_distribution_measurement.md b/docs/src/model_simulation/examples/activation_time_distribution_measurement.md index 39d5dc4324..49a2e762bf 100644 --- a/docs/src/model_simulation/examples/activation_time_distribution_measurement.md +++ b/docs/src/model_simulation/examples/activation_time_distribution_measurement.md @@ -45,10 +45,9 @@ eprob = EnsembleProblem(sprob; output_func, safetycopy = true) esol = solve(eprob, ImplicitEM(); trajectories = 250, callback) nothing # hide ``` -Finally, we can plot the distribution of activation times. For this, we will use the [StatsPlots.jl](https://docs.juliaplots.org/latest/generated/statsplots/) package's `density` function (essentially a smoothed histogram). The input to `density` is the activation times (which our output function has saved to `esol.u`). +Finally, we can plot the distribution of activation times. For this, we will use the [`histogram`](@ref https://docs.juliaplots.org/latest/series_types/histogram/) function (with the `normalize = true` argument to create a probability density function). An alternative we also recommend is [StatsPlots.jl](https://docs.juliaplots.org/latest/generated/statsplots/)'s `density` function (which creates a smoothed histogram that is also easier to combine with other plots). The input to `density` is the activation times (which our output function has saved to `esol.u`). ```@example activation_time_distribution_measurement -using StatsPlots -density(esol.u; label = "Activation time distribution", xlabel = "Activation time") +histogram(esol.u; normalize = true, label = "Activation time distribution", xlabel = "Activation time") ``` Here we that the activation times take some form of long tail distribution (for non-trivial models like this one, it is generally not possible to identify the activation times as any known statistical distribution). From 766722ae59981974ed8d08b4fc5e7ef485a69153 Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Thu, 10 Apr 2025 22:27:04 +0200 Subject: [PATCH 8/8] Update activation_time_distribution_measurement.md --- .../examples/activation_time_distribution_measurement.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/model_simulation/examples/activation_time_distribution_measurement.md b/docs/src/model_simulation/examples/activation_time_distribution_measurement.md index 49a2e762bf..842c4ece35 100644 --- a/docs/src/model_simulation/examples/activation_time_distribution_measurement.md +++ b/docs/src/model_simulation/examples/activation_time_distribution_measurement.md @@ -45,7 +45,7 @@ eprob = EnsembleProblem(sprob; output_func, safetycopy = true) esol = solve(eprob, ImplicitEM(); trajectories = 250, callback) nothing # hide ``` -Finally, we can plot the distribution of activation times. For this, we will use the [`histogram`](@ref https://docs.juliaplots.org/latest/series_types/histogram/) function (with the `normalize = true` argument to create a probability density function). An alternative we also recommend is [StatsPlots.jl](https://docs.juliaplots.org/latest/generated/statsplots/)'s `density` function (which creates a smoothed histogram that is also easier to combine with other plots). The input to `density` is the activation times (which our output function has saved to `esol.u`). +Finally, we can plot the distribution of activation times. For this, we will use the [`histogram`](https://docs.juliaplots.org/latest/series_types/histogram/) function (with the `normalize = true` argument to create a probability density function). An alternative we also recommend is [StatsPlots.jl](https://docs.juliaplots.org/latest/generated/statsplots/)'s `density` function (which creates a smoothed histogram that is also easier to combine with other plots). The input to `density` is the activation times (which our output function has saved to `esol.u`). ```@example activation_time_distribution_measurement histogram(esol.u; normalize = true, label = "Activation time distribution", xlabel = "Activation time") ``` @@ -54,4 +54,4 @@ Here we that the activation times take some form of long tail distribution (for --- ## References [^1]: [David Frigola, Laura Casanellas, José M. Sancho, Marta Ibañes, *Asymmetric Stochastic Switching Driven by Intrinsic Molecular Noise*, PLoS One (2012).](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0031407) -[^2]: [Christian P Schwall et al., *Tunable phenotypic variability through an autoregulatory alternative sigma factor circuit*, Molecular Systems Biology (2021).](https://www.embopress.org/doi/full/10.15252/msb.20209832) \ No newline at end of file +[^2]: [Christian P Schwall et al., *Tunable phenotypic variability through an autoregulatory alternative sigma factor circuit*, Molecular Systems Biology (2021).](https://www.embopress.org/doi/full/10.15252/msb.20209832)