third pass

Author: TorkelE

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).
+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)