Skip to content

Latest commit

 

History

History
147 lines (114 loc) · 5.74 KB

Dynamic_RShell.md

File metadata and controls

147 lines (114 loc) · 5.74 KB

RShell for Quasi-dynamic (level IV) calculations using SimpleBox

Joris Quik 3/11/2021

Introduction

This document describes the code and functions to produce quasi-dynamic (‘levelIV’) solutions of the SimpleBox multimedia fate model.

Prerequisites to run this code is a recent version of SimpleBox (Feb. 2020 or later).

This script should work for both SimpleBox4.0 (conventional) and SimpleBox4nano.

preparing data

Define the location of the simplebox xlsx file. If used in an R project, the default location is the data folder within the project folder.

#For purpose of this example the latest SimpleBox file is downloaded to the data directory:
download.file(url = "https://github.com/rivm-syso/SimpleBox/archive/refs/heads/xl_version.zip"
                                   , destfile = "data/SimpleBox_xl.zip")
## Warning in download.file(url = "https://github.com/rivm-syso/SimpleBox/archive/
## refs/heads/xl_version.zip", : URL https://github.com/rivm-syso/SimpleBox/
## archive/refs/heads/xl_version.zip: cannot open destfile 'data/SimpleBox_xl.zip',
## reason 'Device or resource busy'

## Warning in download.file(url = "https://github.com/rivm-syso/SimpleBox/archive/
## refs/heads/xl_version.zip", : download had nonzero exit status

# unzip the SimpleBox xl file to the data directory
unzip(zipfile = "data/SimpleBox_xl.zip",files = "SimpleBox-xl_version/SimpleBox4.0_web.xlsm",
      junkpaths=TRUE, exdir = "data")

sb4n.loc <- "data/SimpleBox4.0_web.xlsm"

## NOTE ##
## Sometimes the naming of the cells in excel is different based on the version you are using.
## i.e. check that matrix "K" is the correct name, see worksheet "engine". It could be "Kmatrix"

SB.K <- as.matrix(read.xlsx(sb4n.loc,colNames=FALSE, namedRegion ="K")) # matrix of rate constants "k"
SB.m0 <- read.xlsx(sb4n.loc,colNames=FALSE, namedRegion ="m0") # Initial mass of each compartment, could be "m0"
SB.names <- read.xlsx(sb4n.loc,colNames=FALSE, namedRegion ="box_names") #Names for each compartment
SB.v <- read.xlsx(sb4n.loc,colNames=FALSE, namedRegion ="v") #Volumes
SB.e1 <- read.xlsx(sb4n.loc,colNames=FALSE, namedRegion ="emis.I") # Emission rates for first period
SB.e2 <- read.xlsx(sb4n.loc,colNames=FALSE, namedRegion ="emis.II") # Emission rates for second period
SB.tend1 <- as.numeric(read.xlsx(sb4n.loc,colNames=FALSE, namedRegion ="tend.I")) # Emission rates for second period
SB.tend2 <- as.numeric(read.xlsx(sb4n.loc,colNames=FALSE, namedRegion ="tend.II")) # Emission rates for second period
#ChemClass <- as.character(read.xlsx(sb4n.loc,colNames=FALSE, namedRegion ="ChemClass")) 
colnames(SB.v) <- SB.names
colnames(SB.e1) <- SB.names
colnames(SB.e2) <- SB.names

The K matrix, starting masses at t=0, emissions vector, the compartment names and compartment volumes are read from the SimpleBox xls file: data/SimpleBox4.0_web.xlsm.

Functions

The inverse of the matrix of rate constants (K) multiplied with the mass in each compartment (m) plus the emission (e) gives the change in mass of each time step per compartment. This is the main function:

# ODE function of SimpleBox:
SimpleBoxODE <- function(t, m, parms) {
  dm <- with(parms, K %*% m + e)
  return(list(dm))
}

An helper function is applied to calculate the dynamic output for separate periods. These are then added together to get the raw output which is a matrix of masses in each compartment per time point.

# helper function to calculate dynamic output for set time span
Rdyn.base <- function(Parms,    # list of parms needed for SimpleBoxODE function
                      tstart,   # First time point (in seconds)
                      tend,     # Last time point (in seconds)
                      SB.m0,    # Vector of initial masses
                      SB.names, # Vector of compartment names
                      tlength){ # Number of concentrations to output between start and end
  
  tset <- seq(from=tstart, to=tend, length.out=tlength) # make time sequence for which output is wanted
  Rpar.SBdyn <- list(tset=tset, m=as.numeric(SB.m0), parms=Parms)
  
  # Run ODE function to create out1 with masses in each compartment at different time points
  out <- with(Rpar.SBdyn, {
    ode(y=m,times=tset,func=SimpleBoxODE,parms=parms)
  })
  colnames(out)<-c("Time",SB.names)
  
  out
}

Calculation

Calculate the output.

# calculate mass in each compartment for 2 time spans with different emission regimes 
# period 1: 100% emission as input in sb4n.xlsx
# period 2: 0% emissions as input in sb4n.xlsx

Parms1 <- list(e=SB.e1, K=as.matrix(SB.K))
out.t1 <- Rdyn.base(Parms = Parms1,
                    tstart = 0, 
                    tend = SB.tend1*60*60*24*365,
                    SB.m0=SB.m0,
                    SB.names,
                    tlength = 50)

# period 2: 50% emission
Parms2 <- list(e=SB.e2, K=as.matrix(SB.K))
out.t2 <-  Rdyn.base(Parms = Parms2,
                     tstart = SB.tend1*60*60*24*365,
                     tend = SB.tend2*60*60*24*365,
                     SB.m0=out.t1[50,-1],
                     SB.names,
                     tlength = 50)
invMinf <- -as.vector(solve(SB.K)%*%as.numeric(SB.e1))


Mt1 <- cbind(out.t1[,1]/(3600*24*365),t(t(out.t1[,2:(length(out.t1[1,]))])/invMinf))
Mt2 <- cbind(out.t2[,1]/(3600*24*365),t(t(out.t2[,2:(length(out.t2[1,]))])/invMinf))


Mt <- rbind(Mt1,Mt2)
colnames(Mt) <- c("Time",SB.names)

write.xlsx(x=Mt, file = "data/SimpleBox_Dynamic_data.xlsx", sheetName = "dynamicR_data", colNames = TRUE,startRow=1,startCol = 1)

The output is stored in data/SimpleBox_Dynamic_data.xlsx. The output is the fraction of steady state in a compartment at each time step. Please copy the output to the dynamicR sheet in the data/SimpleBox4.0_web.xlsm file.