0 Prerequisites
0 Basic MOnte Carlo practice
Task 1: Uncertainty analysis & Model calibration for simple PK model
Reference

0 Prerequisites

Import the R packages

library(simuloR)
library(tidyverse)
library(pksensi)
library(httk)
theme_set(theme_bw())

0 Basic MOnte Carlo practice

x <- rnorm(n = 28, mean = 10, sd =  2)
x

##  [1] 10.282183 10.052703  9.019463  8.889859 10.924005 10.789138  7.760453
##  [8] 10.136920  6.655561 10.464955 11.137035  7.568546  7.999040 11.867852
## [15]  9.291173  8.983861 10.954530 10.556185 10.505030 12.947875  8.217614
## [22]  5.748279 12.011019  9.966287  8.328393  9.461880  6.958836 11.449811

mean(x)

## [1] 9.604589

sd(x)

## [1] 1.7421

hist(x)

Generating random number

runif(n=5, min=0, max=1)

## [1] 0.8836757 0.4342189 0.4667599 0.6882098 0.3669369

runif(n=5, min=0, max=1)

## [1] 0.6830976 0.1962670 0.5833900 0.3285861 0.9533697

set.seed(1)
runif(n=5, min=0, max=1)

## [1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819

set.seed(1)
runif(n=5, min=0, max=1)

## [1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819

x <- runif(10)
plot(x)

hist(x)

x <- runif(100)
plot(x)

hist(x)

x <- runif(1000)
plot(x)

hist(x)

x <- runif(10000)
plot(x)

hist(x)

Task 1: Uncertainty analysis & Model calibration for simple PK model

In the previous exercise, we find that the predcited result can not used to describe the real cases. Therefore, we need to conduct the uncertainty analysis to figure out how to reset the model parameter.
Use the parameter distributions that we test in uncertainty analysis and conduct MCMC simulation to do model calibration.

`MonteCarlo()` specification

MonteCarlo specification gives general information required for the runs: the output file name, the number of runs to perform, and a starting seed for the random number generator. Its syntax is:

MonteCarlo("<OutputFilename>", <nRuns>, <RandomSeed>);

The character string , enclosed in double quotes, should be a valid filename for your operating system. If a null-string "" is given, the default name simmc.out will be used. The number of runs should be an integer, and is only limited by either your storage space for the output file or the largest (long) integer available on your machine. When using the GNU Scientific Library, the seed of the pseudo-random number generator can be any positive integer (including zero). When using the native code, it can be any positive real number (seeds between 1.0 and 2147483646.0 are used as is, the others are rescaled silently within those bounds). Here is an example of use:

MonteCarlo("percsimmc.out", 50000, 9386);

`Distrib()` specification

Distrib() specifies the name of the parameter to sample, and its sampling distribution. Its syntax is:

Distrib(<parameter identifier>, <distribution-name>, [<shape parameters>]);

For example Normal, takes two reals numbers as parameters: the mean and the standard deviation, the latter being strictly positive. LogNormal, takes two reals numbers as parameters: the geometric mean (exponential of the mean in log-space) and the geometric standard deviation (exponential, strictly superior to 1, of the standard deviation in log-space).

parms <- httk::parameterize_1comp(chem.name = "theophylline")

## Warning in available_rblood2plasma(chem.cas = chem.cas, chem.name =
## chem.name, : Human in vivo Rblood2plasma returned.

Fgutabs <- parms$Fgutabs * parms$hepatic.bioavailability
kgutabs <- parms$kgutabs
kelim <- parms$kelim
Vdist <- parms$Vdist

LL <- 0.5
Fgutabs_LL <- 0.8
kgutabs_LL <- kgutabs * LL
kelim_LL <- kelim * LL
Vdist_LL <- Vdist * LL
c(Fgutabs_LL, kgutabs_LL, kelim_LL, Vdist_LL)

## [1] 0.8000000 1.0900000 0.1326966 0.3900998

UL <- 2
Fgutabs_UL <- 1
kgutabs_UL <- kgutabs * UL
kelim_UL <- kelim * UL
Vdist_UL <- Vdist * UL
c(Fgutabs_UL, kgutabs_UL, kelim_UL, Vdist_UL)

## [1] 1.0000000 4.3600000 0.5307864 1.5603994

set.seed(1)
dist_Fgutabs <- runif(1000, Fgutabs_LL, Fgutabs_UL)
set.seed(1)
dist_kgutabs <- runif(1000, kgutabs_LL, kgutabs_UL)
set.seed(1)
dist_kelim <- runif(1000, kelim_LL, kelim_UL)
set.seed(1)
dist_Vdist <- runif(1000, Vdist_LL, Vdist_UL)

df <- data.frame(dist_Fgutabs, dist_kgutabs, dist_kelim, dist_Vdist)
head(df)

##   dist_Fgutabs dist_kgutabs dist_kelim dist_Vdist
## 1    0.8531017     1.958213  0.2383929  0.7008245
## 2    0.8744248     2.306845  0.2808353  0.8255963
## 3    0.9145707     2.963230  0.3607437  1.0605099
## 4    0.9816416     4.059839  0.4942449  1.4529750
## 5    0.8403364     1.749500  0.2129841  0.6261281
## 6    0.9796779     4.027734  0.4903364  1.4414849

names(df) <- c("Fgutabs", "kgutabs", "kelim", "Vdist")
names(df)

## [1] "Fgutabs" "kgutabs" "kelim"   "Vdist"

par(mfrow=c(1,4))
hist(df$Fgutabs)
hist(df$kgutabs)
hist(df$kelim)
hist(df$Vdist)

ffpk <- function(D, Fa, ka, ke, V, t){
    MW <- 180.167
    A <- (D * Fa * ka)/(V * (ka - ke))
    C <- A * exp(-ke * t) - A * exp(-ka * t)
    C_uM <- C/MW*1000
    return(C_uM)
}

S1 <- subset(Theoph, Subject==1)
t1 <- S1$Time
ffpk(D = 4.02, Fa = df$Fgutabs[1],ka = df$kgutabs[1], 
     ke = df$kelim[1], V = df$Vdist[1], t = t1)

##  [1]  0.00000000 10.18206001 16.86735666 20.22891128 18.51444367
##  [6] 12.42266743  9.16716362  5.78737712  3.57559277  1.71989793
## [11]  0.09273141

m.Conc <- matrix(nrow=1000, ncol=11)
colnames(m.Conc) <- t1

m.Conc[1,] <- ffpk(D = 4.02, Fa = df$Fgutabs[1],ka = df$kgutabs[1], 
                   ke = df$kelim[1], V = df$Vdist[1], t = t1)
head(m.Conc)

##       0     0.25     0.57     1.12     2.02     3.82      5.1     7.03
## [1,]  0 10.18206 16.86736 20.22891 18.51444 12.42267 9.167164 5.787377
## [2,] NA       NA       NA       NA       NA       NA       NA       NA
## [3,] NA       NA       NA       NA       NA       NA       NA       NA
## [4,] NA       NA       NA       NA       NA       NA       NA       NA
## [5,] NA       NA       NA       NA       NA       NA       NA       NA
## [6,] NA       NA       NA       NA       NA       NA       NA       NA
##          9.05    12.12      24.37
## [1,] 3.575593 1.719898 0.09273141
## [2,]       NA       NA         NA
## [3,]       NA       NA         NA
## [4,]       NA       NA         NA
## [5,]       NA       NA         NA
## [6,]       NA       NA         NA

for(i in 2:1000){
  m.Conc[i,] <- ffpk(D = 4.02, Fa = df$Fgutabs[i],ka = df$kgutabs[i], 
                     ke = df$kelim[i], V = df$Vdist[i], t = t1)
}
tail(m.Conc)

##         0      0.25     0.57     1.12      2.02      3.82       5.1
##  [995,] 0 10.449687 18.46668 24.24117 24.612671 18.746916 14.845385
##  [996,] 0  9.189284 12.16230 11.06759  7.638540  3.453179  1.961303
##  [997,] 0  9.454243 13.25610 12.88258  9.490076  4.755451  2.902034
##  [998,] 0 10.148504 16.67826 19.78739 17.896759 11.833968  8.657326
##  [999,] 0 10.332812 17.74770 22.37477 21.665128 15.572100 11.951813
## [1000,] 0 10.182003 16.86704 20.22815 18.513377 12.421643  9.166273
##               7.03      9.05      12.12       24.37
##  [995,] 10.3758745 7.1250664 4.02398214 0.411661584
##  [996,]  0.8358155 0.3422962 0.08813778 0.000392609
##  [997,]  1.3780912 0.6320787 0.19333381 0.001711913
##  [998,]  5.3959568 3.2894790 1.55044486 0.077084163
##  [999,]  7.9906066 5.2406410 2.76025043 0.213758611
## [1000,]  5.7866897 3.5750874 1.71959587 0.092702436

plot(t1, m.Conc[1,], type="l")

plot(t1, m.Conc[1,], type="l", col="grey")
for(i in 2:1000){
  lines(t1, m.Conc[i,], col="grey")
}

max.conc <- max(m.Conc)
plot(t1, m.Conc[1,], type="l", col="grey", ylim = c(0, max.conc))
for(i in 2:1000){
  lines(t1, m.Conc[i,], col="grey")
}

max.conc <- max(m.Conc)
plot(t1, m.Conc[1,], type="l", col="grey", ylim = c(0, max.conc))
for(i in 2:1000){
  lines(t1, m.Conc[i,], col="grey")
}
Conc1 <- S1$conc
points(t1, Conc1)

MCSim

Input file

MonteCarlo ("simmc.out", 
            1000, 
            10101010);

MW       = 180.164; # Molecular weight (g/mole)
IngDose  = 4.02;    # ingested dose (mg)

Distrib (Fgutabs, Uniform, 0.8, 1);
Distrib (kgutabs, Uniform, 1.09, 4.36);
Distrib (kelim, LogNormal, 0.2653932, 2); 
Distrib (Vdist, Uniform, 0.39, 1.56); 

Simulation {
  Print (Ccompartment 0.00   0.25   0.57   1.12   2.02   3.82   5.10   7.03   9.05  12.12  24.37);
}
End.

set.seed(2222)
ffpk.mtc.out <- mcsim("models/ffpk.model", "inputs/ffpk.mtc.in")

## * Creating executable program, pleas wait...

## * Created executable program 'models/mcsim.ffpk.model'.

## Executing...

head(ffpk.mtc.out)

##   Iter  Fgutabs kgutabs     kelim    Vdist Conc_1.1 Conc_1.2 Conc_1.3
## 1    0 0.807641 1.48323 0.0814536 0.918445  1.08347  1.96472  2.70396
## 2    1 0.817178 2.98160 0.5490830 0.820539  1.94910  2.69159  2.47912
## 3    2 0.935946 1.80647 0.1324660 1.138770  1.17955  2.03286  2.60240
## 4    3 0.816768 1.27744 0.2548760 0.837924  1.03609  1.86984  2.50895
## 5    4 0.881173 2.31719 0.0479567 1.525620  1.01427  1.67416  2.07003
## 6    5 0.821408 4.33600 0.1099600 1.124480  1.91213  2.57541  2.64036
##   Conc_1.4 Conc_1.5 Conc_1.6 Conc_1.7  Conc_1.8   Conc_1.9   Conc_1.10
## 1  2.98601 2.727290 2.466990 2.109660 1.7896900 1.39373000 5.13849e-01
## 2  1.60671 0.602373 0.298312 0.103380 0.0340989 0.00631912 7.57662e-06
## 3  2.63562 2.145980 1.813970 1.405010 1.0751600 0.71590900 1.41294e-01
## 4  2.55455 1.811770 1.327020 0.815234 0.4874980 0.22294100 9.82254e-03
## 5  2.13007 1.973730 1.856520 1.692420 1.5361600 1.32586000 7.36820e-01
## 6  2.41233 1.979540 1.719640 1.390820 1.1138000 0.79469200 2.06634e-01

par(mfrow = c(1,4))
hist(ffpk.mtc.out$Fgutabs, main = "Fgutabs", xlab = "", ylab = "")
hist(ffpk.mtc.out$kgutabs, main = "kgutabs", xlab = "", ylab = "")
hist(ffpk.mtc.out$kelim, main = "kelim", xlab = "", ylab = "")
hist(ffpk.mtc.out$Vdist, main = "Vdist", xlab = "", ylab = "")

Conc.str <- which(names(ffpk.mtc.out)=="Conc_1.1")
Conc.end <- which(names(ffpk.mtc.out)=="Conc_1.10")
df.Conc <- ffpk.mtc.out[,Conc.str:Conc.end]
head(df.Conc)

##   Conc_1.1 Conc_1.2 Conc_1.3 Conc_1.4 Conc_1.5 Conc_1.6 Conc_1.7  Conc_1.8
## 1  1.08347  1.96472  2.70396  2.98601 2.727290 2.466990 2.109660 1.7896900
## 2  1.94910  2.69159  2.47912  1.60671 0.602373 0.298312 0.103380 0.0340989
## 3  1.17955  2.03286  2.60240  2.63562 2.145980 1.813970 1.405010 1.0751600
## 4  1.03609  1.86984  2.50895  2.55455 1.811770 1.327020 0.815234 0.4874980
## 5  1.01427  1.67416  2.07003  2.13007 1.973730 1.856520 1.692420 1.5361600
## 6  1.91213  2.57541  2.64036  2.41233 1.979540 1.719640 1.390820 1.1138000
##     Conc_1.9   Conc_1.10
## 1 1.39373000 5.13849e-01
## 2 0.00631912 7.57662e-06
## 3 0.71590900 1.41294e-01
## 4 0.22294100 9.82254e-03
## 5 1.32586000 7.36820e-01
## 6 0.79469200 2.06634e-01

max.Conc <- max(df.Conc)
plot(t1[-1], df.Conc[1,], type="l", col="grey", ylim=c(0,max.Conc))
for(i in 2:1000){
  lines(t1[-1], df.Conc[i,], col="grey")
}
points(t1, Conc1)

pksensi

Set parameter distributions (assign parms, dist, q.qarg)

parameters <- c("Fgutabs",  "kgutabs", "kelim", "Vdist")  
dist <- c("Uniform", "Uniform", "LogNormal", "Uniform") # MCSim definition
q.arg <- list(list(0.8, 1.0),
            list(1.09, 4.36),
            list(0.2653932, 2),
            list(0.39, 1.56))

Set experiment time-points, output variables, and conditions

outputs <- c("Conc")
times <- c(0.00, 0.25, 0.57, 1.12, 2.02, 3.82, 5.10, 7.03, 9.05, 12.12, 24.37)
conditions <- c("MW = 180.164", "Period = 48", "IngDose = 4.02")

Modeling

set.seed(2222)
y<-solve_mcsim(mName = "models/pbtk1cpt.model", params = parameters, vars = outputs, monte_carlo = 1000,
               dist = dist, q.arg = q.arg, time = times, condition = conditions)

## Starting time: 2019-12-06 10:22:46

## * Created input file "sim.in".

## Execute: models/mcsim.pbtk1cpt.model sim.in

## Ending time: 2019-12-06 10:22:47

pksim(y)
points(t1, Conc1)

Calibration

set.seed(1)
ffpk.mcmc.out <- mcsim(model = "models/ffpk.model", input = "inputs/ffpk.mcmc.in")

## * Creating executable program, pleas wait...

## * Created executable program 'models/mcsim.ffpk.model'.

## Executing...

## * Create 'MCMC.check.out' from the last iteration.

par(mfrow = c(2,2))
plot(ffpk.mcmc.out$Fgutabs.1., type = "l")
plot(ffpk.mcmc.out$kgutabs.1., type = "l")
plot(ffpk.mcmc.out$kelim.1., type = "l")
plot(ffpk.mcmc.out$Vdist.1., type = "l")

check.out <- read.delim("MCMC.check.out") 
plot(check.out$Data, check.out$Prediction, xlim=c(0,10), ylim=c(0,10))
abline(0,1)

Use R

IngDose <- 4.02
Fgutabs <- ffpk.mcmc.out$Fgutabs.1.
kgutabs <- ffpk.mcmc.out$kgutabs.1.
kelim <- ffpk.mcmc.out$kelim.1. 
Vdist <- ffpk.mcmc.out$Vdist.1.
MW <- 180.164

t <- Theoph$Time[1:11]
obs<- Theoph$con[1:11]
A <- (IngDose * Fgutabs * kgutabs)/(Vdist * (kgutabs - kelim))

Conc.matrix <- matrix(nrow = length(Fgutabs), ncol = length(t))
for(i in 1:11){
  Conc.matrix[,i] <- A * exp(-kelim * t[i]) - A * exp(-kgutabs * t[i])
}
plot(t, obs, pch = "")
for(i in 501:1000){
  lines(t, Conc.matrix[i,], col="grey")
}
points(t, obs)

Use setpoints

ffpk.setpts.out <- mcsim("models/ffpk.model", input = "inputs/ffpk.setpts.in")

## * Creating executable program, pleas wait...

## * Created executable program 'models/mcsim.ffpk.model'.

## Executing...

names(ffpk.setpts.out)

##  [1] "Iter"      "Fgutabs"   "kgutabs"   "kelim"     "Vdist"    
##  [6] "Conc_1.1"  "Conc_1.2"  "Conc_1.3"  "Conc_1.4"  "Conc_1.5" 
## [11] "Conc_1.6"  "Conc_1.7"  "Conc_1.8"  "Conc_1.9"  "Conc_1.10"
## [16] "Conc_1.11"

t <- Theoph$Time[1:11]
obs<- Theoph$conc[1:11]
str.pt <- which(names(ffpk.setpts.out)=="Conc_1.1") 
end.pt <- which(names(ffpk.setpts.out)=="Conc_1.11") 
par(mfrow = c(1,1))
plot(t, obs, pch = "")
for(i in 501:1000){
  lines(t, ffpk.setpts.out[i,str.pt:end.pt], col="grey")
}
points(t, obs)

Use setpoints - Global evaluation

sessionInfo()

## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: elementary OS 5.0 Juno
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] httk_1.10.1        pksensi_1.1.4.9000 forcats_0.4.0     
##  [4] stringr_1.4.0      dplyr_0.8.3        purrr_0.3.3       
##  [7] readr_1.3.1        tidyr_1.0.0        tibble_2.1.3      
## [10] ggplot2_3.2.1      tidyverse_1.3.0    simuloR_0.0.1     
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.3        lubridate_1.7.4   msm_1.6.7        
##  [4] mvtnorm_1.0-11    lattice_0.20-38   getPass_0.2-2    
##  [7] assertthat_0.2.1  zeallot_0.1.0     digest_0.6.23    
## [10] R6_2.4.1          cellranger_1.1.0  plyr_1.8.4       
## [13] backports_1.1.5   reprex_0.3.0      survey_3.36      
## [16] evaluate_0.14     httr_1.4.1        pillar_1.4.2     
## [19] rlang_0.4.2       lazyeval_0.2.2    readxl_1.3.1     
## [22] rstudioapi_0.10   data.table_1.12.6 Matrix_1.2-17    
## [25] rmarkdown_1.15    splines_3.6.1     munsell_0.5.0    
## [28] broom_0.5.2       compiler_3.6.1    modelr_0.1.5     
## [31] xfun_0.9          pkgconfig_2.0.3   mitools_2.4      
## [34] htmltools_0.3.6   tidyselect_0.2.5  expm_0.999-4     
## [37] reshape_0.8.8     crayon_1.3.4      dbplyr_1.4.2     
## [40] withr_2.1.2       grid_3.6.1        nlme_3.1-142     
## [43] jsonlite_1.6      gtable_0.3.0      lifecycle_0.1.0  
## [46] DBI_1.0.0         magrittr_1.5      scales_1.1.0     
## [49] cli_1.1.0         stringi_1.4.3     fs_1.3.1         
## [52] xml2_1.2.2        generics_0.0.2    vctrs_0.2.0      
## [55] deSolve_1.25      tools_3.6.1       glue_1.3.1       
## [58] hms_0.5.2         survival_2.44-1.1 yaml_2.2.0       
## [61] colorspace_1.4-1  rvest_0.3.5       knitr_1.25       
## [64] haven_2.2.0

Reference

Bois F.Y. 2013. Bayesian Inference. In: Reisfeld B., Mayeno A. (eds) Computational Toxicology. Methods in Molecular Biology (Methods and Protocols), vol 930. Humana Press, Totowa, NJ

Advancing applications (physiologically based pharmacokinetic modeling)

Nan-Hung Hsieh

2019/12/09 (update: 2019-12-06)

0 Prerequisites

0 Basic MOnte Carlo practice

Task 1: Uncertainty analysis & Model calibration for simple PK model

`MonteCarlo()` specification

`Distrib()` specification

Calibration

Reference

Advancing applications (physiologically based pharmacokinetic modeling)

Nan-Hung Hsieh

2019/12/09 (update: 2019-12-06)

0 Prerequisites

0 Basic MOnte Carlo practice

Task 1: Uncertainty analysis & Model calibration for simple PK model

MonteCarlo() specification

Distrib() specification

Calibration

Reference

`MonteCarlo()` specification

`Distrib()` specification