The aim of this section is to reproduce our previous published (Hsieh et al. 2018) result of global sensitivity analysis for acetaminophen PBPK model through pksensi. The model codes are included in this package and can be generated through pbpk_apap_model(). We applied the global sensitivity analysis workflow to the original published model with 21 model parameters (Zurlinden and Reisfeld 2016). The descriptions of each parameter and the sampling ranges are list in Table 1.
Same as the example of one-compartment PBTK model. The model parameter and the corresponding sampling range should be defined to create the parameter matrix. Previously, the probability distributions of model parameters were set to either truncated normal or uniform distribution when the parameters have informative prior information or not. To rapidly reach the acceptance convergence, we apply uniform distribution for all testing parameters. The ranges of informative parameters are set to 1.96-times difference for single side (approximate 54.6 times difference between minimum and maximum) under log-scaled. The nominal values of informative model parameters were defined as:
# Nominal value
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 Generally, The wide range of parameter value might cause the computational error in the solver. One of the effective ways to prevent this problem is to adjust the value of relative and absolute error tolerance to control the error appearance by resetting these parameters in a lower value. The generate_infile() provide the arguments of rtol and atol that can be adjusted to prevent the unwanted error. However, the modification will slow down the computational speed. Therefore, the alternative method to prevent the computational error is to detect the crucial parameter range that causes the problem. Also, setting the maximum number of (internally defined) steps to higher value instead of using the default value (500) can prevent this problem. The maximum number of step is set to 5000 in this case.
In this test case, we adjusted the range of SULT_VmaxC and UGT_VmaxC from U(0, 15) to U(0, 10). The relative and absolute error tolerance were set to 1e-7 and 1e-9, respectively, to prevent the computational error in MCSim,
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))
times <- seq(from = 0.1, to = 12.1, by = 0.2)
set.seed(1234)
x <- rfast99(params = params, n = 512, q = q, q.arg = q.arg, replicate = 10) After creating the pbpk_apap.model in the working directory, the next step is to generate the executable files (mcsim.pbpk_apap) through compile_model().
To improve the computational speed, this case only uses MCSim to estimate the concentration of acetaminophen (APAP) and its metabolites glucuronide (AG) and sulfate (AS) in plasma. The setting oral dose of APAP is 20 mg/kg in this example. Generally, the input dosing method can be defined through the condition argument. Since the unit of the given dose is mg/kg, the mgkg_flag is set to 1 to declare the statement. More definition of input can be found in the section of input functions in GNU MCSim User’s Manual (https://www.gnu.org/software/mcsim/mcsim.html#Input-functions).
vars <- c("lnCPL_APAP_mcgL", "lnCPL_AG_mcgL", "lnCPL_AS_mcgL")
conditions <- c("mgkg_flag = 1",
"OralExp_APAP = NDoses(2, 1, 0, 0, 0.001)",
"OralDose_APAP_mgkg = 20.0")
generate_infile(params = params,
vars = vars,
time = times,
condition = conditions,
rtol = 1e-7, atol = 1e-9)
system.time(y <- solve_mcsim(x, mName = mName,
params = params,
vars = vars,
time = times,
condition = conditions,
generate.infile = F))
tell2(x,y)The plotting function can output the result of time-dependent sensitivity measurement to determine the parameter impact on model output over time (Figure 1).
The uncertainty analysis is a crucial step before model calibration. We can apply uncertainty analysis through the pksim() by the given name of output variables (Figure 2). Through this visualization approach, we can recognize whether the simulated outputs can accurately simulate the same concentration profile as the in-vivo experiment under the setting of parameter ranges. In Figure 2, all experiment points are included in the intervals, representing the acceptable set of the parameter range.
par(mfrow = c(2,2), mar = c(2,2,1,1), oma = c(2,2,1,1))
pksim(y, vars = "lnCPL_APAP_mcgL")
text(1, 15, "APAP",cex = 1.2)
points(APAP$Time, log(APAP$APAP * 1000))
pksim(y, vars = "lnCPL_AG_mcgL", legend = F)
text(1, 15, "AG",cex = 1.2)
points(APAP$Time, log(APAP$AG * 1000))
pksim(y, vars = "lnCPL_AS_mcgL", legend = F)
text(1, 15, "AS",cex = 1.2)
points(APAP$Time, log(APAP$AS * 1000))
mtext("Time", SOUTH<-1, line=0.4, cex=1.2, outer=TRUE)
mtext("Conc.", WEST<-2, line=0.4, cex=1.2, outer=TRUE)In addition, through using the check(), the parameter with sensitivity and convergence indices over the given condition can be easily detected for all output variables.
The check() also provides some feasible argument to specify the target output or change the cut-off value.
Based on our previous study, we proposed the heatmap visualization approach to distinguish “influential” and “non-influential” parameters with a cut-off. Through the given argument order, we can select the specific order of sensitivity measurement that we’re interested in (Figure 3 & 4).
Also, adding the index = "CI" in the function can further investigate the convergence of the sensitivity index. Based on the current setting of sampling number, most parameters cannot reach the acceptable criteria of convergence. Therefore, the higher number of sampling is necessary.
Hsieh, Nan-Hung, Brad Reisfeld, Frederic Y. Bois, and Weihsueh A. Chiu. 2018. “Applying a Global Sensitivity Analysis Workflow to Improve the Computational Efficiencies in Physiologically-Based Pharmacokinetic Modeling.” Frontiers in Pharmacology 9: 588. https://doi.org/10.3389/fphar.2018.00588.
Zurlinden, Todd J., and Brad Reisfeld. 2016. “Physiologically Based Modeling of the Pharmacokinetics of Acetaminophen and Its Major Metabolites in Humans Using a Bayesian Population Approach.” European Journal of Drug Metabolism and Pharmacokinetics 41: 267–80. https://doi.org/10.1007/s13318-015-0253-x.