In kdorheim/MEMC: Microbial Explicit Model Comparison: a R package supporting model exploration

In this vignette

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>", 
  fig.width = 6, fig.height = 5
  )

The 2013 MEND model, originally documented in Wang et. al 2013 is a system of differential equations that describe soil carbon dynamics. There are 8 different carbon pools represented as the circles that are connected by 12 arrows, which represent the fluxes between carbon pools. For more information about the MEND 2013 model we encourage users to read Wang et. al 2013.

\

{width=70%}

\

Running default MEND

library(MEMC)     # MEMC should already be installed, see installation instructions. 
library(ggplot2)  # package used to visualize results 
library(knitr)    # makes nice tables
library(magrittr) # import the %>% pipeline

# set a theme to use in all the plots
theme_set(theme_bw()) 

The MEMC package has pre-built model configurations including one titled MEND_model which refers to the model published by Wang et. al 2013, run help(MEND_model) for more details and examples on how to solve.

All of the model configurations included in the MEMC package can be used directly with solve_model. Before trying to solve the model let's take a look at the model configuration. It is a named list that containing the following

The name of the model configuration.

MEND_model$name

A vector of the initial state values.

MEND_model$state

An environment where the parameters and inital state values are defined, this is made with configure_model.

class(MEND_model$env)

A function defining the carbon pools, the object returned from carbon_pools.

head(MEND_model$carbon_pools_func)

A function defining the carbon fluxes, this is returned by carbon_fluxes and cacn be modified by modify_fluxes_func.

head(MEND_model$carbon_fluxes_func)

Run the default MEND configuration and plot the results.

# Set up a time vector 
t <- seq(from = 1, to = 1e4, by = 10)

# Solve the MEND_model configuration, unless a new parameter table or vector of initial state values are
# specified as arguments the MEND_model configuration will use the default_params and default_inital inputs. 
out1 <- solve_model(mod = MEND_model, time =  t)

Plot results.

ggplot(data = out1) + 
  geom_line(aes(time, value, color = name)) + 
  labs(title = "MEND Output", 
       y = unique(out1$units)) + 
  facet_wrap("variable", scales = "free") + 
  theme(legend.position = 'none') 

Run MEND with a different parameter value

Here is the parameter table of the default values from Wang et. al 2013.

kable(default_params)

For this example let's double the half-saturation constant for decomposition of M (K.m) effects the rate of POM update by microbial biomass.

# Extract a the default value for the km parameter
value <- default_params[default_params$parameter == "K.m", "value"]

# Double default value 
new_km <- value * 2
names(new_km) <- "K.m"

# Make a new parameter table. 
new_param_table <- update_params(new_params = new_km, param_table = default_params)

# Use the new parameter table to solve the model. 
out2 <- solve_model(mod = MEND_model, time = t, params = new_param_table)

# Add identifying information to the output tables. 
out1$params <- "default Km"
out2$params <- "doubble Km"

rbind(out1, out2) %>%  
  ggplot(aes(time, value, color = params)) + 
  geom_line() + 
  facet_wrap("variable", scales = "free") + 
    labs(title = "MEND Output", 
       y = unique(out1$units)) + 
  facet_wrap("variable", scales = "free")

Run MEND using different inital conditions

# save a copy of the default initial conditions to manipulate. 
initial_cond <- default_inital

# Double the microbial bio mass 
initial_cond[["M"]] <- initial_cond[["M"]] * 2

# Solve MEND with the default parameter values, but using the new inital conditions
out3 <- solve_model(mod = MEND_model, time = t, params = default_params, state = initial_cond)

out1$inital <- "default inital"
out3$inital <- "double M inital"

rbind(out1, out3, fill = TRUE) %>% 
    ggplot(aes(time, value, color = inital)) + 
  geom_line() + 
  facet_wrap("variable", scales = "free") + 
    labs(title = "MEND Output", 
       y = unique(out1$units)) + 
  facet_wrap("variable", scales = "free")

kdorheim/MEMC documentation built on Jan. 14, 2022, 9:13 a.m.