Exercise 2.1: Uncertainty analysis on PK model
Nan-Hung Hsieh
2019/12/09 (update: 2019-12-07)
TASK
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.
Basic MOnte Carlo practice
x <- rnorm(n = 28, mean = 10, sd = 2)
x## [1] 11.590906 10.836967 12.615093 10.180441 13.498127 8.633502 9.741254
## [8] 9.544412 9.893518 7.767800 11.445249 8.724177 12.450644 10.911855
## [15] 8.236261 7.768533 10.566323 6.872776 10.941408 7.743053 9.316467
## [22] 8.405518 8.366940 7.035457 9.441074 6.118694 6.984528 11.402068mean(x)## [1] 9.536895sd(x)## [1] 1.915777hist(x)
Generating random number
runif(n=5, min=0, max=1)## [1] 0.6992383 0.7317137 0.9417769 0.9002684 0.1346130runif(n=5, min=0, max=1)## [1] 0.6661019 0.1851698 0.5794496 0.5384452 0.2722011set.seed(1)
runif(n=5, min=0, max=1)## [1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819set.seed(1)
runif(n=5, min=0, max=1)## [1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819x <- 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)
Monte Carlo method in GNU MCSim
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
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).
SOLUTION
Packages
library(simuloR)
library(tidyverse)
library(pksensi)
library(httk)Uncertainty analysis of PK model in R
Model
ffpk <- function(IngDose, Fgutabs, kgutabs, kelim, Vdist, BW, t){
A <- (IngDose * Fgutabs * kgutabs)/(Vdist * BW * (kgutabs - kelim))
Conc <- A * exp(-kelim * t) - A * exp(-kgutabs * t)
return(Conc)
}Parameters
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$VdistMonte Carlo simulation
LL <- 0.2
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.80000000 0.43600000 0.05307864 0.15603994UL <- 5
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.000000 10.900000 1.326966 3.900998n <- 1000
set.seed(1)
dist_Fgutabs <- runif(n, Fgutabs_LL, Fgutabs_UL)
set.seed(1)
dist_kgutabs <- runif(n, kgutabs_LL, kgutabs_UL)
set.seed(1)
dist_kelim <- runif(n, kelim_LL, kelim_UL)
set.seed(1)
dist_Vdist <- runif(n, 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 3.214283 0.3913068 1.1503589
## 2 0.8744248 4.329904 0.5271226 1.5496285
## 3 0.9145707 6.430338 0.7828293 2.3013520
## 4 0.9816416 9.939486 1.2100331 3.5572404
## 5 0.8403364 2.546400 0.3099987 0.9113304
## 6 0.9796779 9.836750 1.1975259 3.5204720names(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)
S1 <- subset(Theoph, Subject==1)
t1 <- S1$Time
Dose1 <- mean(S1$Dose)
BW1 <- mean(S1$Wt)m.Conc <- matrix(nrow=n, ncol=11)
head(m.Conc)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] NA NA NA NA NA NA NA NA NA NA NA
## [2,] NA NA NA NA NA NA NA NA NA NA NA
## [3,] NA NA NA NA NA NA NA NA NA NA NA
## [4,] NA NA NA NA NA NA NA NA NA NA NA
## [5,] NA NA NA NA NA NA NA NA NA NA NA
## [6,] NA NA NA NA NA NA NA NA NA NA NAcolnames(m.Conc) <- t1
head(m.Conc)## 0 0.25 0.57 1.12 2.02 3.82 5.1 7.03 9.05 12.12 24.37
## [1,] NA NA NA NA NA NA NA NA NA NA NA
## [2,] NA NA NA NA NA NA NA NA NA NA NA
## [3,] NA NA NA NA NA NA NA NA NA NA NA
## [4,] NA NA NA NA NA NA NA NA NA NA NA
## [5,] NA NA NA NA NA NA NA NA NA NA NA
## [6,] NA NA NA NA NA NA NA NA NA NA NAdim(m.Conc)## [1] 1000 11m.Conc[1,] <- ffpk(IngDose = Dose1*BW1,
Fgutabs = df$Fgutabs[1], kgutabs = df$kgutabs[1],
kelim = df$kelim[1],
Vdist = df$Vdist[1],
BW = BW1,
t = t1)
head(m.Conc)## 0 0.25 0.57 1.12 2.02 3.82 5.1 7.03
## [1,] 0 1.558323 2.172482 2.097204 1.53474 0.7613451 0.4613856 0.2168081
## [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,] 0.09835398 0.02958481 0.0002450613
## [2,] NA NA NA
## [3,] NA NA NA
## [4,] NA NA NA
## [5,] NA NA NA
## [6,] NA NA NAfor(i in 2:n){
m.Conc[i,] <- ffpk(IngDose = Dose1*BW1,
Fgutabs = df$Fgutabs[i], kgutabs = df$kgutabs[i],
kelim = df$kelim[i],
Vdist = df$Vdist[i],
BW = BW1,
t = t1)
}
tail(m.Conc)## 0 0.25 0.57 1.12 2.02 3.82 5.1
## [995,] 0 1.806089 3.0447556 3.7401290 3.5158665 2.43752291 1.833251849
## [996,] 0 0.931611 0.7765811 0.4436527 0.1736169 0.02657661 0.006996459
## [997,] 0 1.063668 0.9930997 0.6354214 0.2924619 0.06184090 0.020484321
## [998,] 0 1.530190 2.0867733 1.9640102 1.4016315 0.66871352 0.394406361
## [999,] 0 1.692536 2.6177876 2.8679642 2.3834973 1.41555760 0.968489824
## [1000,] 0 1.558275 2.1723333 2.0969694 1.5345017 0.76117664 0.461262443
## 7.03 9.05 12.12 24.37
## [995,] 1.1900928379 0.7569719786 3.805660e-01 2.447643e-02
## [996,] 0.0009352291 0.0001138159 4.634657e-06 1.314802e-11
## [997,] 0.0038716083 0.0006770512 4.783229e-05 1.222831e-09
## [998,] 0.1779102432 0.0773277104 2.179570e-02 1.392796e-04
## [999,] 0.5462534367 0.2999769340 1.206351e-01 3.183336e-03
## [1000,] 0.2167352700 0.0983138348 2.956949e-02 2.448272e-04Visualization
plot(t1, m.Conc[1,], type="l")
max.conc <- max(m.Conc)
plot(t1, m.Conc[1,], type="l", col="grey", ylim = c(0, max.conc))
for(i in 2:n){
lines(t1, m.Conc[i,], col="grey")
}
Conc1 <- S1$conc
plot(t1, Conc1)
max.conc <- max(m.Conc)
plot(t1, m.Conc[1,], type="l", col="grey", ylim = c(0, max.conc))
for(i in 2:n){
lines(t1, m.Conc[i,], col="grey")
}
points(t1, Conc1)
Uncertainty analysis of PK model in GNU MCSim
simuloR
Input file
MonteCarlo ("simmc.out",
1000,
10101010);
# Constant
IngDose = 319.992; # ingested dose (mg)
BW = 79.6; # body weight (kg)
# Distribution
Distrib (Fgutabs, Uniform, 0.8, 1); # bioavailability (-)
Distrib (kgutabs, LogUniform, 0.218, 21.8); # absortion rate (/h)
Distrib (kelim, Uniform, 0.05307864, 1.326966); # elimination rate constant (/h)
Distrib (Vdist, Uniform, 0.15603994, 3.900998); # volumes (L/kg BW)
Simulation {
Print (Conc 0.0 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 0.379285 0.201795 0.489291 0.585159 1.216130 2.03917
## 2 1 0.919264 1.060440 0.628445 0.477700 1.661360 2.896870 3.60328
## 3 2 0.915694 1.186920 0.918978 0.976574 0.859788 1.400660 1.54659
## 4 3 0.895485 0.365263 0.868337 0.470008 0.599851 1.125800 1.59117
## 5 4 0.811464 1.270980 0.570105 1.561480 0.528015 0.901527 1.08805
## 6 5 0.893335 1.618240 1.289530 0.556894 1.814340 2.600850 2.30672
## Conc_1.4 Conc_1.5 Conc_1.6 Conc_1.7 Conc_1.8 Conc_1.9
## 1 2.84207 3.229880 3.0173300 2.44664000 1.825100000 1.08586e+00
## 2 3.10623 1.391020 0.6850830 0.21800400 0.063052700 9.29577e-03
## 3 1.09045 0.319671 0.1146470 0.02214660 0.003719590 2.33483e-04
## 4 1.69653 1.176140 0.7968850 0.41413500 0.201807000 6.63069e-02
## 5 0.90691 0.399686 0.2010870 0.06834660 0.021725600 3.78039e-03
## 6 1.13856 0.164716 0.0359404 0.00330585 0.000257402 5.08044e-06
## Conc_1.10
## 1 1.02361e-01
## 2 4.23833e-06
## 3 3.13170e-09
## 4 7.57347e-04
## 5 3.50433e-06
## 6 7.13627e-13par(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 = "")
hist(log(ffpk.mtc.out$kgutabs), main = "kgutabs", 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
## 1 0.585159 1.216130 2.03917 2.84207 3.229880 3.0173300 2.44664000
## 2 1.661360 2.896870 3.60328 3.10623 1.391020 0.6850830 0.21800400
## 3 0.859788 1.400660 1.54659 1.09045 0.319671 0.1146470 0.02214660
## 4 0.599851 1.125800 1.59117 1.69653 1.176140 0.7968850 0.41413500
## 5 0.528015 0.901527 1.08805 0.90691 0.399686 0.2010870 0.06834660
## 6 1.814340 2.600850 2.30672 1.13856 0.164716 0.0359404 0.00330585
## Conc_1.8 Conc_1.9 Conc_1.10
## 1 1.825100000 1.08586e+00 1.02361e-01
## 2 0.063052700 9.29577e-03 4.23833e-06
## 3 0.003719590 2.33483e-04 3.13170e-09
## 4 0.201807000 6.63069e-02 7.57347e-04
## 5 0.021725600 3.78039e-03 3.50433e-06
## 6 0.000257402 5.08044e-06 7.13627e-13max.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", "LogUniform", "Uniform", "Uniform") # MCSim definition
q.arg <- list(list(0.8, 1.0),
list(0.218, 21.8),
list(0.05307864, 1.326966),
list(0.15603994, 3.900998))- Set experiment time-points, output variables, and conditions
outputs <- c("Ccpt")
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("BW = 79.6", "Period = 48", "IngDose = 319.992")- Modeling
model <- "models/pbtk1cpt.model"makemcsim(model)set.seed(2222)
y<-solve_mcsim(mName = model, params = parameters, vars = outputs, monte_carlo = 10000,
dist = dist, q.arg = q.arg, time = times, condition = conditions)## Starting time: 2019-12-07 20:57:26## * Created input file "sim.in".## Execute: models/mcsim.pbtk1cpt.model sim.in## Ending time: 2019-12-07 20:57:27pksim(y)
points(t1, Conc1)