R/CR_model.R

In Chang-Yu-Chang/MigrationCommunity: Microbial Consumer-Resource Model

#' Simulate consumer-resource dynamics
#' @description
#' This function takes a named vector as initial community, simulates based
#' on CR_ode function, and reports the output result matrix.
#' @name CR_model
#' @param community the initial state of species and resource in a community. A named vector.
#' @param time_limit time limit for ode.
#' @param time_step time step for ode.
#' @return a deSolve matrix
#' @examples
#' # Load example parameters
#' data(example_parameter)
#'
#' # Build regional species pool
#' pool_build()
#'
#' # Create a community, a named vector
#' set.seed(1)
#' species <- sample(1:S_pool, S, replace = FALSE) %>% sort()
#' community <-
#'   c(
#'     setNames(c(rep(1, P)), paste0("R", sprintf("%03d", 1:P))),
#'     setNames(c(rep(.1, length(species))), paste0("X", sprintf("%05d", species)))
#'   )
#'
#' # Run ode
#' result <- CR_model(community = community, time_limit = 100, time_step = 1)
#'
#' # Plot
#' result %>%
#'   as.tibble() %>%
#'   gather("index", "value", -1) %>%
#'   ggplot(aes(x = time, y = value, color = index, lty = substr(index, 1, 1))) +
#'   geom_line() +
#'   theme(legend.position = "none")
#'
#' # Main function for running consumer-resource model
CR_model <- function(
                     community = community, # Named vector with biomass as elements and consumer or resource in names.
                     time_limit = 100, # Time limit
                     time_step = 1) {

  # Consumer
  sp <- grep("X", names(community), value = T) %>%
    substr(2, 6) %>%
    as.numeric()

  # Resource
  re <- grep("R", names(community), value = T) %>%
    substr(2, 4) %>%
    as.numeric()

  # Time sequence
  time <- seq(0, time_limit, by = time_step)

  # Parameters: a named vector
  parameters <- list(sp = sp, re = re, J = J, Rin = Rin, P = P, u = u, w = w, m = m, stoi = stoi)

  # Initial condition: a named vector
  state <- community

  # Run ode
  result <- ode(y = state, times = time, func = CR_ode, parms = parameters)

  # Melt result
  result <-
    result %>%
    as.data.frame() %>%
    melt(id.var = "time") %>%
    as.tibble()

  #
  return(result)
}