R

Use R function to generate the random sample based on the given sampling method. We assume the normal distribution for Vdist and kelim, and uniform distribution for gutabs.

Vdist_dist <- rnorm(n = 1000, mean = Vdist, sd = Vdist * 0.2) # CV = 0.2
kelim_dist <- rnorm(n = 1000, mean = kelim, sd = kelim * 0.2) # CV = 0.2
Fgutabs_dist <- runif(n = 1000, min = 0.8, max = 1.0)
css_dist <- dose * Fgutabs_dist / kelim_dist / Vdist_dist 
summary(css_dist)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.716   1.249   1.503   1.622   1.851   6.250

Visualize the parameter distribution through histogram.

par(mfrow = c(2,2))
hist(Vdist_dist)
hist(kelim_dist)
hist(Fgutabs_dist)
hist(css_dist)

MCSim

Run MC simulation by using model (“pbtk1cpt_css.model.R”) and input (“pbtk1cpt_css_mtc.in.R”) files.

out <- mcsim(model = "pbtk1cpt_css.model.R", input = "pbtk1cpt_css_mtc.in.R", dir = "modeling/pbtk1cpt")

## * Created executable program 'mcsim.pbtk1cpt_css.model.R.exe'.

## Execute: ./mcsim.pbtk1cpt_css.model.R.exe modeling/pbtk1cpt/pbtk1cpt_css_mtc.in.R

head(out)

##   Iter  Fgutabs     kelim   Vdist css_1.1
## 1    0 0.998091 0.0311610 7.01874 1.04181
## 2    1 0.940106 0.0327880 4.98724 1.31247
## 3    2 0.949883 0.0284472 3.82896 1.99085
## 4    3 0.848450 0.0229008 1.81946 4.64861
## 5    4 0.824792 0.0283581 5.67027 1.17099
## 6    5 0.823751 0.0230576 7.01450 1.16272

summary(out$css_1.1)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.7362  1.2525  1.5012  1.6189  1.8686  4.6486

Visualize the parameter distribution through histogram.

par(mfrow = c(2,2))
hist(out$Vdist)
hist(out$kelim)
hist(out$Fgutabs)
hist(out$css_1.1)

Compare the results from R and MCSim

# The example of using pipeline
# plot(density(css_dist))
css_dist %>% density() %>% plot()
out$css_1.1 %>% density() %>% lines(col = "red")

R

Generate the lognormal distribution of Css

Css <- rlnorm(n = 1000, meanlog = log(0.18/1000), sdlog = log(7)) # μg/L to mg/l
exp(sd(log(Css))) # Geometric Standard deviation

## [1] 7.229771

oral_equiv_dist <- Css * kelim_dist * Vdist_dist / Fgutabs_dist # mg/kg-d
summary(oral_equiv_dist)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 7.100e-08 7.085e-06 2.770e-05 1.914e-04 1.031e-04 1.441e-02

Make a histogram for oral equivalent distribution.

hist(oral_equiv_dist, main = "Predicted oral equivalent dose (mg/kg-d)")

Make a box plot

boxplot(oral_equiv_dist, log = "y", main = "Predicted oral equivalent dose (mg/kg-d)")

Plot the relationship between Css and oral equivalent dose. The detection rate and arithmetic mean are used to compare with “real” study data.

plot(Css, oral_equiv_dist, log = "xy", xlab = "Css (uM)", ylab = "oral equivalent dose (mg/kg-d)")
abline(v = 0.31/1000) # Add line of "real" detection limit

x <- subset(Css, Css > 0.31/1000)

length(x)/1000 # detection rate

## [1] 0.402

mean(x) * 1000 # arithmetic mean

## [1] 2.955248

MCSim

In this part, we will apply the same approach through MCSim. First, use mcsim() for the model (“pbtk1cpt_rtk.model.R”) and input (“pbtk1cpt_rtk.in.R”) files. Use set.seed() to create reproducible result

set.seed(1234)
out <- mcsim(model = "pbtk1cpt_rtk.model.R", input = "pbtk1cpt_rtk.in.R", dir = "modeling/pbtk1cpt")

## * Created executable program 'mcsim.pbtk1cpt_rtk.model.R.exe'.

## Execute: ./mcsim.pbtk1cpt_rtk.model.R.exe modeling/pbtk1cpt/pbtk1cpt_rtk.in.R

summary(out$Dose_1.1)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 1.200e-08 6.057e-06 2.323e-05 1.950e-04 8.852e-05 1.857e-02

Then, make histograms of Css and oral equivalent does

par(mfrow = c(1,2))
hist(out$Css)
hist(out$Dose_1.1)

Plot the relationship between Css and oral equivalent dose. The detection rate and arithmetic mean are used to compare with “real” study data.

plot(out$Css, out$Dose_1.1, log = "xy", xlab = "Css (uM)", ylab = "oral equivalent dose (mg/kg-d)")
abline(v = 0.31/1000) # detection limit
x <- subset(out$Css, out$Css > 0.31/1000)
DR_13 <- quantile(x, prob = 0.13)
DR_13 * 1000 # ug/L

##       13% 
## 0.4012078

abline(v = DR_13, col = "red") # 13 % Detection rate

length(x)/1000 # Detection rate

## [1] 0.359

mean(x) * 1000 # arithmetic mean

## [1] 3.323244

Compare the simulations from R and MCSim

boxplot(out$Dose_1.1, oral_equiv_dist, log = "y", names = c("R","MCSim"))

The variation of 10 different simulations

for (i in 1:10)
{
  out <- mcsim("pbtk1cpt_rtk.model.R", "pbtk1cpt_rtk.in.R", dir = "modeling/pbtk1cpt")
  x <- subset(out$Css, out$Css > 0.31/1000)
  Detect.rate <- length(x)/1000 # Detection rate
  Observ.mean <- mean(x) * 1000 # arithmetic mean
  plot(out$Css, out$Dose_1.1, log = "xy", 
       xlab = "Css (uM)", ylab = "oral equivalent dose (mg/kg-d)",
       main = paste("Detection rate: ", round(Detect.rate, digit = 2), 
                    "Observation mean: ", round(Observ.mean, digit = 2)))
  abline(v = 0.31/1000) # detection limit
  DL <- quantile(x, prob = 0.13)
  abline(v = DL, col = "red") # 13 %
  date_time<-Sys.time()
  while((as.numeric(Sys.time()) - as.numeric(date_time))<1){} 
}

eqn_1 <- function (X)
{
  X[, 1] + X[, 2] + X[, 3] # X1 + X2 + X3
}
eqn_2 <- function (X) 
{
  X[, 1] * X[, 2] / X[, 3] # X1 * X2 / X3
}

Morris method for eqn_1()

The Morris method,based on design of experiments, allows to identify the few important factors at a cost of \(r(p+ 1)\) simulations (where \(p\) is the number of factors). \(r\) is the number of repetitions of the design.

Using morris() and setting the arguments. The distribution for X1, X2, and X3 are U(0,1), U(1,4), and U(1,7), respectively.

set.seed(1234)
x <- morris(model = eqn_1, factors = 3, r = 18, 
            design = list(type = "oat", levels = 6, grid.jump = 3), 
            binf = c(0, 1, 1), bsup = c(1, 4, 7))
nrow(x$X) # Model evaluation

## [1] 72

#head(x$X)

Plotting sampling process

par(mfrow = c(3,1), mar = c(2,4,1,1))
plot(x$X[, 1], xlab = "n", ylab = "x1")
plot(x$X[, 2], ylab = "x2")
plot(x$X[, 3], ylab = "x3")

Parameter distributions

par(mfrow = c(1,3))
for(i in 1:3){
  hist(x$X[,i], main = colnames(x$X)[i])
}

Visualize the relationship between parameters and model outputs

par(mfrow = c(1,3))
for(i in 1:3){
  cor <- cor(x$X[,i], x$y)
  plot(x$X[,i], x$y, 
       xlab = colnames(x$X)[i],
       main = paste("r = ", round(cor, digits = 2)))
}

The \(\mu^*\) is the mean absolute value of elementary effect, representing the actual value from the output that was effected by specific parameter. The \(\sigma\) is the standard deviation of elementary effect, which can also represent the interaction of the parameters.

plot(x)

Report

x

## 
## Call:
## morris(model = eqn_1, factors = 3, r = 18, design = list(type = "oat",     levels = 6, grid.jump = 3), binf = c(0, 1, 1), bsup = c(1,     4, 7))
## 
## Model runs: 72 
##    mu mu.star        sigma
## X1  1       1 1.355467e-15
## X2  3       3 9.871551e-16
## X3  6       6 1.160045e-15

Analytical solution (first order)

1-0 # X1

## [1] 1

4-1 # X2

## [1] 3

7-1 # X3

## [1] 6

FAST for eqn_1()

Using fast99() and setting the arguments. This method allows the estimation of first order and total Sobol indices for all the factors (alltogether \(2p\) indices, where \(p\) is the number of factors) at a total cost of \(n \times p\) simulations.

Set parameter distribution

q <- "qunif"
q.arg <- list(list(min = 0,  max = 1),
              list(min = 1,  max = 4),
              list(min = 1,  max = 7))

Using fast99() and setting the arguments

x <- fast99(model = eqn_1, factors = 3, n = 128, # Need n = 1000 to converge 
            q = q, q.arg = q.arg)

Plotting sampling process

par(mfrow = c(3,1), mar = c(2,2,1,1))
plot(x$X[,1], type = "l")
plot(x$X[,2], type = "l")
plot(x$X[,3], type = "l")

Parameter distributions

par(mfrow = c(1,3))
for(i in 1:3){
  hist(x$X[,i], main = colnames(x$X)[i])
}

Visualize the relationship between parameters and model outputs

par(mfrow = c(1,3))
for(i in 1:3){
  cor <- cor(x$X[,i], x$y)
  plot(x$X[,i], x$y, 
       xlab = colnames(x$X)[i],
       main = paste("r = ", round(cor, digits = 2)))
}

Plot the sensitivity measurements

plot(x)

Report

x

## 
## Call:
## fast99(model = eqn_1, factors = 3, n = 128, q = q, q.arg = q.arg)
## 
## Model runs: 384 
## 
## Estimations of the indices:
##    first order total order
## X1  0.01118566   0.0114904
## X2  0.15241837   0.1529183
## X3  0.69199565   0.6934759

Analytical solution

Total <- (1-0)^2 + (4-1)^2 + (7-1)^2

1^2 / Total # S1

## [1] 0.02173913

3^2 / Total # S2

## [1] 0.1956522

6^2 / Total # S3

## [1] 0.7826087

Morris method for eqn_2()

Call morris() and use the same parameter setting as previous case.

set.seed(1234)
x <- morris(model = eqn_2, factors = 3, r = 18, 
            design = list(type = "oat", levels = 6, grid.jump = 3),
            binf = c(0,1,1), bsup = c(1,4,7))
nrow(x$X) # Model evaluation

## [1] 72

Visualize the relationship between parameters and model outputs

par(mfrow = c(1,3))
for(i in 1:3){
  cor <- cor(x$X[,i], x$y)
  plot(x$X[,i], x$y, 
       xlab = colnames(x$X)[i],
       main = paste("r = ", round(cor, digits = 2)))
}

Plot the sensitivity measurements

plot(x, xlim = c(0,2), ylim = c(0,2))

Report

x

## 
## Call:
## morris(model = eqn_2, factors = 3, r = 18, design = list(type = "oat",     levels = 6, grid.jump = 3), binf = c(0, 1, 1), bsup = c(1,     4, 7))
## 
## Model runs: 72 
##            mu   mu.star     sigma
## X1  1.0008280 1.0008280 0.8454144
## X2  0.5051922 0.5051922 0.7424151
## X3 -0.8651525 0.8651525 1.2190095

FAST method for eqn_2()

q <- "qunif"
q.arg <- list(list(min = 0,  max = 1),
              list(min = 1,  max = 4),
              list(min = 1,  max = 7))

Call fast99()

x <- fast99(model = eqn_2, factors = 3, n = 1000, q = q, q.arg = q.arg) #

Visualize the relationship between parameters and model outputs

par(mfrow = c(1,3))
for(i in 1:3){
  cor <- cor(x$X[,i], x$y)
  plot(x$X[,i], x$y, 
       xlab = colnames(x$X)[i],
       main = paste("r = ", round(cor, digits = 2)))
}

Plot the sensitivity measurements

plot(x)

Report

x

## 
## Call:
## fast99(model = eqn_2, factors = 3, n = 1000, q = q, q.arg = q.arg)
## 
## Model runs: 3000 
## 
## Estimations of the indices:
##    first order total order
## X1   0.3220563   0.5063070
## X2   0.1160247   0.2293264
## X3   0.3396511   0.5139026

Morris method of Css

Define the function to estimate the maximum concentration at steady-state (Css)

Css_fun <- function (X) 
{
  Dose <- 0.228 # Ingestion dose (uM)
  Dose * X[, 1] / X[, 2] / X[, 3] # Fgutabs_dist / kelim_dist / Vdist_dist 
}

Set parameter range,

# The order of binf and bsup are c(Fgutabs, kelim, Vdist)
binf <- c(0.8, kelim * 0.5, Vdist * 0.5)
bsup <- c(1, kelim * 2, Vdist * 2)

Call morris()

set.seed(1234)
x <- morris(model = Css_fun, factors = c("Fgutabs", "kelim", "Vdist"), r = 32, 
            design = list(type = "oat", levels = 6, grid.jump = 3), 
            binf = binf, bsup = bsup)

Visualize the relationship between parameters and model outputs

par(mfrow = c(1,3))
for(i in 1:3){
  cor <- cor(x$X[,i], x$y)
  plot(x$X[,i], x$y, 
       xlab = colnames(x$X)[i],
       main = paste("r = ", round(cor, digits = 2)))
}

Plot the sensitivity measurements

par(mfrow = c(1,1))
plot(x, xlim = c(0, 5), ylim = c(0, 5))
abline(0,1) # non-linear and/or non-monotonic
abline(0,0.5, lty = 2) # monotonic
abline(0,0.1, lty = 3) # almost linear
legend("topleft", legend = c("non-linear and/or non-monotonic",
                             "monotonic", "linear"), lty = c(1:3))

The stability testing

for (i in 1:10)
{
  x <- morris(model = Css_fun, factors = c("Fgutabs", "kelim", "Vdist"), r = 32, # test r = 32
              design = list(type = "oat", levels = 6, grid.jump = 3), 
              binf = binf, bsup = bsup)
  plot(x, xlim = c(0, 6), ylim = c(0, 6))
  abline(0,1) # non-linear and/or non-monotonic
  abline(0,0.5, lty = 2) # monotonic
  abline(0,0.1, lty = 3) # almost linear
  legend("topleft", legend = c("non-linear and/or non-monotonic",
                               "monotonic", "linear"), lty = c(1:3))
  date_time<-Sys.time()
  while((as.numeric(Sys.time()) - as.numeric(date_time))<1){} 
}

The example of convergence diagnostics

for (i in 3:11){
  x <- morris(model = Css_fun, factors = c("Fgutabs", "kelim", "Vdist"), r = 2^i, # test r = 32
              design = list(type = "oat", levels = 6, grid.jump = 3), 
              binf = binf, bsup = bsup)
  if (i == 3){ X <- apply(abs(x$ee), 2, mean) } else X <- rbind(X, apply(abs(x$ee), 2, mean))
}
X

##     Fgutabs     kelim    Vdist
## X 0.2240959 1.1615334 1.143272
##   0.2137596 0.9145757 1.374492
##   0.3056444 2.0466927 2.010046
##   0.3404735 1.8931413 1.970940
##   0.3028373 1.9423194 1.924475
##   0.3253824 1.9860593 1.988660
##   0.3041595 1.8543542 1.867326
##   0.3212376 1.8882862 1.942005
##   0.3243742 1.9074727 1.896238

ylim <- range(X)
plot(X[,1], ylim = ylim, type = "b", xaxt = "n")
axis(1, at = seq(1,9,1), labels = 2^seq(3,11,1))
lines(X[,2], type = "b")
lines(X[,3], type = "b")

FAST method of Css

Set parameter distribution

q <- "qunif"
q.arg <- list(list(min = 0.8,  max = 1.0),
              list(min = kelim * 0.5,  max = kelim * 2),
              list(min = Vdist * 0.5,  max = Vdist * 2))

Call fast99()

x <- fast99(model = Css_fun, factors = c("Fgutabs", "kelim", "Vdist"), 
            n = 500, q = q, q.arg = q.arg)

Visualize the relationship between parameters and model outputs

par(mfrow = c(1,3))
for(i in 1:3){
  cor <- cor(x$X[,i], x$y)
  plot(x$X[,i], x$y, 
       xlab = colnames(x$X)[i],
       main = paste("r = ", round(cor, digits = 2)))
}

Plot the sensitivity measurements

plot(x)

Report

x

## 
## Call:
## fast99(model = Css_fun, factors = c("Fgutabs", "kelim", "Vdist"),     n = 500, q = q, q.arg = q.arg)
## 
## Model runs: 1500 
## 
## Estimations of the indices:
##         first order total order
## Fgutabs  0.01211322  0.02427179
## kelim    0.44375555  0.53829078
## Vdist    0.44375555  0.53829078

The example of convergence diagnostics

for (i in 1:8){
  x <- fast99(model = Css_fun, factors = c("Fgutabs", "kelim", "Vdist"), n = 2^(i+6),
            q = q, q.arg = q.arg)
  if (i == 1){ X <- 1 - x$Dt/x$V } else X <- rbind(X, 1 - x$Dt/x$V)
}
ylim <- range(X)
X

##         [,1]      [,2]      [,3]
## X 0.03065326 0.6140697 0.6140697
##   0.02860018 0.5356477 0.5356477
##   0.01497975 0.5571496 0.5571496
##   0.02325690 0.5420919 0.5420919
##   0.02366698 0.5397771 0.5397771
##   0.02330789 0.5394178 0.5394178
##   0.02348126 0.5394023 0.5394023
##   0.02352585 0.5393984 0.5393984

plot(X[,1], ylim = ylim, type = "b", xaxt = "n")
axis(1, at = seq(1,8,1), labels = 2^seq(7,14,1))
lines(X[,2], type = "b")
lines(X[,3], type = "b")

Setpoints method (Morris)

Call Morris()

# The order of binf and bsup are c("Vdist", "kelim", "kgutabs", "Fgutabs")
binf <- c(Vdist * 0.5, kelim * 0.5, kgutabs * 0.5, 0.8)
bsup <- c(Vdist * 2, kelim * 2, kgutabs * 2, 1.0)
set.seed(1234)
x <- morris(model = NULL, factors = c("Vdist", "kelim", "kgutabs", "Fgutabs"), r = 32, 
            design = list(type = "oat", levels = 6, grid.jump = 3), 
            binf = binf, bsup = bsup)
head(x$X)

##          Vdist      kelim kgutabs Fgutabs
## [1,]  4.909793 0.01826586   2.398    0.84
## [2,] 10.433310 0.01826586   2.398    0.84
## [3,] 10.433310 0.01826586   4.360    0.84
## [4,] 10.433310 0.01826586   4.360    0.96
## [5,] 10.433310 0.03881496   4.360    0.96
## [6,]  4.909793 0.01141616   4.360    0.88

nrow(x$X)

## [1] 160

Create the setpoints file setpts.out

X <- cbind(1, x$X)
write.table(X, file = "setpts.out", row.names = F, sep = "\t")

Modeling

out <- mcsim("pbtk1cpt.model.R", "pbtk1cpt_setpts.in.R", dir = "modeling/pbtk1cpt")

## Execute: ./mcsim.pbtk1cpt.model.R.exe modeling/pbtk1cpt/pbtk1cpt_setpts.in.R

head(out)

##   Iter    Vdist     kelim kgutabs Fgutabs Ccompartment_1.1
## 1    0  4.90979 0.0182659   2.398    0.84                0
## 2    1 10.43330 0.0182659   2.398    0.84                0
## 3    2 10.43330 0.0182659   4.360    0.84                0
## 4    3 10.43330 0.0182659   4.360    0.96                0
## 5    4 10.43330 0.0388150   4.360    0.96                0
## 6    5  4.90979 0.0114162   4.360    0.88                0
##   Ccompartment_1.2 Ccompartment_1.3 Ccompartment_1.4 Ccompartment_1.5
## 1         0.672862         0.721846         0.714340         0.701916
## 2         0.316641         0.339692         0.336160         0.330313
## 3         0.343218         0.341391         0.335268         0.329201
## 4         0.392249         0.390162         0.383164         0.376229
## 5         0.385993         0.376230         0.361969         0.348189
## 6         0.768177         0.769274         0.760667         0.752034
##   Ccompartment_1.6 Ccompartment_1.7 Ccompartment_1.8 Ccompartment_1.9
## 1         0.689257         0.676786         0.664536         0.652508
## 2         0.324356         0.318487         0.312723         0.307063
## 3         0.323242         0.317391         0.311647         0.306006
## 4         0.369419         0.362733         0.356167         0.349721
## 5         0.334933         0.322182         0.309916         0.298117
## 6         0.743497         0.735058         0.726714         0.718465
##   Ccompartment_1.10 Ccompartment_1.11 Ccompartment_1.12 Ccompartment_1.13
## 1          0.640698          0.629101          0.617714          0.606534
## 2          0.301505          0.296048          0.290689          0.285428
## 3          0.300467          0.295029          0.289688          0.284445
## 4          0.343391          0.337175          0.331073          0.325080
## 5          0.286767          0.275850          0.265348          0.255246
## 6          0.710309          0.702246          0.694275          0.686394
##   Ccompartment_1.14 Ccompartment_1.15 Ccompartment_1.16 Ccompartment_1.17
## 1          0.595555          0.584776          0.574191          0.563798
## 2          0.280261          0.275189          0.270208          0.265317
## 3          0.279297          0.274241          0.269278          0.264404
## 4          0.319196          0.313419          0.307746          0.302176
## 5          0.245528          0.236181          0.227189          0.218540
## 6          0.678603          0.670900          0.663284          0.655755
##   Ccompartment_1.18 Ccompartment_1.19 Ccompartment_1.20 Ccompartment_1.21
## 1          0.553594          0.543574          0.533735          0.524074
## 2          0.260515          0.255799          0.251169          0.246623
## 3          0.259618          0.254919          0.250305          0.245774
## 4          0.296706          0.291336          0.286063          0.280885
## 5          0.210219          0.202216          0.194517          0.187112
## 6          0.648311          0.640952          0.633677          0.626484
##   Ccompartment_1.22 Ccompartment_1.23 Ccompartment_1.24
## 1          0.514588          0.505274          0.496129
## 2          0.242159          0.237776          0.233472
## 3          0.241326          0.236958          0.232669
## 4          0.275801          0.270809          0.265907
## 5          0.179988          0.173136          0.166544
## 6          0.619372          0.612342          0.605391

Check the time-dependent sensitivity measures

{r, eval=TF for(i in 2:24){ tell(x, out[,i+5]) plot(x, main = paste("Time = ", i-1, "hr"), xlim = c(0, 1), ylim = c(0, 1)) abline(0,1) # non-linear and/or non-monotonic abline(0,0.5, lty = 2) # monotonic abline(0,0.1, lty = 3) # almost linear date_time<-Sys.time() while((as.numeric(Sys.time()) - as.numeric(date_time))<0.5){} }

Data manipulate (Transfer to tidydata format)

index <- which(names(out) == "Ccompartment_1.1" | names(out) == "Ccompartment_1.24")
X <- apply(out[,index[1]:index[2]], 2, quantile,  c(0.5, 0.025, 0.975))
dat <- t(X)
colnames(dat) <- c("median", "LCL", "UCL")
df <- as.data.frame(dat)
df$time <- seq(0, 23, 1)
head(df)

##                     median       LCL      UCL time
## Ccompartment_1.1 0.0000000 0.0000000 0.000000    0
## Ccompartment_1.2 0.4717685 0.2386066 1.237713    1
## Ccompartment_1.3 0.4795500 0.2719481 1.319638    2
## Ccompartment_1.4 0.4759115 0.2715558 1.306399    3
## Ccompartment_1.5 0.4696055 0.2654703 1.279288    4
## Ccompartment_1.6 0.4575730 0.2544477 1.256138    5

Visualization

ggplot(df, aes(x = time, y = median)) +
  geom_ribbon(aes(ymin = LCL, ymax = UCL), fill = "grey70", alpha = 0.5) + 
  geom_line()