Tutorial 2: Uncertainty and Sensitivity Analysis (Supplementary)
Nan-Hung Hsieh
2019/05/16 (update: 2019-05-29)
4 Example 2: PERC-PBPK model
The list of R packages should be installed first to do the following testing.
library(tidyverse)
library(sensitivity)
library(pksensi)4.1 Monte Carlo simulation in PERC-PBPK model
Prepare the model (“perc.model.R”) and input (“perc_mtc.in.R”) files. There are two simulations, which are inhalation of 72 and 144 ppm with 4 hour exposure.
Modeling
out <- mcsim(model = "perc.model.R", input = "perc_mtc.in.R", dir = "modeling/perc")## * Created executable program 'mcsim.perc.model.R.exe'.## Execute: ./mcsim.perc.model.R.exe modeling/perc/perc_mtc.in.R- Data manipulate
vars <- names(out)
sim1.1 <- which(vars == "C_exh_ug_1.1" | vars == "C_exh_ug_1.10")
sim1.2 <- which(vars == "C_ven_1.1" | vars == "C_ven_1.8")
sim2.1 <- which(vars == "C_exh_ug_2.1" | vars == "C_exh_ug_2.10")
sim2.2 <- which(vars == "C_ven_2.1" | vars == "C_ven_2.8")- Visualize
last.param <- which(names(out) == "Flow_pul")
par(mfrow = c(4, 4), mar = c(2,2,4,1))
for (i in 2:last.param){
hist(out[,i], main = names(out[i]), xlab = "")
}
ggPK <- function(data, index, time, ...){
X <- apply(data[,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 <- time
ggplot(df, aes(x = time, y = median)) +
geom_ribbon(aes(ymin = LCL, ymax = UCL), fill = "grey70", alpha = 0.5) +
geom_line() +
labs(...) +
scale_y_log10()
}t_exh_ug <- c(239.9, 245, 270, 360, 1320, 2760, 4260, 5700, 8580, 10020)
t_Ven <- c(239.9, 360, 1320, 2760, 4260, 5700, 8580, 10020)
ggPK(data = out, index = sim1.1, time = t_exh_ug, x = "Time (min)", y = "mg/l in the exhaled air")
ggPK(data = out, index = sim1.2, time = t_Ven, x = "Time (min)", y = "mg/l in the venous blood")
ggPK(data = out, index = sim2.1, time = t_exh_ug, x = "Time (min)", y = "mg/l in the exhaled air")
ggPK(data = out, index = sim2.2, time = t_Ven, x = "Time (min)", y = "mg/l in the venous blood")
4.2 Application of pksensi in Monte Carlo simulation
- Set parameter distribution (assign
parms,dist,q.qarg)
parameters <- c("LeanBodyWt","Pct_M_fat", "Pct_LM_liv", "Pct_LM_wp",
"Pct_Flow_fat", "Pct_Flow_liv", "Pct_Flow_pp",
"PC_fat", "PC_liv", "PC_wp", "PC_pp", "PC_art",
"Vent_Perf", "sc_Vmax", "Km", "Flow_pul")
dist <- rep("Uniform", 16) # MCSim definition
q.arg<-list(list(50, 70),
list(0.2, 0.3),
list(0.03, 0.04),
list(0.25, 0.3),
list(0.06, 0.08),
list(0.2, 0.3),
list(0.2, 0.25),
list(110, 150),
list(5.0, 8.0),
list(5.0, 8.5),
list(1.6, 1.8),
list(12, 15),
list(0.8, 1.3),
list(0.04, 0.06),
list(7, 13),
list(7.4, 7.6))- Set experiment time-points, output variables, and conditions
times <- c(239.9, 245, 270, 360, 1320, 2760, 4260, 5700, 8580, 10020)
outputs <- c("C_exh_ug", "C_ven")
conditions <- c("InhMag = 72", "Period = 1e10", "Exposure = 240")- Modeling
set.seed(2222)
y<-solve_mcsim(mName = "perc.model.R", params = parameters, vars = outputs,
monte_carlo = 1000, dist = dist, q.arg = q.arg, time = times,
condition = conditions)## Starting time: 2019-05-29 12:02:47## * Created input file "sim.in".## Execute: ./mcsim.perc.model.R.exe sim.in## Ending time: 2019-05-29 12:02:48- Visualization
par(mfrow = c(1, 2), mar = c(2,3,1,1))
pksim(y, vars = "C_exh_ug", log = T)
pksim(y, vars = "C_ven", log = T)
4.3 Application of pksensi in sensitivity analysis
- Set parameter distribution (assign
parms,dist,q.qarg)
parameters <- c("LeanBodyWt",
"Pct_M_fat", "Pct_LM_liv", "Pct_LM_wp","Pct_Flow_fat", "Pct_Flow_liv", "Pct_Flow_pp",
"PC_fat", "PC_liv", "PC_wp", "PC_pp", "PC_art",
"Vent_Perf",
"sc_Vmax", "Km",
"Flow_pul")
q <- "qunif"
q.arg <- list(list(50, 70),
list(0.2, 0.3),
list(0.03, 0.04),
list(0.25, 0.3),
list(0.06, 0.08),
list(0.2, 0.3),
list(0.2, 0.25),
list(110, 150),
list(5.0, 8.0),
list(5.0, 8.5),
list(1.6, 1.8),
list(12, 15),
list(0.8, 1.3),
list(0.04, 0.06),
list(7, 13),
list(7.4, 7.6))- Generate parameter space
x <- rfast99(params = parameters, n = 2000, q = q, q.arg = q.arg, replicate = 1)
dim(x$a)## [1] 32000 1 16- Set experiment time-points, output variables, and conditions
times <- c(239.9, 245, 270, 360, 1320, 2760, 4260, 5700, 8580, 10020)
outputs <- c("C_exh_ug", "C_ven")
conditions <- c("InhMag = 72", "Period = 1e10", "Exposure = 240")- Modeling
y <- solve_mcsim(x, mName = "perc.model.R", params = parameters, vars = outputs,
time = times, condition = conditions)## Starting time: 2019-05-29 12:02:48## * Created input file "sim.in".## Execute: ./mcsim.perc.model.R.exe sim.in## Ending time: 2019-05-29 12:02:57- Visualization & Decision
par(mfrow = c(4, 4), mar = c(2,2,4,1))
for(i in 1:16){
hist(x$a[,,i])
}
plot(y)
plot(y, vars = 2)
heat_check(y) 
check(y)##
## Sensitivity check ( Index > 0.05 )
## ----------------------------------
## First order:
## LeanBodyWt Pct_M_fat Pct_Flow_fat Pct_Flow_liv PC_fat PC_liv PC_wp PC_art Vent_Perf sc_Vmax Km
##
## Interaction:
## Vent_Perf
##
## Total order:
## LeanBodyWt Pct_M_fat Pct_Flow_fat Pct_Flow_liv PC_fat PC_liv PC_wp PC_art Vent_Perf sc_Vmax Km
##
## Unselected factors in total order:
## Pct_LM_liv Pct_LM_wp Pct_Flow_pp PC_pp Flow_pul
##
##
## Convergence check ( Index > 0.05 )
## ----------------------------------
## First order:
##
##
## Interaction:
##
##
## Total order:
## 5 Exercise
In the exercise 1 to 4, we’ll use the EB.model.R to conduct Monte Carlo simulation and sensitivity analysis.
5.1 Monte Carlo simulation
Create the MCSim’s input file (“EB_mtc.in.R”) and run Monte Carlo simulation to conduct uncertainty analysis of EB concentration in blood (0 - 6 hr), the exposure condition is 100 ppm EB exposure for 4 hours.
out <- mcsim(model = "EB.model.R", input = "EB_mtc.in.R", dir = "modeling/EB")## * Created executable program 'mcsim.EB.model.R.exe'.## Execute: ./mcsim.EB.model.R.exe modeling/EB/EB_mtc.in.Rhead(out)## Iter BW Pb Pl Pf Pm Pvrg Ppu Pbr
## 1 0 0.0536129 46.4148 1.33159 68.6502 0.866053 3.09945 1.49902 1.09595
## 2 1 0.0421327 85.2850 3.39734 25.1517 1.110290 2.06726 2.39587 1.22651
## 3 2 0.0413123 53.6559 1.61733 40.7100 0.353186 3.17564 3.11890 2.09262
## 4 3 0.0481246 36.1252 1.19474 51.8125 0.875638 1.14822 2.91475 1.96220
## 5 4 0.0541359 77.5347 1.29443 45.5223 0.638040 2.80191 1.69476 1.08215
## 6 5 0.0456809 66.3341 1.01652 60.6285 0.646264 3.73141 3.66144 3.61670
## VmaxC VmaxClu VmaxCvr Cvtot_1.1 Cvtot_1.2 Cvtot_1.3 Cvtot_1.4
## 1 22.32630 49.1016 32.4981 0 1.62551e-06 1.95642e-06 2.16516e-06
## 2 6.24365 41.0164 17.1240 0 3.14474e-06 3.78492e-06 4.26062e-06
## 3 13.26260 41.4444 27.0524 0 2.49120e-06 2.84300e-06 2.99543e-06
## 4 23.64120 12.1425 56.9746 0 1.21366e-06 1.51105e-06 1.69948e-06
## 5 24.74920 17.9539 38.6639 0 1.62319e-06 1.97096e-06 2.16689e-06
## 6 4.49791 18.9275 16.3190 0 3.80701e-06 4.51663e-06 4.92893e-06
## Cvtot_1.5 Cvtot_1.6 Cvtot_1.7 Cvtot_1.8 Cvtot_1.9 Cvtot_1.10
## 1 2.30293e-06 2.39916e-06 2.47074e-06 2.52729e-06 2.57441e-06 2.61532e-06
## 2 4.62151e-06 4.90094e-06 5.12133e-06 5.29812e-06 5.44205e-06 5.56067e-06
## 3 3.09379e-06 3.17359e-06 3.24352e-06 3.30624e-06 3.36284e-06 3.41400e-06
## 4 1.82517e-06 1.91424e-06 1.98137e-06 2.03488e-06 2.07957e-06 2.11822e-06
## 5 2.28903e-06 2.37434e-06 2.43984e-06 2.49420e-06 2.54142e-06 2.58346e-06
## 6 5.18987e-06 5.37258e-06 5.51348e-06 5.63085e-06 5.73386e-06 5.82718e-06
## Cvtot_1.11 Cvtot_1.12 Cvtot_1.13 Cvtot_1.14 Cvtot_1.15 Cvtot_1.16
## 1 2.65196e-06 2.68539e-06 2.71632e-06 2.74514e-06 2.77217e-06 2.79760e-06
## 2 5.65942e-06 5.74229e-06 5.81227e-06 5.87166e-06 5.92224e-06 5.96545e-06
## 3 3.46028e-06 3.50214e-06 3.53999e-06 3.57423e-06 3.60520e-06 3.63321e-06
## 4 2.15245e-06 2.18328e-06 2.21134e-06 2.23702e-06 2.26062e-06 2.28237e-06
## 5 2.62143e-06 2.65603e-06 2.68768e-06 2.71669e-06 2.74330e-06 2.76773e-06
## 6 5.91327e-06 5.99350e-06 6.06867e-06 6.13928e-06 6.20571e-06 6.26826e-06
## Cvtot_1.17 Cvtot_1.18 Cvtot_1.19 Cvtot_1.20 Cvtot_1.21 Cvtot_1.22
## 1 2.82155e-06 2.84415e-06 2.86549e-06 2.88564e-06 2.90468e-06 2.92267e-06
## 2 6.00244e-06 6.03415e-06 6.06137e-06 6.08477e-06 6.10488e-06 6.12218e-06
## 3 3.65854e-06 3.68145e-06 3.70217e-06 3.72090e-06 3.73785e-06 3.75317e-06
## 4 2.30244e-06 2.32098e-06 2.33811e-06 2.35396e-06 2.36862e-06 2.38217e-06
## 5 2.79016e-06 2.81076e-06 2.82968e-06 2.84705e-06 2.86301e-06 2.87766e-06
## 6 6.32716e-06 6.38264e-06 6.43491e-06 6.48415e-06 6.53053e-06 6.57422e-06
## Cvtot_1.23 Cvtot_1.24 Cvtot_1.25 Cvtot_1.26 Cvtot_1.27 Cvtot_1.28
## 1 2.93967e-06 2.95573e-06 2.97091e-06 2.98526e-06 2.99881e-06 3.01162e-06
## 2 6.13706e-06 6.14988e-06 6.16091e-06 6.17041e-06 6.17858e-06 6.18563e-06
## 3 3.76703e-06 3.77956e-06 3.79089e-06 3.80113e-06 3.81040e-06 3.81878e-06
## 4 2.39470e-06 2.40629e-06 2.41701e-06 2.42692e-06 2.43609e-06 2.44458e-06
## 5 2.89111e-06 2.90347e-06 2.91482e-06 2.92524e-06 2.93482e-06 2.94362e-06
## 6 6.61537e-06 6.65414e-06 6.69065e-06 6.72504e-06 6.75743e-06 6.78793e-06
## Cvtot_1.29 Cvtot_1.30 Cvtot_1.31 Cvtot_1.32 Cvtot_1.33 Cvtot_1.34
## 1 3.02373e-06 3.03517e-06 3.04598e-06 3.05620e-06 3.06586e-06 3.07498e-06
## 2 6.19169e-06 6.19692e-06 6.20141e-06 6.20529e-06 6.20862e-06 6.21150e-06
## 3 3.82635e-06 3.83320e-06 3.83940e-06 3.84500e-06 3.85006e-06 3.85464e-06
## 4 2.45243e-06 2.45969e-06 2.46641e-06 2.47262e-06 2.47836e-06 2.48368e-06
## 5 2.95170e-06 2.95911e-06 2.96592e-06 2.97217e-06 2.97791e-06 2.98318e-06
## 6 6.81665e-06 6.84371e-06 6.86918e-06 6.89317e-06 6.91576e-06 6.93703e-06
## Cvtot_1.35 Cvtot_1.36 Cvtot_1.37 Cvtot_1.38 Cvtot_1.39 Cvtot_1.40
## 1 3.08361e-06 3.09176e-06 3.09946e-06 3.10674e-06 3.11361e-06 3.12011e-06
## 2 6.21397e-06 6.21611e-06 6.21794e-06 6.21952e-06 6.22088e-06 6.22206e-06
## 3 3.85879e-06 3.86253e-06 3.86592e-06 3.86898e-06 3.87175e-06 3.87425e-06
## 4 2.48859e-06 2.49314e-06 2.49734e-06 2.50123e-06 2.50483e-06 2.50816e-06
## 5 2.98802e-06 2.99247e-06 2.99655e-06 3.00030e-06 3.00375e-06 3.00691e-06
## 6 6.95706e-06 6.97592e-06 6.99368e-06 7.01040e-06 7.02614e-06 7.04096e-06
## Cvtot_1.41 Cvtot_1.42 Cvtot_1.43 Cvtot_1.44 Cvtot_1.45 Cvtot_1.46
## 1 3.12626e-06 1.45767e-06 1.13800e-06 9.39173e-07 8.09987e-07 7.21233e-07
## 2 6.22307e-06 2.72125e-06 2.10441e-06 1.65902e-06 1.32958e-06 1.08046e-06
## 3 3.87651e-06 1.30261e-06 9.68009e-07 8.25445e-07 7.33882e-07 6.60114e-07
## 4 2.51124e-06 1.27761e-06 9.85733e-07 8.02095e-07 6.80405e-07 5.95159e-07
## 5 3.00981e-06 1.34927e-06 1.00983e-06 8.20488e-07 7.03562e-07 6.22679e-07
## 6 7.05492e-06 2.76990e-06 2.11032e-06 1.74768e-06 1.52851e-06 1.38040e-06
## Cvtot_1.47 Cvtot_1.48 Cvtot_1.49 Cvtot_1.50 Cvtot_1.51 Cvtot_1.52
## 1 6.56281e-07 6.05699e-07 5.64044e-07 5.28178e-07 4.96284e-07 4.67284e-07
## 2 8.88039e-07 7.36527e-07 6.15232e-07 5.16778e-07 4.35938e-07 3.68939e-07
## 3 5.95325e-07 5.37368e-07 4.85170e-07 4.38061e-07 3.95544e-07 3.57163e-07
## 4 5.31067e-07 4.80369e-07 4.38210e-07 4.01842e-07 3.69710e-07 3.40819e-07
## 5 5.60788e-07 5.09713e-07 4.65512e-07 4.26155e-07 3.90617e-07 3.58275e-07
## 6 1.26889e-06 1.17738e-06 1.09782e-06 1.02631e-06 9.60763e-07 9.00016e-07
## Cvtot_1.53 Cvtot_1.54 Cvtot_1.55 Cvtot_1.56 Cvtot_1.57 Cvtot_1.58
## 1 4.40533e-07 4.15629e-07 3.92315e-07 3.70412e-07 3.49791e-07 3.30351e-07
## 2 3.13006e-07 2.66049e-07 2.26455e-07 1.92955e-07 1.64541e-07 1.40394e-07
## 3 3.22509e-07 2.91222e-07 2.62972e-07 2.37466e-07 2.14435e-07 1.93639e-07
## 4 3.14570e-07 2.90550e-07 2.68484e-07 2.48162e-07 2.29416e-07 2.12109e-07
## 5 3.28719e-07 3.01652e-07 2.76839e-07 2.54079e-07 2.33196e-07 2.14033e-07
## 6 8.43414e-07 7.90527e-07 7.41041e-07 6.94701e-07 6.51289e-07 6.10611e-07
## Cvtot_1.59 Cvtot_1.60 Cvtot_1.61
## 1 3.12011e-07 2.94701e-07 2.78358e-07
## 2 1.19843e-07 1.02335e-07 8.74066e-08
## 3 1.74862e-07 1.57906e-07 1.42596e-07
## 4 1.96118e-07 1.81341e-07 1.67681e-07
## 5 1.96446e-07 1.80305e-07 1.65492e-07
## 6 5.72488e-07 5.36758e-07 5.03268e-07Data manipulate
vars <- names(out)
index <- which(vars == "Cvtot_1.2" | vars == "Cvtot_1.61")
time <- seq(0.1, 6, 0.1)
MW <- 106.16 # g/mol
X <- apply(out[,index[1]:index[2]], 2, quantile, c(0.5, 0.01, 0.99))
dat <- t(X) * MW * 1000
colnames(dat) <- c("median", "LCL", "UCL")
df <- as.data.frame(dat)
df$time <- time
head(df)## median LCL UCL time
## Cvtot_1.2 0.2097000 0.1217075 0.5732374 0.1
## Cvtot_1.3 0.2510530 0.1528508 0.6943341 0.2
## Cvtot_1.4 0.2745473 0.1722578 0.7726854 0.3
## Cvtot_1.5 0.2902409 0.1858845 0.8295683 0.4
## Cvtot_1.6 0.3023230 0.1955990 0.8751173 0.5
## Cvtot_1.7 0.3118078 0.2029496 0.9226145 0.6Plot
ggplot(df, aes(x = time, y = median)) +
geom_ribbon(aes(ymin = LCL, ymax = UCL), fill = "grey70", alpha = 0.5) +
geom_line() 
5.2 Morris screening
At the same exposure condition, conduct Morris elementary effects screening method to find the influential parameters for blood concentration during the 2-hr post exposure.
- Parameter setting
BW <- 0.043
Pb <- 42.7
Pl <- 1.96
Pf <- 36.4
Pm <- 0.609
Pvrg <- 1.96
Ppu <- 1.96
Pbr <- 1.96
VmaxC <- 6.39
VmaxClu <- 13.4
VmaxCvr <- 17.4
factors <- c("BW", "Pb", "Pl", "Pf", "Pm", "Pvrg", "Ppu", "Pbr", "VmaxC", "VmaxClu", "VmaxCvr")
baseline <- c(BW, Pb, Pl, Pf, Pm, Pvrg, Ppu, Pbr, VmaxC, VmaxClu, VmaxCvr)
limit <- c(1.2, rep(2, 7), rep(4, 3))
binf <- baseline/limit
bsup <- baseline*limit
binf## [1] 0.03583333 21.35000000 0.98000000 18.20000000 0.30450000
## [6] 0.98000000 0.98000000 0.98000000 1.59750000 3.35000000
## [11] 4.35000000bsup## [1] 0.0516 85.4000 3.9200 72.8000 1.2180 3.9200 3.9200 3.9200
## [9] 25.5600 53.6000 69.6000- Generate parameter matrix
set.seed(1234)
x <- morris(model = NULL, factors = factors, r = 512,
design = list(type = "oat", levels = 5, grid.jump = 3),
binf = binf, bsup = bsup)
X <- cbind(1, x$X)
write.table(X, file = "setpts.out", row.names = F, sep = "\t")
head(X)## BW Pb Pl Pf Pm Pvrg Ppu Pbr VmaxC VmaxClu
## [1,] 1 0.039775 37.3625 0.98 31.85 1.218 1.715 3.185 3.92 19.56937 3.3500
## [2,] 1 0.039775 37.3625 0.98 72.80 1.218 1.715 3.185 3.92 19.56937 3.3500
## [3,] 1 0.039775 85.4000 0.98 72.80 1.218 1.715 3.185 3.92 19.56937 3.3500
## [4,] 1 0.039775 85.4000 0.98 72.80 1.218 1.715 0.980 3.92 19.56937 3.3500
## [5,] 1 0.039775 85.4000 0.98 72.80 1.218 3.920 0.980 3.92 19.56937 3.3500
## [6,] 1 0.039775 85.4000 0.98 72.80 1.218 3.920 0.980 3.92 19.56937 41.0375
## VmaxCvr
## [1,] 4.35
## [2,] 4.35
## [3,] 4.35
## [4,] 4.35
## [5,] 4.35
## [6,] 4.35- Modeling
out <- mcsim(model = "EB.model.R", input = "EB_setpts.in.R", dir = "modeling/EB")## Execute: ./mcsim.EB.model.R.exe modeling/EB/EB_setpts.in.Rhead(out)## Iter BW Pb Pl Pf Pm Pvrg Ppu Pbr VmaxC VmaxClu
## 1 0 0.039775 37.3625 0.98 31.85 1.218 1.715 3.185 3.92 19.5694 3.3500
## 2 1 0.039775 37.3625 0.98 72.80 1.218 1.715 3.185 3.92 19.5694 3.3500
## 3 2 0.039775 85.4000 0.98 72.80 1.218 1.715 3.185 3.92 19.5694 3.3500
## 4 3 0.039775 85.4000 0.98 72.80 1.218 1.715 0.980 3.92 19.5694 3.3500
## 5 4 0.039775 85.4000 0.98 72.80 1.218 3.920 0.980 3.92 19.5694 3.3500
## 6 5 0.039775 85.4000 0.98 72.80 1.218 3.920 0.980 3.92 19.5694 41.0375
## VmaxCvr Cvtot_1.1 Cvtot_1.2 Cvtot_1.3 Cvtot_1.4 Cvtot_1.5
## 1 4.35 9.50553e-06 2.16779e-06 1.06171e-06 5.66928e-07 3.10056e-07
## 2 4.35 9.08168e-06 2.36341e-06 1.57338e-06 1.17212e-06 8.93469e-07
## 3 4.35 9.60993e-06 2.55974e-06 1.70937e-06 1.27599e-06 9.74858e-07
## 4 4.35 9.61063e-06 2.55525e-06 1.70731e-06 1.27467e-06 9.73874e-07
## 5 4.35 9.60570e-06 2.58665e-06 1.72075e-06 1.28325e-06 9.80242e-07
## 6 4.35 8.36660e-06 2.12817e-06 1.40701e-06 1.04477e-06 7.94037e-07- Visualization
str <- which(names(out) == "Cvtot_1.1")
end <- which(names(out) == "Cvtot_1.5")
par(mfrow = c(2,3), mar = c(2,2,2,1))
for (i in str:end){
tell(x, out[,i])
print(x)
plot(x,
main = paste("Time = ",4 + (i - str)*0.5, "hr"))
}##
## Call:
## morris(model = NULL, factors = factors, r = 512, design = list(type = "oat", levels = 5, grid.jump = 3), binf = binf, bsup = bsup)
##
## Model runs: 6144
## mu mu.star sigma
## BW -4.539138e-08 4.539190e-08 1.099540e-07
## Pb 9.393399e-07 9.393399e-07 2.066289e-06
## Pl -5.843099e-09 5.862786e-09 2.975201e-08
## Pf -4.930911e-07 4.930911e-07 1.189831e-06
## Pm -3.653966e-08 3.654039e-08 1.526393e-07
## Pvrg -4.954167e-09 4.963958e-09 2.572916e-08
## Ppu -1.114974e-09 1.133307e-09 5.132925e-09
## Pbr -5.919010e-10 6.068490e-10 3.723951e-09
## VmaxC -5.155248e-06 5.155248e-06 7.147079e-06
## VmaxClu -1.249138e-06 1.249138e-06 3.248143e-06
## VmaxCvr -8.187726e-06 8.187726e-06 8.061333e-06##
## Call:
## morris(model = NULL, factors = factors, r = 512, design = list(type = "oat", levels = 5, grid.jump = 3), binf = binf, bsup = bsup)
##
## Model runs: 6144
## mu mu.star sigma
## BW 5.946135e-08 6.134332e-08 8.811453e-08
## Pb 3.116289e-07 3.116289e-07 8.342762e-07
## Pl 1.904573e-08 1.907047e-08 7.443155e-08
## Pf 3.730996e-07 4.237563e-07 4.421595e-07
## Pm 4.324676e-07 4.324676e-07 5.534734e-07
## Pvrg 2.157932e-08 2.159228e-08 6.618055e-08
## Ppu 5.215727e-09 5.218378e-09 1.356543e-08
## Pbr 1.380672e-09 1.458161e-09 6.762572e-09
## VmaxC -1.297719e-06 1.297719e-06 2.431121e-06
## VmaxClu -4.100508e-07 4.100508e-07 1.259435e-06
## VmaxCvr -1.855217e-06 1.855217e-06 2.390441e-06##
## Call:
## morris(model = NULL, factors = factors, r = 512, design = list(type = "oat", levels = 5, grid.jump = 3), binf = binf, bsup = bsup)
##
## Model runs: 6144
## mu mu.star sigma
## BW 5.351556e-08 5.359379e-08 7.759246e-08
## Pb 1.811355e-07 1.811355e-07 4.592240e-07
## Pl 8.564219e-09 8.568033e-09 3.469396e-08
## Pf 6.574977e-07 6.603309e-07 5.892706e-07
## Pm 1.354224e-07 1.354224e-07 2.353699e-07
## Pvrg 9.607430e-09 9.609444e-09 2.791916e-08
## Ppu 2.412334e-09 2.414574e-09 6.423419e-09
## Pbr 6.327591e-10 6.580232e-10 3.214738e-09
## VmaxC -7.503134e-07 7.503134e-07 1.367775e-06
## VmaxClu -2.365413e-07 2.365413e-07 7.057592e-07
## VmaxCvr -1.081020e-06 1.081020e-06 1.453326e-06##
## Call:
## morris(model = NULL, factors = factors, r = 512, design = list(type = "oat", levels = 5, grid.jump = 3), binf = binf, bsup = bsup)
##
## Model runs: 6144
## mu mu.star sigma
## BW 4.666649e-08 4.666649e-08 5.838176e-08
## Pb 1.171406e-07 1.171406e-07 2.927991e-07
## Pl 4.625805e-09 4.630551e-09 1.734766e-08
## Pf 6.628573e-07 6.628573e-07 6.317125e-07
## Pm 6.473649e-08 6.473649e-08 1.174730e-07
## Pvrg 5.530547e-09 5.531018e-09 1.498474e-08
## Ppu 1.346258e-09 1.349947e-09 3.277857e-09
## Pbr 3.573281e-10 3.813115e-10 1.763603e-09
## VmaxC -4.832033e-07 4.832033e-07 9.048960e-07
## VmaxClu -1.503851e-07 1.503851e-07 4.506583e-07
## VmaxCvr -7.059492e-07 7.059492e-07 1.021888e-06##
## Call:
## morris(model = NULL, factors = factors, r = 512, design = list(type = "oat", levels = 5, grid.jump = 3), binf = binf, bsup = bsup)
##
## Model runs: 6144
## mu mu.star sigma
## BW 3.982414e-08 3.982414e-08 4.584417e-08
## Pb 8.042267e-08 8.042267e-08 2.069331e-07
## Pl 2.734992e-09 2.738771e-09 9.522463e-09
## Pf 5.750717e-07 5.750717e-07 5.738846e-07
## Pm 3.823143e-08 3.823143e-08 6.838545e-08
## Pvrg 3.502187e-09 3.504246e-09 9.243791e-09
## Ppu 8.244039e-10 8.244963e-10 1.833713e-09
## Pbr 2.184616e-10 2.313429e-10 1.029069e-09
## VmaxC -3.293028e-07 3.293028e-07 6.520604e-07
## VmaxClu -1.012795e-07 1.012795e-07 3.142924e-07
## VmaxCvr -4.867221e-07 4.867221e-07 7.606528e-07
Convergence check
# Create data frame
for (i in 3:11){
x <- morris(model = NULL, factors = factors, r = 2^i,
design = list(type = "oat", levels = 6, grid.jump = 3),
binf = binf, bsup = bsup)
dat <- cbind(1, x$X)
write.table(dat, file = "setpts.out", row.names = F, sep = "\t")
out <- mcsim(model = "EB.model.R", input = "EB_setpts.in.R", dir = "modeling/EB")
# Specific time point
y <- out$Cvtot_1.5
tell(x, y)
if (i == 3){ X <- apply(abs(x$ee), 2, mean) } else X <- rbind(X, apply(abs(x$ee), 2, mean))
}## Execute: ./mcsim.EB.model.R.exe modeling/EB/EB_setpts.in.R
## Execute: ./mcsim.EB.model.R.exe modeling/EB/EB_setpts.in.R
## Execute: ./mcsim.EB.model.R.exe modeling/EB/EB_setpts.in.R
## Execute: ./mcsim.EB.model.R.exe modeling/EB/EB_setpts.in.R
## Execute: ./mcsim.EB.model.R.exe modeling/EB/EB_setpts.in.R
## Execute: ./mcsim.EB.model.R.exe modeling/EB/EB_setpts.in.R
## Execute: ./mcsim.EB.model.R.exe modeling/EB/EB_setpts.in.R
## Execute: ./mcsim.EB.model.R.exe modeling/EB/EB_setpts.in.R
## Execute: ./mcsim.EB.model.R.exe modeling/EB/EB_setpts.in.RX## BW Pb Pl Pf Pm
## X 5.790646e-08 1.063699e-07 1.067521e-09 5.992944e-07 1.009211e-07
## 3.222659e-08 3.580549e-08 4.402062e-10 4.710151e-07 2.021413e-08
## 2.900282e-08 2.821310e-08 6.238573e-10 4.187986e-07 1.956760e-08
## 3.101394e-08 2.556600e-08 5.149719e-10 3.932252e-07 2.463100e-08
## 3.105886e-08 4.301121e-08 5.595169e-10 4.351427e-07 2.186930e-08
## 3.200012e-08 4.208949e-08 1.060956e-09 4.378947e-07 2.774631e-08
## 3.268151e-08 5.337593e-08 1.443239e-09 4.656815e-07 2.955017e-08
## 3.399517e-08 4.524042e-08 1.322251e-09 4.750906e-07 2.745396e-08
## 3.338046e-08 4.573180e-08 1.074968e-09 4.651997e-07 2.717093e-08
## Pvrg Ppu Pbr VmaxC VmaxClu
## X 1.277958e-08 7.252917e-10 7.350000e-11 2.305959e-07 2.808915e-07
## 2.697657e-09 3.541354e-10 3.041667e-11 3.033048e-07 8.732804e-08
## 1.080122e-09 4.603260e-10 1.005755e-10 2.357721e-07 4.600499e-08
## 1.089266e-09 4.337273e-10 5.998437e-11 1.532712e-07 3.673301e-08
## 2.477879e-09 4.114302e-10 1.320648e-10 1.760887e-07 6.387713e-08
## 1.898244e-09 5.271701e-10 1.239016e-10 2.159678e-07 4.604850e-08
## 2.543221e-09 4.968879e-10 1.036403e-10 2.247941e-07 4.922868e-08
## 2.247678e-09 5.572310e-10 1.068919e-10 2.269863e-07 6.268561e-08
## 2.202031e-09 5.558069e-10 1.103888e-10 2.192794e-07 5.743721e-08
## VmaxCvr
## X 4.702080e-07
## 3.476158e-07
## 2.843540e-07
## 2.921703e-07
## 3.619865e-07
## 3.548697e-07
## 3.702963e-07
## 3.964481e-07
## 3.708041e-07# Plot
ylim <- range(X)
plot(X[,1], ylim = ylim, type = "b", xaxt = "n")
axis(1, at = seq(3,11,2), labels = 2^seq(5,13,2))
for(i in 2 :11){
lines(X[,i], type = "b")
}
5.3 Fourier amplitude sensitivity testing
Use pksensi and conduct FAST method to find the non-influential parameter for blood concentrations.
params <- c("BW", "Pb", "Pl", "Pf", "Pm", "Pvrg", "Ppu", "Pbr", "VmaxC", "VmaxClu", "VmaxCvr")
q <- "qunif"
q.arg <- list(list(min = 0.0358, max = 0.0516),
list(min = 21.35, max = 85.4),
list(min = 0.98, max = 3.92),
list(min = 18.2, max = 72.8),
list(min = 0.3045, max = 1.218),
list(min = 0.98, max = 3.92),
list(min = 0.98, max = 3.92),
list(min = 0.98, max = 3.92),
list(min = 1.5975, max = 25.56),
list(min = 3.35, max = 53.6),
list(min = 4.35, max = 69.6)) - Generate parameter space
set.seed(1234)
x <- rfast99(params, n = 8000, q = q, q.arg = q.arg, replicate = 1) # convergence sizedim(x$a) # c(Evaluation, replication, parameters)## [1] 88000 1 11- Set experiment time-points, output variables, and conditions
var <- c("Cvtot")
times <- c(4.0, 4.5, 5, 5.5, 6)
condition <- c("Cppm = NDoses (2, 100, 0, 0, 4 )")- Modeling
Note: There is a problem to define rtol = 1e-11 and atol = 1e-13. Update pksensi to the latest version by devtools::install_github("nanhung/pksensi")
fast <- solve_mcsim(x, mName = "EB.model.R", params = params, time = times, vars = var,
condition = condition, rtol = 1e-11,atol = 1e-13)## Starting time: 2019-05-29 12:03:25## * Created input file "sim.in".## Execute: ./mcsim.EB.model.R.exe sim.in## Ending time: 2019-05-29 12:03:59- Visualization and decision
plot(fast)
pksim(fast)
check(fast, SI.cutoff = 0.05)##
## Sensitivity check ( Index > 0.05 )
## ----------------------------------
## First order:
## Pf VmaxC VmaxCvr
##
## Interaction:
## BW Pb Pl Pf Pm Pvrg Ppu Pbr VmaxC VmaxClu VmaxCvr
##
## Total order:
## BW Pb Pl Pf Pm Pvrg Ppu Pbr VmaxC VmaxClu VmaxCvr
##
## Unselected factors in total order:
##
##
##
## Convergence check ( Index > 0.05 )
## ----------------------------------
## First order:
##
##
## Interaction:
##
##
## Total order:
## heat_check(fast, SI.cutoff = c(0.05, 0.09, 0.15))
5.4 Morris vs. FAST
Compare the sensitivity measures (first order, interaction, total order) for Morris (exercise 2) and FAST (exercise 4).
set.seed(1234)
morrx <- morris(model = NULL, factors = factors, r = 512,
design = list(type = "oat", levels = 5, grid.jump = 3),
binf = binf, bsup = bsup)
X <- cbind(1, morrx$X)
write.table(X, file = "setpts.out", row.names = F, sep = "\t")
out <- mcsim(model = "EB.model.R", input = "EB_setpts.in.R", dir = "modeling/EB")## Execute: ./mcsim.EB.model.R.exe modeling/EB/EB_setpts.in.Ry <- "Cvtot_1.2"
str <- which(names(out) == y)
tell(morrx, out[,str])
mu_star <- apply(abs(morrx$ee), 2, mean)
mu_star## BW Pb Pl Pf Pm
## 6.134332e-08 3.116289e-07 1.907047e-08 4.237563e-07 4.324676e-07
## Pvrg Ppu Pbr VmaxC VmaxClu
## 2.159228e-08 5.218378e-09 1.458161e-09 1.297719e-06 4.100508e-07
## VmaxCvr
## 1.855217e-06fast$mSI## , , Cvtot
##
## BW Pb Pl Pf Pm
## 4 6.115264e-06 0.002355767 2.800047e-06 0.000675948 1.541967e-06
## 4.5 8.600197e-04 0.003288462 7.005449e-06 0.058171149 2.962126e-02
## 5 1.650442e-03 0.002388506 5.650310e-06 0.243576101 4.386352e-03
## 5.5 2.531170e-03 0.001816455 4.734416e-06 0.384911155 1.670921e-03
## 6 3.296907e-03 0.001476888 4.175366e-06 0.469768580 1.046634e-03
## Pvrg Ppu Pbr VmaxC VmaxClu
## 4 2.460485e-06 5.261785e-07 2.108597e-06 0.16761944 0.002916109
## 4.5 1.876890e-05 3.148733e-06 7.108806e-06 0.12862622 0.004714960
## 5 1.361004e-05 3.469226e-06 7.547569e-06 0.10110815 0.003870680
## 5.5 1.300877e-05 4.238494e-06 8.318794e-06 0.07682976 0.003033988
## 6 1.312728e-05 5.035698e-06 9.471238e-06 0.06106666 0.002464525
## VmaxCvr
## 4 0.7057013
## 4.5 0.5995002
## 5 0.4757105
## 5.5 0.3645130
## 6 0.2919153dim(fast$mSI)## [1] 5 11 1plot(mu_star, fast$mSI[2,,1], ylim=c(0,1))
abline(0.05, 0)
5.5 Acetaminophen-PBPK model
Reproduce the published result of acetaminophen-PBPK model by following the vignette in pksensi.
https://nanhung.rbind.io/pksensi/articles/pbpk_apap.html
- Create & compile model file
pksensi::pbpk_apap_model()
makemcsim(model = "pbpk_apap.model")## * The 'pbpk_apap.model' has been moved to modeling folder.## * Created executable program 'mcsim.pbpk_apap.model.exe'.- Define arguments
Tg <- log(0.23)
Tp <- log(0.033)
CYP_Km <- log(130)
SULT_Km_apap <- log(300)
SULT_Ki <- log(526)
SULT_Km_paps <- log(0.5)
UGT_Km <- log(6.0e3)
UGT_Ki <- log(5.8e4)
UGT_Km_GA <-log(0.5)
Km_AG <- log(1.99e4)
Km_AS <- log(2.29e4)
rng <- 1.96
# Parameter distribution
params <- c("lnTg", "lnTp", "lnCYP_Km","lnCYP_VmaxC",
"lnSULT_Km_apap","lnSULT_Ki","lnSULT_Km_paps","lnSULT_VmaxC",
"lnUGT_Km","lnUGT_Ki","lnUGT_Km_GA","lnUGT_VmaxC",
"lnKm_AG","lnVmax_AG","lnKm_AS","lnVmax_AS",
"lnkGA_syn","lnkPAPS_syn", "lnCLC_APAP","lnCLC_AG","lnCLC_AS")
q <- "qunif"
q.arg <-list(list(Tg-rng, Tg+rng),
list(Tp-rng, Tp+rng),
list(CYP_Km-rng, CYP_Km+rng),
list(-2., 5.),
list(SULT_Km_apap-rng, SULT_Km_apap+rng),
list(SULT_Ki-rng, SULT_Ki+rng),
list(SULT_Km_paps-rng, SULT_Km_paps+rng),
list(0, 10),
list(UGT_Km-rng, UGT_Km+rng),
list(UGT_Ki-rng, UGT_Ki+rng),
list(UGT_Km_GA-rng, UGT_Km_GA+rng),
list(0, 10),
list(Km_AG-rng, Km_AG+rng),
list(7., 15),
list(Km_AS-rng, Km_AS+rng),
list(7., 15),
list(0., 13),
list(0., 13),
list(-6., 1),
list(-6., 1),
list(-6., 1))
# Mutivariate
vars <- c("lnCPL_APAP_mcgL", "lnCPL_AG_mcgL", "lnCPL_AS_mcgL")
times <- seq(from = 0.1, to = 12.1, by = 0.2)
# Conditions
conditions <- c("mgkg_flag = 1",
"OralExp_APAP = NDoses(2, 1, 0, 0, 0.001)",
"OralDose_APAP_mgkg = 20.0")- Generate Parameter space
set.seed(1234)
x <- rfast99(params = params, n = 1000, q = q, q.arg = q.arg, replicate = 1) - Modeling
system.time(out <- solve_mcsim(x, mName = "pbpk_apap.model",
params = params,
vars = vars,
time = times,
condition = conditions))## Starting time: 2019-05-29 12:04:11## * Created input file "sim.in".## Execute: ./mcsim.pbpk_apap.model.exe sim.in## Ending time: 2019-05-29 12:05:00## user system elapsed
## 48.606 0.418 48.964- Visualization & Decision
plot(out, vars = "lnCPL_AG_mcgL")
heat_check(out)
check(out)##
## Sensitivity check ( Index > 0.05 )
## ----------------------------------
## First order:
## lnTg lnTp lnSULT_VmaxC lnUGT_Km lnUGT_VmaxC lnKm_AG lnVmax_AG lnKm_AS lnVmax_AS lnkGA_syn lnkPAPS_syn lnCLC_APAP lnCLC_AG lnCLC_AS
##
## Interaction:
## lnTg lnSULT_Km_apap lnSULT_VmaxC lnUGT_VmaxC lnVmax_AG lnVmax_AS lnkGA_syn lnkPAPS_syn lnCLC_APAP lnCLC_AG lnCLC_AS
##
## Total order:
## lnTg lnTp lnSULT_Km_apap lnSULT_VmaxC lnUGT_Km lnUGT_VmaxC lnKm_AG lnVmax_AG lnKm_AS lnVmax_AS lnkGA_syn lnkPAPS_syn lnCLC_APAP lnCLC_AG lnCLC_AS
##
## Unselected factors in total order:
## lnCYP_Km lnCYP_VmaxC lnSULT_Ki lnSULT_Km_paps lnUGT_Ki lnUGT_Km_GA
##
##
## Convergence check ( Index > 0.05 )
## ----------------------------------
## First order:
##
##
## Interaction:
##
##
## Total order:
##