R/biotic_resource_competition.R

In gauravsk/ecoevoapps: Simulate Dynamics of Ecology/Evolution Models

#' Internal function for running the biotic resource
#' competition between two consumers and a shared resource species
#' @param time vector of time units over which to run model
#' @param init initial population size of the two consumers P1 and P2,
#' and the resource species H
#' @param params vector of parameters:
#' r = per capita growth rate of the resource species
#' q = 1/prey carrying capacity of the resource species H
#' T_h1, T_h2 = handling time for consumer species P1 and P2
#' a1, a2 = attack rate of consumer species P1 and P2
#' e1, e2 = conversion efficiency of consumer species P1 and P2
#' d1, d2 = predator death rate of consumer species P1 and P2
#' @keywords internal
biotic_competition <- function(time, init, params) {
  with (as.list(c(time, init, params)), {
    # description of parameters:
    dH_dt = r*H*(1 - q*H) - (a1*H*P1)/(1 + a1*T_h1*H) - (a2*H*P2)/(1 + a2*T_h2*H)
    dP1_dt = e1*(a1*H*P1)/(1 + a1*T_h1*H) - d1*P1
    dP2_dt = e2*(a2*H*P2)/(1 + a2*T_h2*H) - d2*P2
    return(list(c(dH = dH_dt, dP1 = dP1_dt, dP2 = dP2_dt)))

  })
}

#' Run a model of competition for a biotic resource
#'
#' @param params vector of model parameters:
#' c(r, q, a1, T_h1, e1, d1, a2, T_h2, e2, d2)
#' @param time vector of time units over which to run model, starting from 0.
#'   `time` can also be supplied as just the total length of the simulation
#'   (i.e. tmax)
#' @param init vector of initial population sizes for both consumer species,
#' with names P1 and P2, and resource species H
#' @import deSolve
#' @examples
#' run_biotic_comp_model(
#' time = seq(0,10),
#' init = c(H = 30, P1 = 25, P2 = 25),
#' params = c(r = 0.2, q = .0066,
#' a1 = .02, T_h1 = 0.1, e1 = 0.4, d1 = 0.1,
#' a2 = .02, T_h2 = 0.1, e2 = 0.39, d2 = 0.1))
#' @seealso [plot_biotic_comp_time()] for plots of the population dynamics over
#'   time

#' @export
run_biotic_comp_model <- function(time = seq(0,10),
                                  init = c(P1 = 25, P2 = 25, H = 30),
                                  params =  c(r = 0.2, q = .0066,
                                              a1 = .02, T_h1 = 0.1,
                                              e1 = 0.4, d1 = 0.1,
                                              a2 = .02, T_h2 = 0.1,
                                              e2 = 0.39, d2 = 0.1)) {

  if(length(time) == 1) {
    tmax <- time
    time <- seq(0, tmax, 0.1)
  } else if(time[1] != 0) {
    stop("The time vector should start at 0.")
  }

  # Check that init has been defined properly
  if(!(is.numeric(init))) stop("init should be a numeric vector of length 3, e.g. c(P1 = 25, P2 = 25, H = 30)")
  if(length(init) != 3) stop("init should be a numeric vector of length 3, e.g. c(P1 = 25, P2 = 25, H = 30)")
  if(!(all(names(init) %in% c("P1","P2", "H")))) stop("init should be a numeric vector of length 4, e.g. c(P1 = 25, P2 = 25, H = 30)")
  # Check that params is correctly defined
  if(!(is.numeric(params))) stop("params should be a numeric vector")
  if(!(all(c("r", "q", "a1", "T_h1", "e1", "d1", "a2", "T_h2", "e2", "d2") %in%
           names(params)))) {
    stop("Please provide a complete parameter vector (see ?run_biotic_comp_model for details)")
  }


  ode(func = biotic_competition, y = init, times = time, parms = params)
}



#'Plot population size over time for the boitic resource competition model
#'
#' @param sim_df matrix or data frame of simulations generated by
#'   [ecoevoapps::run_biotic_comp_model]
#' @examples
#' sim_df_biotic_rc <- run_biotic_comp_model(
#' time = seq(0,10),
#' init = c(P1 = 25, P2 = 25, H = 30),
#' params = c(r = 0.2, q = .0066,
#' a1 = .02, T_h1 = 0.1, e1 = 0.4, d1 = 0.1,
#' a2 = .02, T_h2 = 0.1, e2 = 0.39, d2 = 0.1))
#' plot_biotic_comp_time(sim_df_biotic_rc)
#' @import ggplot2
#' @import tidyr
#' @import dplyr
#' @seealso [run_biotic_comp_model()]  for simulating the biotic resource
#' competition model
#' @export
plot_biotic_comp_time <- function(sim_df) {

  # To suppress CMD Check
  H <- P1 <- P2 <- time <- value <- Population <- NULL

  sim_df <- data.frame(sim_df)

  # Reshape the data for plotting
  sim_df_long <- pivot_longer(sim_df, cols = c(H, P1, P2), names_to = "Population")

  # Make plot
  ggplot(sim_df_long) +
    geom_line(aes(x = time, y = value, color = Population), size = 2) +
    scale_fill_manual(values = c("#77AADD", "#FFAABB", "#AAAA00")) +
                        ylab("Population size") +
                        scale_x_continuous(expand = c(0, 0)) +
                        scale_y_continuous(expand = c(0, 0)) +
                        theme_apps()
}

#' Plot population size over time for the biotic resource competition model
#'
#' @param params vector of parameters describing the consumer species'
#' attack, consumption, and mortality rates, as well as their handling times;
#' and the resource species' intrinsic growth and density-dependent mortality rates
#' @examples
#' params <- c(r = 0.2, q = .0066,
#' a1 = .02, T_h1 = 0.1, e1 = 0.4, d1 = 0.1,
#' a2 = .02, T_h2 = 0.1, e2 = 0.39, d2 = 0.1)
#' plot_functional_responses(params)
#' @import ggplot2
#' @import tidyr
#' @import dplyr
#' @seealso [run_biotic_comp_model()]  for simulating the biotic resource
#' competition model
#' @export
plot_functional_responses <- function(params) {
  # To suppress CMD Check
  y1 <- y2 <- prey <- NULL

  fr_df <- data.frame(prey = 1:2000)
  fr_df <- fr_df %>%
    mutate(y1 = (params["a1"] * prey)/(1 + params["a1"] * params["T_h1"] * prey),
           y2 = (params["a2"] * prey)/(1 + params["a2"] * params["T_h2"] * prey))

  plot_fxnrep_sp1 <-
    ggplot(fr_df, aes(x = prey, y = y1)) +
    geom_path(col = "#77AADD", size = 2) +
    xlab("Prey Density") +
    ylab("Prey Consumed \nper Predator") +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0)) +
    theme_apps()

  plot_fxnrep_sp2 <-
    ggplot(fr_df, aes(x = prey, y = y2)) +
    geom_path(col = "#AAAA00", size = 2) +
    xlab("Prey Density") +
    ylab("Prey Consumed \nper Predator") +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0)) +
    theme_apps()

  to_return <- list(plot_fxnrep_sp1, plot_fxnrep_sp2)

  return(to_return)

}