R/CR_ode.R

Defines functions CR_ode

Documented in CR_ode

#' 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 intial (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()
#' @export
#' @import tidyverse
#' @import deSolve
#' @import data.table



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


  
Chang-Yu-Chang/MigrationCommunity documentation built on April 22, 2019, 7:36 p.m.