R/model.R

In meyer-lab/cell-cycle-growth: Analysis of Cell Cycle Reporter Data

#' Run model to time t
#'
#' @param x The unknown parameters
#' @param start The starting point of the model.
#' @param t The time to simulate.
#'
#' @return The results of the simulation at \code{t}.
#' @export
#'
#' @examples
#' cycleModelSim(c(0.1, 0.2, 0.0, 0.01), time = 10.0)
#'
cycleModelSim <- function(x, start = c(1.0, 1.0, 0.0), time) {
  A <- t(matrix(c(-x[2] - x[3],         2*x[1], 0,
                  x[2],   -x[1] - x[4], 0,
                  x[3],           x[4], 0), 3, 3))
  start <- as.vector(start)

  outt <- rexpokit::expokit_wrapalldgexpv_tvals(A, time, force_list_if_1_tval = T)

  outt <- lapply(outt, function(x) x %*% start)

  return(t(do.call(cbind, outt)))
}

#' Model wherein cells follow the control condition up to a delay
#'
#' @param xBefore The parameters before the delay.
#' @param xAfter The parameters after the delay.
#' @param delay The delay time.
#' @param start The starting point.
#' @param t A vector of time values for which to solve.
#'
#' @return Matrix with the model simulation.
#' @export
#'
#' @examples
cycleModelDelay <- function(xBefore, xAfter, delay, start = c(1.0, 1.0, 0.0), t) {
  midpoint <- cycleModelSim(xBefore, start = start, time = delay)

  tBefore <- subset(t, t <= delay)
  tAfter <- subset(t, t > delay)

  if (length(tBefore) > 0) {
    modelBefore <- cycleModelSim(xBefore, start = start, time = tBefore)
  } else {
    modelBefore <- double()
  }

  if (length(tAfter) > 0) {
    modelAfter <- cycleModelSim(xAfter, start = midpoint, time = tAfter - delay)
  } else {
    modelAfter <- double()
  }

  return(rbind(modelBefore, modelAfter))
}



#' Calculates the steady-state growth rate and ratio of cell cycle phases.
#'
#' @param x The unknown parameters of the model.
#'
#' @return
#' @export
#'
#' @examples
getSS <- function(x) {
  b <- x[4] + x[1] - x[3] - x[2]

  R <- (b + sqrt(b**2 + 8 * x[1] * x[2]))/2/x[2]

  growthRate <- (x[1] - x[4] - x[3] * R) / (1 + R)

  return(c(growthRate, R))
}


#' Title
#'
#' @param growthGoal
#' @param Rgoal
#'
#' @return
#' @export
#'
#' @examples
solveFromOutcome <- function(growthGoal, Rgoal) {
  beta <- 0.03

  model <- function(x) {
    b <- x[2] - beta
    F1 <- (b + sqrt(b**2 + 8 * x[2] * beta))/2/beta - Rgoal
    F2 <- (x[2] - x[1] - x[1] * Rgoal) / (1 + Rgoal) - growthGoal
    c(F1, F2)
  }

  # first solution
  ss <- rootSolve::multiroot(f = model, start = c(0.0001, 0.5), positive = T)

  return(ss$root)
}