Skip to content

Commit d624a6b

Browse files
committed
second pass
1 parent 09576a6 commit d624a6b

File tree

2 files changed

+42
-38
lines changed

2 files changed

+42
-38
lines changed

docs/pages.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@ pages = Any[
2222
"Examples" => Any[
2323
"model_creation/examples/basic_CRN_library.md",
2424
"model_creation/examples/programmatic_generative_linear_pathway.md",
25+
"model_creation/examples/noise_modelling_approaches.md",
2526
"model_creation/examples/hodgkin_huxley_equation.md",
2627
"model_creation/examples/smoluchowski_coagulation_equation.md"
2728
]
Lines changed: 41 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -1,27 +1,27 @@
11
# [Approaches for modelling system noise](@id noise_modelling_approaches)
2-
Catalyst's primary tools for modelling stochasticity includes the creation of `SDEProblem`s or `JumpProblem`s from reaction network models. However, other approaches for incorporating noise in a model exists, some of which will be discussed here. We will first consider *intrinsic* and *extrinsic* noise. These are well-established terms, both of which we will describe in detail below (however, to our knowledge, no generally agreed upon definition of these exists)[^1]. The final approach, the utilisation of a noisy input process to an otherwise deterministic system, is less frequently used. However, as it has been used in the literature, we will demonstrate it here as well.
2+
Catalyst's primary tools for modelling stochasticity include the creation of `SDEProblem`s or `JumpProblem`s from reaction network models. However, other approaches for incorporating model noise exist, some of which will be discussed here. We will first consider *intrinsic* and *extrinsic* noise. These are well-established terms, both of which we will describe below (however, to our knowledge, no generally agreed-upon definition of these terms exists)[^1]. Finally, we will demonstrate a third approach, the utilisation of a noisy input process to an otherwise deterministic system. This approach is infrequently used, however, as it is encountered in the literature, we will demonstrate it here as well.
33

4-
Finally, we note that these approaches can all be combined. E.g. modelling intrinsic noise (using an SDE) can be combined extrinsic noise (using randomised parameter values) while also feeding a noisy input process into the system.
4+
We note that these approaches can all be combined. E.g. an intrinsic noise model (using an SDE) can be combined with extrinsic noise (using randomised parameter values), while also feeding a noisy input process into the system.
55

66
!!! note
7-
Here we use intrinsic and extrinsic noise as description of two of our modelling approaches. It should be noted that while these are established terminologies for noisy biological systems[^1], our use of these terms to describe different approaches for modelling noise is only inspired by this terminology, and nothing that is established in the field.
7+
Here we use intrinsic and extrinsic noise as descriptions of two of our modelling approaches. It should be noted that while these are established terminologies for noisy biological systems[^1], our use of these terms to describe different approaches for modelling noise is only inspired by this terminology, and nothing that is established in the field. Please consider the [references](@ref noise_modelling_approaches_references) for more information on intrinsic and extrinsic noise.
88

99
## [The repressilator model](@id noise_modelling_approaches_model_intro)
1010
For this tutorial we will use the oscillating [repressilator](@ref basic_CRN_library_repressilator) model.
1111
```@example noise_modelling_approaches
1212
using Catalyst
1313
repressilator = @reaction_network begin
14-
hillr(Z,v,K,n), ∅ --> X
15-
hillr(X,v,K,n), ∅ --> Y
16-
hillr(Y,v,K,n), ∅ --> Z
14+
    hillr(Z,v,K,n), ∅ --> X
15+
    hillr(X,v,K,n), ∅ --> Y
16+
    hillr(Y,v,K,n), ∅ --> Z
1717
d, (X, Y, Z) --> ∅
1818
end
1919
```
2020

2121
## [Using intrinsic noise](@id noise_modelling_approaches_model_intrinsic)
22-
Generally, intrinsic noise is randomness inherent to a system itself. This means that it cannot be controlled for, or filtered out by, experimental settings. Low-copy number cellular systems, were reaction occurs due to the encounters of molecules due to random diffusion, is an example of intrinsic noise. In practise, this can be precisely modelled through [SDE](@ref simulation_intro_SDEs) (chemical Langevin equations) or [jump](@ref simulation_intro_jumps) (stochastic chemical kinetics) simulations.
22+
Generally, intrinsic noise is randomness inherent to a system itself. This means that it cannot be controlled for, or filtered out by, experimental settings. Low-copy number cellular systems, were reaction occurs due to the encounters of molecules due to random diffusion, is an example of intrinsic noise. In practise, this can be modelled exactly through [SDE](@ref simulation_intro_SDEs) (chemical Langevin equations) or [jump](@ref simulation_intro_jumps) (stochastic chemical kinetics) simulations.
2323

24-
In Catalyst, intrinsic noise is accounted for whenever a `SDEProblem` or `JumpProblem` is created and simulated. Here we will model intrinsic noise through SDEs, which means creating a `SDEProblem` using the standard approach.
24+
In Catalyst, intrinsic noise is accounted for whenever an `SDEProblem` or `JumpProblem` is created and simulated. Here we will model intrinsic noise through SDEs, which means creating an `SDEProblem` using the standard approach.
2525
```@example noise_modelling_approaches
2626
u0 = [:X => 45.0, :Y => 20.0, :Z => 20.0]
2727
tend = 200.0
@@ -39,66 +39,64 @@ plot(sol_intrinsic; idxs = :X)
3939
Here, each simulation is performed from the same system using the same settings. Despite this, due to the noise, the individual trajectories are different.
4040

4141
## [Using extrinsic noise](@id noise_modelling_approaches_model_extrinsic)
42-
Next, we consider extrinsic noise. This is randomness caused by stochasticity external to, yet affecting, a system. Examples could be different bacteria experiencing different microenvironments or cells being in different part of the cell cycle. Here, what noise is intrinsic and extrinsic to a system may depend on how one defines the system itself (which is a reason why exact definitions of these terms is difficult, please consider the [references](@ref noise_modelling_approaches_references) for more information).
42+
Next, we consider extrinsic noise. This is randomness caused by stochasticity external to, yet affecting, a system. Examples could be different bacteria experiencing different microenvironments or cells being in different parts of the cell cycle. This is noise which (in theory) can be controlled for experimentally (e.g. by ensuring a uniform environment). Whenever a specific source of noise is intrinsic and extrinsic to a system may depend on how one defines the system itself (this is a reason why giving an exact definition of these terms is difficult).
4343

44-
In Catalyst we can modelled extrinsic noise by letting the model parameters be probability distributions. Here, at the beginning of each simulation, random parameter values are drawn from their distributions. Let us imagine that our repressilator circuit was inserted into an bacterial population. Here, while each bacteria would have the same circuit, their individual number of e.g. ribosomes (which will be random) might affect the production rates (which while constant within each bacteria, might differ between the individuals).
44+
In Catalyst we can model extrinsic noise by letting the model parameters be probability distributions. Here, at the beginning of each simulation, random parameter values are drawn from their distributions. Let us imagine that our repressilator circuit was inserted into a bacterial population. Here, while each bacteria would have the same circuit, their individual number of e.g. ribosomes (which will be random) might affect the production rates (which while constant within each bacteria, might differ between the individuals).
4545

46-
Again we will perform ensemble simulation. Instead of creating an `SDEProblem`, we will create an `ODEProblem`, as well as a [problem function](@ref ensemble_simulations_varying_conditions) which draws random parameter values for each simulations (here we have implemented the parameter's probability distributions using the [Distributions.jl](https://github.com/JuliaStats/Distributions.jl) package).
46+
Again we will perform ensemble simulation. Instead of creating an `SDEProblem`, we will create an `ODEProblem`, as well as a [problem function](@ref ensemble_simulations_varying_conditions) which draws random parameter values for each simulation. Here we have implemented the parameter's probability distributions as [normal distributions](https://en.wikipedia.org/wiki/Normal_distribution) using the [Distributions.jl](https://github.com/JuliaStats/Distributions.jl) package.
4747
```@example noise_modelling_approaches
4848
using Distributions
49-
p_dists = Dict([:v => Normal(10.0, 2.0), :K => Normal(20.0, 5.0),
50-
:n => Normal(3, 0.2), :d => Normal(0.1, 0.02)])
49+
p_dists = Dict([:v => Normal(10.0, 2.0), :K => Normal(20.0, 5.0), :n => Normal(3, 0.2), :d => Normal(0.1, 0.02)])
5150
function prob_func(prob, i, repeat)
5251
p = [par => rand(p_dists[par]) for par in keys(p_dists)]
5352
return remake(prob; p)
5453
end
5554
nothing # hide
5655
```
57-
Next, we again performs 4 simulations. While the individual trajectories are performed using deterministic simulations, the randomised parameter values creates heterogeneity across the ensemble.
56+
Next, we again perform 4 simulations. While the individual trajectories are performed using deterministic simulations, the randomised parameter values create heterogeneity across the ensemble.
5857
```@example noise_modelling_approaches
5958
using OrdinaryDiffEqDefault
6059
oprob = ODEProblem(repressilator, u0, tend, ps)
6160
eprob_extrinsic = EnsembleProblem(oprob; prob_func)
6261
sol_extrinsic = solve(eprob_extrinsic; trajectories = 4)
6362
plot(sol_extrinsic; idxs = :X)
6463
```
65-
We note that a similar approach can be used to also randomise the initial conditions. Finally, in a very detailed model the parameter values could fluctuate across the simulations, something which could be implemented using the approach from the next section.
64+
We note that a similar approach can be used to also randomise the initial conditions. In a very detailed model, the parameter values could fluctuate during a single simulation, something which could be implemented using the approach from the next section.
6665

6766
## [Using a noisy input process](@id noise_modelling_approaches_model_input_noise)
68-
Finally, we will consider the case where we have a deterministic system, but which is exposed to a some noisy input process. One example could be a [light sensitive system, where the amount of experienced sunlight is stochastic due to e.g. variable cloud cover](@ref functional_parameters_circ_rhythm). Practically, this can be considered extrinsic noise, but the model is created using a different approach from teh previous section. Here, we will pre-simulate a random process in time, which we then feed into the system as a function, time-dependent, parameter. The technical aspects of this approach is described in more details [here](@ref time_dependent_parameters).
67+
Finally, we will consider the case where we have a deterministic system, but which is exposed to a noisy input process. One example could be a [light sensitive system, where the amount of experienced sunlight is stochastic due to e.g. variable cloud cover](@ref functional_parameters_circ_rhythm). Practically, this can be considered as extrinsic noise, however, we will generate the noise using a different approach from in the previous section. Here, we pre-simulate a random process in time, which we then feed into the system as a functional, time-dependent, parameter. A more detailed description of functional parameters can be found [here](@ref time_dependent_parameters).
6968

70-
He we assume that our repressilator have an input, which correspond to the $K$ value controlling $X$'s production. First we create a function, `make_K_series`, which creates a randomised time series representing $K$'s value over time.
69+
We assume that our repressilator has an input, which corresponds to the $K$ value that controls $X$'s production. First we create a function, `make_K_series`, which creates a randomised time series representing $K$'s value over time.
7170
```@example noise_modelling_approaches
7271
using DataInterpolations
73-
function make_K_series(;K_mean = 20.0, n = 500, θ = 0.01)
72+
function make_K_series(; K_mean = 20.0, n = 500, θ = 0.01)
7473
t_samples = range(0.0, stop = tend, length = n)
7574
K_series = fill(K_mean, n)
7675
for i = 2:n
77-
K_series[i] = K_series[i-1] + (rand()-0.5) - θ*(K_series[i-1] - K_mean)
78-
end
79-
return LinearInterpolation(K_series, t_samples)
76+
K_series[i] = K_series[i-1] + (rand() - 0.5) - θ*(K_series[i-1] - K_mean)
77+
    end
78+
    return LinearInterpolation(K_series, t_samples)
8079
end
8180
plot(make_K_series())
8281
```
8382
Next, we create an updated repressilator model, where the input $K$ value is modelled as a time-dependent parameter.
8483
```@example noise_modelling_approaches
8584
@parameters (K_in::typeof(make_K_series()))(..)
8685
K_in = K_in(default_t())
87-
repressilator = @reaction_network begin
88-
hillr(Z,v,$K_in,n), ∅ --> X
89-
hillr(X,v,K,n), ∅ --> Y
90-
hillr(Y,v,K,n), ∅ --> Z
86+
repressilator_Kin = @reaction_network begin
87+
    hillr(Z,v,$K_in,n), ∅ --> X
88+
    hillr(X,v,K,n), ∅ --> Y
89+
    hillr(Y,v,K,n), ∅ --> Z
9190
d, (X, Y, Z) --> ∅
9291
end
9392
```
9493
Finally, we will again perform ensemble simulations of our model. This time, at the beginning of each simulation, we will use `make_K_series` to generate a new $K$, and set this as the `K_in` parameter's value.
9594
```@example noise_modelling_approaches
9695
function prob_func_Kin(prob, i, repeat)
97-
p = [ps; :K_in => make_K_series()]
98-
# return remake(prob; p)
99-
return ODEProblem(repressilator, prob.u0, prob.tspan, p)
96+
p = [ps; :K_in => make_K_series()]    
97+
    return ODEProblem(repressilator_Kin, prob.u0, prob.tspan, p)
10098
end
101-
oprob = ODEProblem(repressilator, u0, tend, [ps; :K_in => make_K_series()])
99+
oprob = ODEProblem(repressilator_Kin, u0, tend, [ps; :K_in => make_K_series()])
102100
eprob_inputnoise = EnsembleProblem(oprob; prob_func = prob_func_Kin)
103101
sol_inputnoise = solve(eprob_inputnoise; trajectories = 4)
104102
plot(sol_inputnoise; idxs = :X)
@@ -107,21 +105,26 @@ Like in the previous two cases, this generates heterogeneous trajectories across
107105

108106

109107
## [Investigating the mean of noisy oscillations](@id noise_modelling_approaches_model_noisy_oscillation_mean)
110-
Finally, we will consider what happens to the mean activity of $X$ as the heterogeneous trajectories diverge from their initial condition. First we create ensemble simulations with a larger number of trajectories.
108+
Finally, we will observe an interesting phenomenon for ensembles of stochastic oscillators. First, we create ensemble simulations with a larger number of trajectories.
111109
```@example noise_modelling_approaches
112-
sol_intrinsic = solve(eprob_intrinsic, ImplicitEM(); trajectories = 100)
113-
sol_extrinsic = solve(eprob_extrinsic; trajectories = 100)
110+
sol_intrinsic = solve(eprob_intrinsic, ImplicitEM(); trajectories = 200)
111+
sol_extrinsic = solve(eprob_extrinsic; trajectories = 200)
114112
```
115-
Next we can use the `EnsembleSummary` interface to plot each ensemble's mean activity (as well as 5% and 95% quantiles) over time:
116-
113+
Next, we can use the `EnsembleSummary` interface to plot each ensemble's mean activity (as well as 5% and 95% quantiles) over time:
117114
```@example noise_modelling_approaches
118115
e_sumary_intrinsic = EnsembleAnalysis.EnsembleSummary(sols, 0.0:1.0:tend)
119116
e_sumary_extrinsic = EnsembleAnalysis.EnsembleSummary(sols, 0.0:1.0:tend)
120-
plot(e_sumary_intrinsic; label = "Intrinsic noise")
121-
plot(e_sumary_extrinsic; label = "Extrinsic noise")
117+
plot(e_sumary_intrinsic; label = "Intrinsic noise", idxs = 1)
118+
plot!(e_sumary_extrinsic; label = "Extrinsic noise", idxs = 1)
122119
```
123-
Here we can see that while the individual simulations (as seen in the previous sections) maintain their amplitudes, the mean amplitude across the ensemble tends towards zero. This is a well-known phenomenon in circadian biology (and a reason for considering single-cell observations, i.e. by only observing the mean it is impossible to tell whether the individual trajectories experience de-synchronisation or dampening).
120+
Here we can see that, over time, the systems' mean $X$ activity reaches a constant level around $30$.
121+
122+
This is a well-known phenomenon (especially in circadian biology[^2]). Here, as stochastic oscillators evolve from a common initial condition the mean behaves as a damped oscillator. This can be caused by two different phenomena:
123+
- The individual trajectories are themselves damped.
124+
- The individual trajectories's phases get de-synchronised.
125+
However, if we only observe the mean behaviour (and not the single trajectories), we cannot know which of these cases we are encountering. Here, by checking the single-trajectory plots from the previous sections, we note that this is due to trajectory de-synchronisation. Stochastic oscillators have often been cited as a reason for the importance to study cellular systems at the *single-cell* level, and not just in bulk.
124126

125127
---
126128
## [References](@id noise_modelling_approaches_references)
127-
[^1]: [Michael B. Elowitz, Arnold J. Levine, Eric D. Siggia, Peter S. Swain, *Stochastic Gene Expression in a Single Cell*, Science (2002).](https://www.science.org/doi/10.1126/science.1070919)
129+
[^1]: [Michael B. Elowitz, Arnold J. Levine, Eric D. Siggia, Peter S. Swain, *Stochastic Gene Expression in a Single Cell*, Science (2002).](https://www.science.org/doi/10.1126/science.1070919)
130+
[^2]: [Qiong Yang, Bernardo F. Pando, Guogang Dong, Susan S. Golden, Alexander van Oudenaarden, *Circadian Gating of the Cell Cycle Revealed in Single Cyanobacterial Cells*, Science (2010).](https://www.science.org/doi/10.1126/science.1181759)

0 commit comments

Comments
 (0)