inst/doc/noisy-LG-model.R

In clockSim: Streamlined Simulation of Circadian Gene Networks

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  out.width = "80%",
  dpi = 300
)

## ----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)

## -----------------------------------------------------------------------------
# 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
  ))
}

## -----------------------------------------------------------------------------
# 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)

## -----------------------------------------------------------------------------
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)