-
Notifications
You must be signed in to change notification settings - Fork 1
/
Dynamic_RShell.Rmd
149 lines (108 loc) · 6.05 KB
/
Dynamic_RShell.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
---
title: "RShell for Quasi-dynamic (level IV) calculations using SimpleBox"
author: "Joris Quik"
date: "3/11/2021"
output: github_document
editor_options:
chunk_output_type: console
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(openxlsx)
library(deSolve)
```
## Introduction
This document describes the code and functions to produce quasi-dynamic ('levelIV') solutions of the [SimpleBox](https://www.rivm.nl/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.
```{r read in data, echo=TRUE}
#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")
# 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,sheet="dynamicR",colNames=FALSE, rows=15,cols=c(7:39)) # Emission rates for first period
SB.e2 <- read.xlsx(sb4n.loc,sheet="dynamicR",colNames=FALSE, rows=16,cols=c(7:39)) # Emission rates for second period
SB.tend1 <- as.numeric(read.xlsx(sb4n.loc,sheet="dynamicR",colNames=FALSE, rows=11,cols=6)) # Emission rates for second period
SB.tend2 <- as.numeric(read.xlsx(sb4n.loc,sheet="dynamicR",colNames=FALSE, rows=11,cols=7)) # Emission rates for second period
# Somehow these namedRegion's are not read in anymore (openxlsx bug?)
# 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: `r sb4n.loc`.
## 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:
```{r SimpleBoxODE}
# 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.
```{r calculate raw output}
# 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.
```{r}
# 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=as.data.frame(Mt), file = "data/SimpleBox_Dynamic_data.xlsx", sheetName = "dynamicR_data", asTable=TRUE)
```
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 *`r paste(sb4n.loc)`* file.