# !diagnostics off knitr::opts_chunk$set(cache = T, echo = T)
Install this package from github using devtools
.
if (FALSE) { install.packages("devtools") library(devtools) install_github("Chang-Yu-Chang/micro.crm") }
library(micro.crm)
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$.
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)
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
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.
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)
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) }
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")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.