Acetaminophen-PBPK model

Prepare and compile the model file

The model code needs to be prepared in the following global SA workflow. After creating the pbpk_apap.model file in the working directory, the next step is to generate the executable program (mcsim.pbpk_apap) through compile_model() function.

mName <- "pbpk_apap"
pbpk_apap_model()
compile_model(mName, application = "mcsim")

Define the parameter and its distribution

The 21 testing model parameters are defined in this part, including parameter name, probability distribution, and distributed parameter value. To prevent the computing error, the range of SULT_VmaxC and UGT_VmaxC need to adjust from \(U(0, 15)\) (Zurlinden and Reisfeld 2016) to \(U(0, 10)\) (Hsieh et al. 2018). The objects q and dist are set to the type of distribution that will use to generate the parameter matrix in GNU MCSim (for uncertainty analysis) and R (for SA).

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")
dist <- rep("Uniform", 21)
q <- rep("qunif", 21)
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))

Define additional input condition and output time and variables

To optimize the computing speed, this case only uses GNU MCSim to estimate the concentration of APAP and its metabolites glucuronide (APAP-G) and sulfate (APAP-S) 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. More definition of input schedule functions 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).

conditions <- c("mgkg_flag = 1",
                "OralExp_APAP = NDoses(2, 1, 0, 0, 0.001)",
                "OralDose_APAP_mgkg = 20.0")
vars <- c("lnCPL_APAP_mcgL", "lnCPL_AG_mcgL", "lnCPL_AS_mcgL")
times <- seq(0.1, 0.5, 1, 2, 3, 4, 6, 8, 12)

Uncertainty analysis

We apply uncertainty analysis through the solve_mcsim() and visualize the result by pksim() function. Some example data are included in the pksensi with experiment time (h) and concentration (mg/L).

head(APAP)

In the setting condition of simulation, The relative and absolute error tolerance (rtol & atol) were set to 1e-7 and 1e-9, respectively, to prevent the computing error. The Monte Carlo simulation is run for 1000 iteration as the assignment of monte_carlo. The input file (‘sim.in’) and output file (‘simmc.out’) will be generated under the standard ASCII format.

set.seed(1111)
out <- solve_mcsim(mName = mName, params = params, vars = vars,
                   monte_carlo = 1000, dist = dist, q.arg = q.arg, 
                   time = times, condition = conditions, 
                   rtol = 1e-7, atol = 1e-9)

par(mfrow = c(1,3), mar = c(4,4,1,1))
pksim(out, xlab = "Time (h)", ylab = "Conc. (ug/L)", main = "APAP")
points(APAP$Time, log(APAP$APAP * 1000))
pksim(out, vars = "lnCPL_AG_mcgL", xlab = "Time (h)", main = "APAP-G", 
      ylab = " ", legend = FALSE)
points(APAP$Time, log(APAP$AG * 1000))
pksim(out, vars = "lnCPL_AS_mcgL", xlab = "Time (h)", main = "APAP-S", 
      ylab = " ", legend = FALSE)
points(APAP$Time, log(APAP$AS * 1000))

Here shows the coverage checks of prior PBPK model predictions with calibrated APAP data. For parent compound, all data points are located in the simulated interval of 25-75%. Through this result, we can determine that the simulated outputs can accurately generate the same concentration profile as the in-vivo experiment under the setting of parameter ranges for APAP. The simulated result of metabolites APAP-G shows the different pharmacokinetic profile with experiment data. However, all data points are located in the simulated interval.

Generate parameter matrix

In global SA, we have to additionally generate the parameter matrix from the eFAST method. The current setting uses 512 sample size with 10 replication.

set.seed(1234)
x <- rfast99(params = params, n = 512, q = q, q.arg = q.arg, replicate = 10) 

Conduct the global SA

To conduct the global SA with GNU MCSim and pksensi, the input file with given “setpoint” condition should be generated before modeling. The file can create by generate_infile function. The solve_mcsim can also automatically create the input file and compute the output.

out <- solve_mcsim(x, mName = mName,
                   params = params, 
                   time = times, 
                   vars = vars,
                   condition = conditions, 
                   rtol = 1e-7, atol = 1e-9)

Visualization and decision

The plotting function can create the result of time-dependent sensitivity measurement to determine the parameter impact on model output over time.

plot(out, vars = "lnCPL_APAP_mcgL")

In addition, through using the check, the parameter with sensitivity and convergence indices over the given condition can be preliminary detected for all output variables. Based on our previous study, we proposed the heatmap visualization approach heat_check to distinguish “influential” and “non-influential” parameters with a “cut-off” point. Through the given argument order, we can select the specific order of sensitivity measurement that we’re interested in.

heat_check(out, order = "total order", show.all = T)

In the default setting, the heat_check can only show the influential parameters. The argument show.all is used to show all results. Adding the index = "CI" in the function can further investigate the convergence index. Based on the current setting of sampling size, most parameters cannot reach the acceptable criteria of convergence. Therefore, a higher number of sampling is necessary. The sample size of convergence in the current PBPK model is 8,192 (Hsieh et al. 2018). However, based on the current sample size we still can find 6 parameters that can be an important parameter for the plasma APAP concentration.

heat_check(out, index = "CI", order = "total order")