R/CR_ode.R

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

#' ODE of consumer-resource model
#' @description The ordinary differential equations that describe the
#' consumer and resource interactions.
#' @name CR_ode
#' @param time time sequence for which output is wanted; the first value of times must be the initial time.
#' @param state the initial (state) values for the ODE, a vector.
#' @param parameters parameters passed to ODE.
#' @return a matrix of class deSolve.
#' @examples
#' # Load packages and example parameters
#' library(tidyverse)
#' library(deSolve)
#' data(example_parameter)
#'
#' # Create species pool
#' pool_build()
#'
#' # Time sequence
#' time <- seq(0, 100, by = 1)
#'
#' # Parameters: a named vector
#' sp <- 1:S
#' re <- 1:P
#' parameters <- list(sp = sp, re = re, J = J, Rin = Rin, P = P, S = S, u = u, w = w, m = m, stoi = stoi)
#'
#' # Initial condition: a named vector
#' state <- c(
#'   setNames(rep(.1, P), paste0("R", sprintf("%03d", 1:P))),
#'   setNames(rep(1, length(sp)), paste0("X", sprintf("%05d", sp)))
#' )
#'
#' # Run ode
#' result <- ode(y = state, times = time, func = CR_ode, parms = parameters)
#'
#' # Plot
#' result %>%
#'   as.tibble() %>%
#'   gather("index", "value", -1) %>%
#'   ggplot(aes(x = time, y = value, color = index, lty = substr(index, 1, 1))) +
#'   geom_line()
CR_ode <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    # Create a vector of R
    RR <- sapply(paste0("R", sprintf("%03d", re)), function(x) eval(parse(text = x)))

    # Create a vector of X
    XX <- sapply(paste0("X", sprintf("%05d", sp)), function(x) eval(parse(text = x)))

    # Resource derivatives
    dR <- J * (Rin[re] - RR) # Supplied resources
    for (i in 1:length(sp)) dR <- dR + stoi[, , sp[i]] %*% (u[, sp[i]] * RR) * XX[i] # Secreted metabolites; loop through species

    # Consumer derivatives
    dX <- XX * (t(u[, sp] * w[, sp]) %*% RR - m[sp])

    # Return dR and dX
    return(list(c(dR, dX)))
  })
}