Scanning parameter grid to understand circadian clock model"

In clockSim: Streamlined Simulation of Circadian Gene Networks

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

This vignette documents the grid_scan() function provided in this package which allows high-performance repeated simulation runs over a grid of parameters and/or initial states. Here, we take the classical LG circadian clock model and demonstrate the following:

1. Under a wide range of initial states, LG model converges to a determined limit cycle attractor reflecting the circadian clock.
2. Such stability holds even when model parameters changes (albeit shape/period of the attractor will change).

Workflow shown here shall be generally useful for interpretation of wet-lab experiment results as study of clock mutants typically ends up with hypotheses where certain parameter(s) of the circadian clock is/are affected.

grid_scan() is generally useful for cases where you have an Odin model too; this function is decoupled from specifics of circadian clock models.

library(clockSim)
library(matrixStats)
library(dplyr)

Baseline: a limit cycle oscillator

Under wildtype condition the circadian clock is a limit cycle oscillator with \~24hr period. The LG model reproduces such behavior. To see this we run simulation with default setting:

model_gen <- getOdinGen()$continuous_LG
model <- model_gen$new()

sim_hours <- seq(from = 0, to = 2400, by = 1)
res <- model$run(sim_hours) |> as.data.frame()
res$time <- res$t
plot(plot_phase(res, M_T, C_N))
plot(plot_timeSeries(res, 0, 240, 1, 6, M_T, C_N))
print(compute_period(res$M_T |> tail(n = 240), method = "lomb"))

As shown above:

1. A stable limit cycle attractor is observed.
2. First 10 days show stabilization against the all-0 initial state.
3. Last 10 days show a 23.9hr period reflecting a circadian clock.

While promising, however, the LG model is a complex 10-state system. We next explore different properties of the system:

1. Would different initial states always lead to the clock limit cycle?
2. Would increase/decrease of parameters disrupt the attractor?

Initial state grid scan

10-variable makes it impossible to do a extensive grid on initial states. To see how fast our simulation runs:

run_eta(model, sim_hours)

With each simulation at \~10ms, we should be able to do 1) a very coarse grid on all 10 states with \~3 for each state, or 2) an extensive grid on selected states (say TIM/PER RNA):

$$ 3^{10}10\mbox{ms}\sim 250^{2}10\mbox{ms}\sim600\mbox{sec} $$

With the CRAN build time limitation, in this vignette we will not run the \~10min grid above but instead deal with a smaller grid on TIM/PER RNA.

Grid generation

First, we define the grid. We compute summary of all states in the baseline simulation and use the min/mean/max to define the grid. For each state, we scan N multiples of mean+N*spread of the baseline simulation, where spread = (max-min)/2 and N = 2,3, ..., Nmax. For the CRAN grid Nmax=5.

# Compute summary
summary <- res |>
  select(-t, -time) |>
  apply(2, summary)
# Only keep min/mean/max
summary <- summary[c(1,4,6),]
# Add on mean+N*spread, N=2,3,...,N_max
N_max <- 3 # For larger scan increase this. CRAN=3
get_multiples <- function(s, k) {
  # Extract components
  min <- s[1, ]
  mean <- s[2, ]
  delta <- s[3, ] - min

  # Create new rows using vectorized operations
  multiples <- outer(k, delta) |> sweep(2, mean, "+")
  attr(multiples, "original") <- s

  # Return
  multiples
}
summary <- get_multiples(summary, 2:N_max)
summary <- summary[,c("M_T", "M_P")] # Only RNA states
# Create grid
grid <- expand.grid(
  summary |> as.data.frame(), KEEP.OUT.ATTRS = FALSE)
#   User variables for initial state start with setUserInitial_
names(grid) <- paste0("setUserInitial_",names(grid))

Running grid_scan()

Now perform grid scan of initial states. grid_scan() is designed for efficient repeated model execution on parameter grid. Refer to its help for more details, and here I summarize its design concept:

1. Allow scanning arbitrary grid generated by, e.g., expand.grid().
2. Parallel computing based on base R parallel.
3. Return raw run data or allow user to apply a summary function (to reduce memory load).

We want the following statistics:

1. Did the simulation converge (i.e., successful integration)?
2. Did the simulation fall into the "default" attractor (i.e., the solution when initial states are all-zero default)? For this, we compute max NRMSE (among states) and min cosine similarity (among simulation time) to capture the largest deviation. For explanation of these metrics, refer to ?compute_rmse and ?compute_cosine.

These statistics mean we apply the following function to each run data:

default_attractor <- model_gen$new()$run(sim_hours)
default_attractor <- 
  default_attractor[(length(sim_hours)-240):length(sim_hours),]
stat.fn <- function(raw_run, reference = default_attractor){
  # Return code == 2 means successful integration (at least for lsoda)
  succ <- attr(raw_run, "istate")[1] == 2
  # Subset only the last 240 time points - should be stabilized
  raw_run <- raw_run[(nrow(raw_run)-240):nrow(raw_run),]
  # Compute normalized RMSE
  nrmse <- compute_rmse(raw_run, reference, normalize = "range")
  nrmse <- max(nrmse)
  # Compute cosine similarity
  cos <- compute_cosine(raw_run, reference)
  cos <- min(cos)
  # Return
  c(converged = succ, nrmse = nrmse, cos = cos)
}

Before running, let's see how this statistics function will slow down simulation. As shown, we add about \~10% time and \~40% memory footprint, which is minimal.

print(bench::mark(stat.fn(model$run(sim_hours))))
print(bench::mark(model$run(sim_hours)))

Now let's run the grid scan of initial states. We run grid_scan, and then wrangle the results into a data frame. As below, all converged to the same attractor.

scan <- 
  grid_scan(model_gen, grid, apply.fn = stat.fn, 
            n.core = 1, custom.export = "default_attractor",
            sim_hours)

process_scan <- function(){
  .scanDF <- scan |> unlist(use.names = FALSE) |> matrix(ncol = 3, byrow=TRUE)
  colnames(.scanDF) <- names(scan[[1]])
  result <- cbind(grid, .scanDF |> as.data.frame())
  summary(
    result |> select(converged, nrmse, cos)
    )
}

process_scan()

Plot two phase diagrams as an example (note that I rerun grid_scan() with default apply_fn = identity to get the raw data for plotting):

# Rerun scan
scan <- 
  grid_scan(model_gen, grid, apply.fn = identity,
            n.core = 1, custom.export = "default_attractor",
            sim_hours)
# Show first and last of grid
first <- grid[1,]
last <- grid[nrow(grid),]
print(first)
print(last)
plot(plot_phase(scan[[1]] |> as.data.frame(), M_T, C_N))
plot(plot_phase(scan[[nrow(grid)]] |> as.data.frame(), M_T, C_N))

How about modifying one parameter and see if that affects initial state stability with the same grid? Just add that parameter new value as a constant to the grid!

Remember that our goal is to see whether parameter change causes convergence property change - so our reference attractor should be computed with the new parameter too. As below, just as an example, doubling translation rate of TIM mRNA does not affect the behavior - all initial states still converged to the clock attractor.

new_grid <- grid
new_grid$k_sT <- 1.8 # 2X of original, check model$content()
new_model <- model_gen$new()
new_model$set_user(k_sT = 1.8)
default_attractor <- new_model$run(sim_hours)
default_attractor <- 
  default_attractor[(length(sim_hours)-240):length(sim_hours),]
# Then, repeat the same code above.
scan <- 
  grid_scan(model_gen, new_grid, apply.fn = stat.fn, 
            n.core = 1, custom.export = "default_attractor",
            sim_hours)

process_scan()

clockSim documentation built on April 11, 2025, 5:40 p.m.