---
title: "Acetaminophen-PBPK model"
author: "Nan-Hung Hsieh"
date: "`r Sys.Date()`"
bibliography: references.bib
output: 
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{Acetaminophen-PBPK model}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
library(pksensi)
#mcsim_install(version = "6.1.0")

knitr::opts_chunk$set(
  eval = F,
  dev = "png",
  dpi = 200,
  fig.align = "center",
  fig.width = 7,
  out.width = "90%",
  fig.height = 4,
  comment = "#>"
)
```

The aim of this vignette is to reproduce our previous published [@fphar201800588] 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 [@s13318-015-0253-x]. The descriptions of each parameter and the sampling ranges are list in Table 1.

```{r}
mcsim_install(version = "6.1.0")
```

```{r, echo=F}
#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)

r <- 1.96 # exp(1.96)/exp(-1.96) ~ 50

x <- c("Tg", "Tp", "CYP_Km", "CYP_VmaxC",
       "SULT_Km_apap","SULT_Ki","SULT_Km_paps","SULT_VmaxC",
       "UGT_Km","UGT_Ki","UGT_Km_GA","UGT_VmaxC",
       "Km_AG","Vmax_AG","Km_AS","Vmax_AS",
       "kGA_syn","PAPS_syn", "CLC_APAP","CLC_AG","CLC_AS")
y <- c("Gatric emptying time constant",
       "GI perfusion time constant",
       "Cytochrome P450 metabolism, Km",
       "Cytochrome P450 metabolism, VMax",
       "Sulfation pathway acetaminophen, Km",
       "Sulfation pathway substrate inhibition, Ki",
       "Sulfation pathway PAPS, Km",
       "Sulfation pathway acetaminophen, Vmax",
       "Glucronidation pathway acetaminophen, Km",
       "Glucronidation pathway substrate inhibition, Ki",
       "Glucronidation pathway GA, Km",
       "Glucronidation pathway acetaminophen, Vmax",
       "APAP-G hepatic transporter, Km",
       "APAP-G hepatic transporter, Vmax",
       "APAP-S hepatic transporter, Km",
       "APAP-S hepatic transporter, Vmax",
       "UDPGA synthesis",
       "PAPS synthesis",
       "APAP clearance",
       "APAP-G clearance",
       "APAP-S clearance")
z <- c("$h$", "$h$", "$\\mu{M}$", "$\\mu{mole}/h\\cdot{BW}^{0.75}$",
       "$\\mu{M}$", "$\\mu{M}$", "$-$", "$\\mu{mole}/h\\cdot{BW}^{0.75}$",
       "$\\mu{M}$", "$\\mu{M}$", "$-$", "$\\mu{mole}/h\\cdot{BW}^{0.75}$",
       "$\\mu{M}$", "$\\mu{mole}/h$", "$\\mu{M}$", "$\\mu{mole}/h$",
       "$1/h$", "$1/h$",
       "$L/h\\cdot{BW}^{0.75}$", "$L/h\\cdot{BW}^{0.75}$", "$L/h\\cdot{BW}^{0.75}$")
min <- c(round(Tg-r, 3), round(Tp-r, 3), round(CYP_Km-r), round(log(0.14), 3),
         round(SULT_Km_apap-r, 3), round(SULT_Ki-r, 3), round(SULT_Km_paps-r), log(1),
         round(UGT_Km-r, 3), round(UGT_Ki-r, 3), round(UGT_Km_GA-r), log(1),
         round(Km_AG-r, 3), round(log(1.09e3), 3), round(Km_AS-r), round(log(1.09e3), 3),
         log(1), log(1), round(log(2.48e-3), 3), round(log(2.48e-3), 3), round(log(2.48e-3), 3))
max <- c(round(Tg+r, 3), round(Tp+r, 3), round(CYP_Km+r), round(log(2900), 3),
         round(SULT_Km_apap+r, 3), round(SULT_Ki+r, 3), round(SULT_Km_paps+r), round(log(22026), 3),
         round(UGT_Km+r, 3), round(UGT_Ki+r, 3), round(UGT_Km_GA+r), round(log(22026), 3),
         round(Km_AG+r, 3), round(log(3.26e6), 3), round(Km_AS+r), round(log(3.26e6), 3),
         round(log(4.43e5), 3), round(log(4.43e5),3), 
         round(log(2.718), 3), round(log(2.718), 3), round(log(2.718), 3))

df <- data.frame(x, y, z, min, max)
names(df) <- c("Parameter","Description", "Unit", "Min", "Max")

#if (require(kableExtra)) {
  knitr::kable(df, format = 'html', align=c("l","l","l", "c", "c"), 
               caption = "Table 1 Description of sampling range of model parameter") 
  #%>% kableExtra::add_footnote(c("The parameter valur are showed in log-transformed scale."), notation = "number")    
#}

```

Same as the example of one-compartment PK 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 under log-scaled (approximate 54.6 times difference between minimum and maximum in natural scaled). The nominal values of informative model parameters were defined as:

```{r}
# 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, wide range of parameter value might cause the computing error when solving the differential equation. 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()` and `solve_mcsim()` provide the arguments of `rtol` and `atol` that adjust the error tolerance to prevent the unwanted error. However, the modification will decrease the computing speed. Therefore, the alternative method to prevent this issue is to detect the crucial parameter range that causes the problem. Also, setting the maximum number of steps to higher value instead of using the default value (500) in **GNU MCSim** can prevent this problem (internally defined). The maximum number of step is set to 5000 in this case. Here we separate the global SA of APAP-PBPK model process to several steps.

## 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. 

```{r}
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)$ [@s13318-015-0253-x] to $U(0, 10)$ [@fphar201800588]. 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).

```{r}
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).

```{r}
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).

```{r}
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. 

```{r}
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)
```

```{r, fig.height=3.5, fig.width=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. 

```{r}
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.

```{r}
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. 

```{r, fig.height=8, fig.width=8}
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. 

```{r, fig.height=5}
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 [@fphar201800588]. However, based on the current sample size we still can find 6 parameters that can be an important parameter for the plasma APAP concentration.

```{r, fig.height=5}
heat_check(out, index = "CI", order = "total order")
```


## References
