---
title: "Noisy LG model: incorporating noise to the classical LG ODE clock model"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{noisy-LG-model}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
out.width = "80%",
dpi = 300
)
```
This vignette documents a "noisy LG model" provided in this package which
incorporates Gaussian noise in the classical LG model that was initially
defined as an ODE model. We take the classical LG circadian clock ODE model
and provide a stochastic differential equation (SDE) form.
## A SDE LG model
Here I summarize key concepts behind the `noisy_LG` model. Full definition is
available in the `odin/noisy_leloup-goldbeter.R` Odin definition file.
1. Converting the ODE model to difference equations
$\vec{X}(t+\mbox{d}t)=\vec{X}(t)+\vec{X'}(t)\mbox{d}t$. This is the
`discrete_LG` model.
2. Converting all concentration (states and parameters, in nM for LG)
to molecule count. This enables more physical addition of noise as noise
is inherently about molecule counts instead of concentration (mass statistics).
3. Addition of noise terms st $\vec{X}(t+\mbox{d}t)=\vec{X}(t)+\vec{X'}(t)\mbox{d}t\ +\ b(X,t)\mbox{d}W$ (Wiener process) by considering randomness of the underlying
reaction processes (see below "Noise Scaling" section).
The final difference equation update will be Euler-Maruyama where $\vec{X'}$ is
taken from the initial ODE:
$$
\vec{X}(t+\Delta t)=\vec{X}(t)+\vec{X'}(t)\Delta t+\vec{\sigma}(X,t)\Delta W,\ \Delta W=X\sqrt{\Delta t},\ X\sim\mathcal{N}(0, 1)
$$
### When is the model valid?
At molecular level, reactions can be modeled as Poisson processes. Given a
reaction with rate $\alpha$ [molecules/$\Delta t$], within one time step we
turnover $\alpha$ molecules (mean). If $\alpha>>1$, we can think of the process
as a sum of many iid processes (basically central limit theorem).
In this case, Wiener process would be good approximation.
### Noising scaling
Intuitively, faster process -> larger expectation -> larger variance. This
comes natural from the Poisson model. With the following update equation of
RNA level $M$:
$$
M(t+\Delta t) = M(t) + (\mbox{transcription}(t)\ -\ \mbox{degradation}(t))\Delta t
$$
At each update step, transcription and degradation are independent
Poisson. We therefore have the noise-incorporated model:
$$
M(t+\Delta t) = M(t) + (\mbox{txn}(t)\ -\ \mbox{deg}(t))\Delta t\ +\ \sqrt{\mbox{txn}(t)+\mbox{deg}(t)}\ \Delta W
$$
**The noisy LG model allows you to "turn on/off" noise terms for updating of each of the 10 states in the LG model.**
## Loading the noisy LG model
Below we load the noisy LG model. Note that `getOdinGen()$noisy_LG` is a list:
- `$gen`: model generator
- `$count2nM`: multiplier to convert simulation result (count) to concentration (nM).
In this model, we assumed the following:
- Cell dimension `cellDimUM` takes default `10um` (cube of 10um sides).
- Parameters are converted from concentration (`_conc` postfix user variables) to count using volume.
- By default, only `M_T` (timeless RNA) noise is turned on (`NoiseVariance_M_T = -1`).
- Simulation is performed in count and hour units
If you would like to still analyze simulation data in concentration unit, multiply
results by the `$count2nM` multiplier (see example below).
To turn on/off noise for each state, change the `NoiseVariance` prefix user
variables. Setting to -1 will add Poisson-scaled noise while setting to a positive
value (in nM concentration) will add a constant variance (=given variance) noise.
```{r setup}
library(clockSim)
library(ggplot2)
library(dplyr)
mg <- getOdinGen()$noisy_LG
model <- mg$gen$new()
# Time steps
interval_hours <- 0.001
total_hours <- 2400
model$set_user(STEP_HOURS = interval_hours)
steps <- seq(from = 1, to = total_hours / interval_hours)
```
## Example: noise in tim RNA `M_T` turnover
`M_T` noise is turned on by default. Here, we compare results with this noise
off and on (refer to other vignettes in the package for details on these analyses):
1. Phase portrait plot
2. Periodogram
```{r}
# A helper function for running simulation and convert unit
runSim <- function(model, steps, interval_hours){
res <- model$run(steps)
res <- res |> as.data.frame() |> mutate(time = step*interval_hours)
res <- res |> filter(time %% 1 == 0) # Get results every hour
plt1 <- plot_phase(
res |>
select(step, time, C, C_N) |>
# Convert units from count to nM
mutate(across(-c(time, step), .fns = \(x) x*mg$count2nM)),
C, C_N)
res.per <- res |> tail(n=240) # Only use last 240 hours for periodogram
per <- compute_period(
# Use M_T tim RNA data
res.per |> pull(M_T),
# Refer to the lomb::lsp().
# Only consider 18-30hour period, ofac=2 gives around ~1hour precision
method = "lomb", from = 18, to = 30, type = "period", ofac = 2)
return(list(
plt.phase = plt1, lomb = per, res = res
))
}
```
### Phase portrait and periodogram
```{r}
# Turn noise off
model$set_user(NoiseVariance_M_T = 0)
run <- runSim(model, steps, interval_hours)
plot(run$plt.phase)
print(run$lomb)
# Turn noise on
model$set_user(NoiseVariance_M_T = -1)
run <- runSim(model, steps, interval_hours)
plot(run$plt.phase)
print(run$lomb)
```
How about 1.5X k_sT (tim translation rate)?
```{r}
model$set_user(k_sT = 1.8)
# Turn noise off
model$set_user(NoiseVariance_M_T = 0)
run <- runSim(model, steps, interval_hours)
plot(run$plt.phase)
print(run$lomb)
# Turn noise on
model$set_user(NoiseVariance_M_T = -1)
run <- runSim(model, steps, interval_hours)
plot(run$plt.phase)
print(run$lomb)
```