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

# !diagnostics off
knitr::opts_chunk$set(cache = T, echo = T)

Install package

Install this package from github using devtools.

if (FALSE) {
  install.packages("devtools")
  library(devtools)
  install_github("Chang-Yu-Chang/micro.crm")
}

library(micro.crm)

Model description

This model is conceptually based on the MacArthur's consumer-resource model but the relation between consumers and resources is modified to better describe the scenarios in the microbial interactions. In particular, the stoichiometric metabolism, uptaking rate, and resource-using efficiency of a consumer are included.

The model includes two sets of equations: consumers and resources. Given that total number of resource type $P$ and total number of consumers $S$ in this system, the growth rate of consumer $i$ is defined by

$$\frac{dX_i}{dt}=\sum^P_{j=1} R_jX_iw_j^iu_j^i-m_iX_i \quad i = 1,2,...,S$$

where $X_i$ is the biomass of consumer $i$, $R_j$ is the biomass of resource $j$, $m_i$ is the density-independent loss rate of consumers $i$, $w_j^i$ is the conversion efficiency of resource $j$ converted into biomass of consumer $i$, and $u_j^i$ denotes the uptaking rate for consumer $i$ to use resource $j$.

Resource is supplied in chemostat fashion. Let the dynamics of resource $j$ be given by

$$\frac{dR_j}{dt}=J(R_{supply,j}-R_j)+\sum_{i=1}^{S}\sum_{k=1}^{P} R_kX_iu_{k}^{i}w_k^{i} D_{kj}^i$$

The first term in this equation defines the resource supplement from environment to this system, while the second term denotes the supplement from consumers via cross-feeding. In the first term, resource is supplied at flow rate $J$ and inflow concentration $R_{supply,j}$. In the second term, $D_{jk}^i$ is a stoichiometric matrix that describes the metabolism of consumer $i$ which uptakes and transforms resource $k$ into secreted resource $j$.

Examples originally from README.md

The first step is to call the example parameters for the model. The parameters will be saved in global environment.

data(example_parameter)

To generate a species regional pool, use the function pool_build. The parameters will be save in a Rdata file, data/species_pool.RData.

pool_build(save_data = TRUE, data_names = "species_pool")

From this regional pool, subset micrbial communities by given species abundance. The output community compostion is a named vector.

com1 <- community_generate(I = .01, threshold = 1e-3, r = 1)

Now let's take this generated community com1 as intial community. Because there are only consumers in the generated community, another named vector as the initial resource condition has to be concatenated.

resource <- setNames(rep(1, P), paste0("R", sprintf("%03d", 1:P)))
result <- CR_model(c(resource, com1), time_limit = 100, time_step = 1)

result %>%
  as.data.frame() %>%
  gather("index", "value", -1) %>%
  mutate(type = substr(index, 1, 1)) %>%
  ggplot(aes(x = time, y = value, color = index, lty = type)) +
  geom_line() + 
  theme_bw() +
  theme(legend.position="none") +
  NULL

In the cases of migration, co-assembly, or community coalescence, the communities (named vectors) can be coalesced using function community_coalesce().

com2 <- community_generate(I = .01, threshold = 1e-3, r = 2)
community_coalesce(com1, com2)

Self-assembly

Run the following scripts in bash environment.

# Self-assembly 
Rscript self-assembly.R self-assembly001

# Isolation
Rscript self-assembly-isolate.R self-assembly001

# Pairwise competition
Rscript self-assembly-pair.R self-assembly001 0.5 0.5

Examples

This section will cover two simple cases: consumer-resource dynamics and migration. To simply the simulation, the resource is only limited to one single supplied reosurce while growth on cross-feeding is allowed.

Consumer-resource dynamics

The first step is to call the example parameters for the model. The parameters will be saved in global environment.

data(example_parameter)

pool_build generates a species regional pool that includes the characteristics for consumers. The parameters will be saved in a Rdata file, data/species_pool.RData.

pool_build(save_data = TRUE, data_name = "species_pool")

From this regional pool, subset microbial communities by given species abundance. The output community composition is a named vector.

com1 <- community_generate(I = .01, threshold = 1e-3)
com1

Now let's take this sampled community com1 as the initial community. Because there are only consumers in the generated community, another named vector as the initial resource condition has to be concatenated.

resource <- setNames(rep(1, P), paste0("R", sprintf("%03d", 1:P)))
resource

Run the consumer-resource dynamics simulation. CR_model takes a single vector of initial conditions so we have to concatenate the resource and com1. The output is a melted dataframe.

result <- CR_model(c(resource, com1), time_limit = 100, time_step = 1)
head(result)

Plot the dynamics.

result %>%
  mutate(type = substr(variable, 1, 1)) %>%
  ggplot(aes(x = time, y = value, color = variable)) +
  geom_line() +
  theme_bw() +
  facet_grid(type ~ .) +
  theme(legend.position = "none") +
  NULL

# Write result file
parameter_write(result = result, parameter_list = parameter_list, path = "data/", treatment_name = "example", r = 1)

Run replicates

Save the code chunch below to another R script example001.txt and run it in bash.

# Make regional species pool
parameter_load(parameter_text = "parameter_script/example001.txt")
pool_build(save_data = TRUE, data_name = "species_pool.Rdata")

# Make resource environment
resource <- setNames(c(1, rep(0, P - 1)), paste0("R", sprintf("%03d", 1:P)))
assign("Rin", c(Rin1, rep(0, times = P - 1)), envir = .GlobalEnv) # Supply for each resource. Set to 0

# Run
for (r in 1:20) {
  # Simulation
  result <- community_generate(pool = regional_pool, I = I, threshold = threshold, seed = r) %>%
    c(resource, .) %>%
    CR_model(time_limit = time_limit, time_step = 1)

  # Write result file to txt
  parameter_write(
    result = result,
    parameter_list = parameter_list,
    file_name = paste0("simulation_result/example001_", sprintf("%03d", r), ".txt")
  )

  # Report progress
  print(r)
}

CR model plot

Read the data for example001_001.txt to example001_020.txt.

rm(list = ls())
# Load parameters
parameter_load(parameter_text = "parameter_script/example001.txt")

result <- result_read(dir = "simulation_result/", pattern = "example001", threshold = 1e-5)

Resource dynamics

result %>%
  filter(type == "R") %>%
  ggplot(aes(x = time, y = value, color = variable, lty = type), lwd = .1) +
  geom_line() +
  facet_wrap(. ~ replicate) +
  theme_bw() +
  theme(legend.position = "none") +
  labs(x = "time") +
  ggtitle("Resource")

Consumer dynamics

result %>%
  filter(type == "X") %>%
  ggplot(aes(x = time, y = value, color = variable, lty = type), lwd = .1) +
  geom_line() +
  facet_wrap(. ~ replicate) +
  theme_bw() +
  theme(legend.position = "none") +
  labs(x = "time") +
  ggtitle("Consumer")

Final community composition at equilibrium.

result %>%
  filter(type == "X", time == max(time)) %>%
  arrange(value) %>%
  ggplot(aes(x = replicate, y = value, fill = function_group)) +
  geom_bar(stat = "identity", position = "fill", color = "black") +
  facet_grid(. ~ treatment) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  #  scale_color_discrete("black") +
  theme_bw() +
  guides(fill = guide_legend(title = "Family"), color = F) +
  labs(x = "replicate", y = "relative abundancee")