R/pool_build.R

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

#' Build the regional species pool
#' @description Regional species pool
#' @name pool_build
#' @param seed the seed for sampling. Default = 1.
#' @param save_data whether to save species characteristic in a RData file. Logical.
#' @param data_name if `save_data` == TRUE, specify the name of the RData file. String.
#' @return ...
#' @examples
#' # Load packages and example parameters
#' library(tidyverse)
#' data(example_parameter)
#'
#' # Build regional species pool
#' pool_build(save_data = FALSE)
#'
#' # Function that defines species abilities: m, u, w, and stoi
pool_build <- function(
                       seed = 1, # Setting seed
                       save_data = TRUE,
                       data_name = "species_pool") {
  # Seed ----
  set.seed(seed)

  # Death rate; m ----
  ## A numeric vector with length of S
  m <- rep(mi, S_pool)

  # Uptaking rate; u ----
  ## A matrix with P rows and S_pool columns

  if (function_group == 0) {
    # Randomly drawn from uniform distribution, where u sum is 1
    f <- rep(1, S_pool)
    u <- matrix(runif((P - 1) * S_pool), nrow = P - 1, ncol = S_pool) %>%
      rbind(0, 1) %>%
      apply(2, sort)
    u <- matrix(u[-1, ] - u[-(P + 1), ], P, S_pool)
  } else if (is.numeric(function_group)) {
    # Funtional groups that have preference in uptaking resources
    # See "Emergent property SI" for how to conduct phylogenetic sampling
    suppressMessages(library(MCMCpack)) # Load this package for function that draws from dirichlet distribution

    # Empty matrix of uptaking rates
    u <- matrix(NA, nrow = P, ncol = S_pool)

    # Groups for species
    f <- rep(1:function_group, each = floor(S_pool / function_group))[1:S_pool]

    # Group: concentration parameter for resource in a group
    a_family <- rep(NA, S_pool)
    for (i in 1:S_pool) {
      a_preferred <- rnorm(1, group_mu, group_sigma) # Preferred resource
      a_preferred <- min(a_preferred, 1) # Set the maximum to 1
      a_other <- runif(P - 1) # Other resources

      u[f[i], i] <- a_preferred
      u[-f[i], i] <- (1 - a_preferred) * a_other / sum(a_other)
      a_family[i] <- a_preferred
    }

    # Species: total resource capacity
    t_species <- rnorm(S_pool, 1, 0.01)

    # Species: proportion of fixed budget of uptake resources
    K_family <- 100 # High K: species are very similar
    for (i in 1:S_pool) u[, i] <- rdirichlet(1, K_family * u[, i])

    # Uptaking rate
    u <- t_species * u * 0.9
  }

  # Probability of being a generalist; d_general ----
  u <- u * matrix(runif(S_pool * P) <= d_general, nrow = P, ncol = S_pool)

  # Conversion efficiency; w ----
  w <- matrix(wi, P, S_pool)

  # Stoichiometic matrix; stoi ----
  ## One available resource, no cross_feeding
  if (P == 1) {
    stoi <- array(0, dim = c(P, P, S_pool))
    for (i in 1:S_pool) diag(stoi[, , i]) <- -1

    ## Multiple resources, cross-feeding, different metabolism
  } else if (cross_feeding == TRUE & same_stoi == FALSE) {
    # Cross-feeding
    stoi <- array(NA, dim = c(P, P, S_pool))
    for (k in 1:S_pool) {
      for (i in 1:P) {
        w_i <- w[, k] / w[i, k] # denoted as ki in the note
        mu <- runif(1)
        d <- runif(P)
        stoi[i, , k] <- (mu * d) / c(d %*% w_i)
      }
    }
    for (i in 1:S_pool) diag(stoi[, , i]) <- -1

    ## Multiple resources, cross-feeding, same metabolism
  } else if (cross_feeding == TRUE & same_stoi == TRUE) {

    # Cross-feeding
    stoi <- array(NA, dim = c(P, P, S_pool))
    for (i in 1:P) {
      w_i <- w[, 1] / w[i, 1] # denoted as ki in the note
      mu <- runif(1)
      d <- runif(P)
      stoi[i, , 1] <- (mu * d) / c(d %*% w_i)
    }
    diag(stoi[, , 1]) <- -1
    stoi[, , -1] <- stoi[, , 1]
  }

  # Species-abundance distribution; trade-off ----
  if (trade_off == TRUE) {
    abundance_species <-
      rlnorm(n = S_pool, meanlog = 1, sdlog = 2) %>%
      round(2) %>%
      sort()

    a_preferred <- rep(NA, S_pool)
    for (i in 1:S_pool) a_preferred[i] <- u[f[i], i]
    abundance_species[rank(a_preferred)] <- abundance_species
  } else if (trade_off == FALSE) {
    if (species_abundance == "log-normal") {
      abundance_species <- rlnorm(n = S_pool, meanlog = 1, sdlog = 2) %>% round(2)
    } else if (species_abundance == "uniform") {
      abundance_species <- rep(1, S_pool)
    }
  }
  abundance_species <- abundance_species / sum(abundance_species)

  # Save the parameter in a RData ----
  if (save_data == TRUE) {
    save(m, f, u, w, stoi, abundance_species,
      file = paste0("data/", data_name)
    )
  }

  # Return the output to global environment
  m <<- m
  f <<- f
  u <<- u
  w <<- w
  stoi <<- stoi
  abundance_species <<- abundance_species

  # Regional pool with abundance and species names
  regional_pool <<- setNames(abundance_species, paste0("X", sprintf("%05d", 1:S_pool)))
}