output/script/self-assembly-batch.R

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

library(tidyverse)
library(deSolve)
library(data.table)
library(micro.crm)

#rm(list= ls())
# Augument passed from command line
args = commandArgs(trailingOnly=TRUE)
treatment_name = args[1]

# Make regional species pool
#treatment_name <- "self-assembly-batch001"
parameter_load(parameter_text = paste0("parameter_script/", treatment_name, ".txt"))
pool_build(save_data = F)

# Make resource environment
resource <- setNames(c(R1_initial, 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

# Modify the `w` and `stoi` of fermenters.
w[1,1:(S_pool/function_group)] <- w_fermenter_P1
stoi[2,1,1:(S_pool/function_group)] <- stoi_fermenter_P1_P2

# Run
for (r in 1:replicate) {
  result_list <- rep(list(NA), migration_event) 

  # First transfer
  comm <- NULL
  result <- community_generate(pool = regional_pool, I = I, threshold = threshold, seed = 1*r) %>% # Draw from pool
    community_coalesce(resource, comm / sum(comm) * I)  %>% # Consumers and resource from previous transfer
    CR_model(time_limit = migration_time_limit, time_step = 1)

  # Save result to the main list
  result_list[[1]] <- result
  
  # Take the final time step
  temp <- filter(result, time == max(time), value >= threshold) %>% dplyr::select(-time)
  comm <- setNames(temp$value, temp$variable)
  
  cat(paste0(" ", 1))
  
  # Later transfer
  for (i in 2:migration_event) {
    # Simulation
    result <- community_coalesce(resource, comm / sum(comm) * I) %>% # Consumers and resource from previous transfer
      CR_model(time_limit = migration_time_limit, time_step = 1)
    
    # Time
    result$time <- result$time + migration_time_limit*(i-1)
    
    # Save result to the main list
    result_list[[i]] <- result
    
    # Take the final time step
    temp <- filter(result, time == max(time), value >= threshold) %>% dplyr::select(-time)
    comm <- setNames(temp$value, temp$variable)
    
    cat(paste0(" ", i))
  }

  # Merge all migration events
  result_merged <- rbindlist(result_list)
  
  # Write result file to txt
  parameter_write(result = result_merged, 
                  parameter_list = parameter_list, 
                  file_name = paste0("simulation_result/", treatment_name,"_", sprintf("%03d", r),".txt"))
  
  # Report progress
  print(r)
}






if (FALSE) {
result_merged <- rbindlist(result_list[1:3])

mutate(result_merged, type = substr(variable, 1, 1)) %>%
ggplot(aes(x = time, y = value, color = variable)) +
  geom_line() +
  facet_grid(type~., scale = "free_y") +
  guides(color = F) +
  theme_bw()

}