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

# !diagnostics off
knitr::opts_chunk$set(cache = T, echo = F)
library(micro.crm)

This script includes the following steps:

1. Self-assembly. The simple scenario only has two functional groups: fermenters and non-fermenters in the species pool.

2. Isolation. This step reads the result from self-assebled communities and subset the communities at final time step.

3. Pairwise competition within communities. For each self-assembled communities, conduct pairwise competition.

Self-assembly

Example

Parameters for the two functional groups.Set up species pool and initial resource environment.

rm(list=ls())

# Make regional species pool
parameter_load(parameter_text = "self-assembly002.txt")
#save(list = ls(), file = "data/example_parameter.RData")
pool_build(save_data = TRUE, data_name = "species_pool.Rdata")

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

Modify the w and stoi of fermenters.

w[1,1:(S_pool/function_group)] <- .5
stoi[2,1,1:(S_pool/function_group)] <- 0.9

Run one simulation

result <-  community_generate(pool = regional_pool, I = I, threshold = threshold) %>%
  c(resource, .) %>%
  CR_model(time_limit = time_limit, time_step = 1)

Plot the dynamics

result %>%
  mutate(type = substr(variable, 1, 1)) %>%
  ggplot(aes(x = time, y = value, color = variable)) +
  geom_line() +
#  scale_x_continuous(limits = c(0, 50)) +
  theme_bw() +
  facet_grid(type~., scale = "free_y") +
  theme(legend.position="none") +
  NULL

Run self-assembly

Run the following chunk in bash. The example code above is in a R script self-assembly.R

```{bash eval = F} Rscript self-assembly.R self-assembly001

Read result from `self-assembly001.R`

```r
rm(list=ls())
parameter_load(parameter_text = "self-assembly001.txt")
result = result_read(dir = "simulation_result/", pattern = "self-assembly001", threshold = threshold)

Plot

Plot dynamics

result %>%
  filter(type == "R") %>%
  ggplot(aes(x = time, y = value, color = variable)) +
  geom_line() +
  scale_x_continuous(limits = c(0, 50)) +
  theme_bw() +
  facet_wrap(.~replicate) +
  theme(legend.position="none") +
  ggtitle("Resource") 

result %>%
  filter(type == "X") %>%
  ggplot(aes(x = time, y = value, color = variable)) +
  geom_line() +
#  scale_x_continuous(limits = c(0, 50)) +
  theme_bw() +
  facet_wrap(.~replicate) +
  theme(legend.position="none") +
  ggtitle("Consumer")

Plot final composition

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

Isolation

Run the following chunk in bashSave the chunk below to a R script self-assembly002-isolate.R.

```{bash eval = F} Rscript self-assembly-isolate.R self-assembly001

Read result

```r
rm(list=ls())
parameter_load(parameter_text = "self-assembly001.txt")
result_isolate = result_isolate_read(dir = "simulation_result/isolates/", pattern = "self-assembly001", threshold = threshold)

result_isolate %>% filter(time == max(time), type == "X")

Pairwise competition

```{bash eval = F} Rscript self-assembly-pair.R self-assembly001 0.5 0.5

Read result
```r
rm(list=ls())
parameter_load(parameter_text = "parameter_script/self-assembly001.txt")
result_pair <- result_pair_read(dir = "simulation_result/pairs/", pattern = "self-assembly001", threshold = threshold)

Plot dynamics in one community/replicate

result_pair %>%
  filter(type == "X", replicate == 1) %>%
  ggplot(aes(x = time, y = value, group = variable, color = function_group)) +
  geom_line() +
  scale_x_continuous(limits = c(0, 50)) +
  theme_bw() +
  facet_wrap(.~pair) +
  guides(color = guide_legend(title = "Family")) +
  ggtitle("Consumer") 

apply(u[,c(122, 140)], 2, function(x)order(x,decreasing = T))
apply(u[,c(122, 214)], 2, function(x)order(x,decreasing = T))

u[,c(122, 140)] %>% 
  as.data.frame() %>% setNames(c("X122", "X140")) %>%
  mutate(resource = 1:10) %>%
  gather("consumer", "uptaking", 1:2) %>%
  ggplot(aes(x = resource, color = consumer, y = uptaking)) +
  geom_point() +
  scale_x_continuous(breaks = 1:10) +
  scale_y_continuous(limits = c(0,1)) +
  theme_bw()

All pairwise interactions

pairs <- result_pair %>%
  filter(time == max(time), type == "X") %>%
  dplyr::select(-time) %>%
  group_by(treatment, replicate, mixing_ratio, pair) %>%
  summarize(interaction = ifelse(any(value==0), "exclusion", "coexistence"))

# Isolates
pairs$isolate1 <- sapply(strsplit(pairs$pair, "_"), function(x)`[`(x,1))
pairs$isolate2 <- sapply(strsplit(pairs$pair, "_"), function(x)`[`(x,2))

# Function groups of isolates
pairs <- pairs %>%
  mutate(family1 = cut(as.numeric(substr(isolate1, 2, 6)), breaks = seq(0, S_pool, S_pool/function_group), labels = 1:function_group),
         family2 = cut(as.numeric(substr(isolate2, 2, 6)), breaks = seq(0, S_pool, S_pool/function_group), labels = 1:function_group))

# Ratio of coexistence
pairs$two_families <- ifelse(pairs$family1 != pairs$family2, "F-N",
                          ifelse(pairs$family1 == 1, "F-F", "N-N"))

pairs %>% 
  group_by(treatment, mixing_ratio, interaction, two_families) %>%
  ggplot(aes(x = two_families, group = interaction, fill = interaction)) +
  geom_bar(position = position_dodge(preserve = "single"), col = "black") +
  theme_bw() + 
  theme(panel.grid = element_blank(),
        legend.position = "top",
        legend.title = element_blank()) +
  labs(x = "") +
  ggtitle("")

Chang-Yu-Chang/MigrationCommunity documentation built on Aug. 13, 2019, 9:41 p.m.