diff --git a/NEWS.md b/NEWS.md index 448348f..ae8b6d7 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,12 +1,13 @@ # pksensi 1.1.2 -* Add timer in solve_fun -* Add default MCSim version = '6.1.0' in compile_model() -* Update pksim(): Can customize xlab and ylab -* Update pbtk1cpt model: Add body weight -* Revise the mcsim_version() for Windows OS -* Fix bug: Create .dll in win10 through compile_model() -* Fix bug: plot() cannot display variable name when using character +* Add timer in `solve_fun()` +* Add default MCSim `version = '6.1.0'` in `compile_model()` +* Update `pksim()`: Can customize xlab and ylab +* Update `pbtk1cpt model`: Add body weight +* Update vignettes +* Fix bug: Create `.dll` in win10 through `compile_model()` +* Fix bug: `plot()` cannot display variable name when using character +* Revise the `mcsim_version()` for Windows OS # pksensi 1.1.1 diff --git a/cran-comments.md b/cran-comments.md index 48fbab5..a17848b 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,7 +1,7 @@ ## Test environments -* local OS X install, R 3.6.0 +* local macOS Mojave 10.14.5, R 3.6.0 * local Ubuntu 18.04, R 3.6.0 -* local Windows 7, R 3.5.0 +* local Windows 10, R 3.6.0 * Ubuntu 14.04.5 LTS(on travis-ci), R 3.6.0 * Windows Server 2012 R2 x64 (appveyor), R 3.6.0 @@ -10,9 +10,14 @@ There were no ERRORs, WARNINGs and NOTE. ## r-hub builder results -### pksensi 1.1.1: OK +### pksensi 1.1.2: OK -Build ID: pksensi_1.1.1.tar.gz-8c1b6d8236c349c3919c7864b8cdaae0 +Build ID: pksensi_1.1.2.tar.gz-bfe9da0c549d4fe08ced23fff281ca5a Platform: Windows Server 2008 R2 SP1, R-devel, 32/64 bit -Submitted: 3 minutes 47.3 seconds ago -Build time: 3 minutes 21.5 seconds +Submitted: 3 minutes 59.2 seconds ago +Build time: 3 minutes 34.6 seconds + +Build ID: pksensi_1.1.2.tar.gz-03b57ebe8ac646eba1595ddc6c80b67e +Platform: Ubuntu Linux 16.04 LTS, R-release, GCC +Submitted: 20 minutes 6.4 seconds ago +Build time: 20 minutes 0.8 seconds diff --git a/vignettes/pbpk_apap.Rmd b/vignettes/pbpk_apap.Rmd index 5347346..0721b73 100644 --- a/vignettes/pbpk_apap.Rmd +++ b/vignettes/pbpk_apap.Rmd @@ -22,15 +22,12 @@ knitr::opts_chunk$set( ) ``` -## Uncertainty and sensitivity analysis - The aim of this section 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(mxstep = 5000) # Be sure to set the mxstep to 5000 ``` - ```{r, echo=F, eval=F} #Nominal value Tg <- log(0.23) @@ -102,7 +99,7 @@ names(df) <- c("Parameter","Description", "Unit", "Min", "Max") ``` -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: +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, eval=F} # Nominal value @@ -121,9 +118,21 @@ 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. +Generally, wide range of parameter value might cause the computing error when solving the differentiial 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 seperate 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, eval=F} +mName <- "pbpk_apap" +pbpk_apap_model() +compile_model(mName, application = "mcsim") +``` + +### Define the parameter and its distribution -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, +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, eval=F} params <- c("lnTg", "lnTp", "lnCYP_Km","lnCYP_VmaxC", @@ -131,112 +140,109 @@ params <- c("lnTg", "lnTp", "lnCYP_Km","lnCYP_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.), +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(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 = 1024, q = q, q.arg = q.arg, replicate = 5) # Used n = 8192 and rep = 10 in published papaer + 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)) ``` -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()`. - -```{r, eval=F} -mName <- "pbpk_apap" -pbpk_apap_model() -compile_model(mName, application = "mcsim", version = "6.1.0") -``` +### Define additional input condition and output time and variables -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). - -**NOTE:** The `mxstep` was set at `5000` to prevent the computational error. In addition, use `rtol = 1e-7` and `atol = 1e-9` in error tolerance. +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, eval=F} -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") -system.time(out <- solve_mcsim(x, mName = mName, - params = params, - vars = vars, - time = times, - condition = conditions)) +vars <- c("lnCPL_APAP_mcgL", "lnCPL_AG_mcgL", "lnCPL_AS_mcgL") +times <- seq(from = 0.1, to = 12.1, by = 0.2) ``` +### 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). -The plotting function can output the result of time-dependent sensitivity measurement to determine the parameter impact on model output over time (Figure 1). - -```{r, eval=F, fig.height=8, fig.width=8, fig.cap = 'Figure 1. '} -plot(out, vars = "lnCPL_AG_mcgL") +```{r, eval=F} +head(APAP) ``` -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. +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 computating 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, eval=F} +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, eval=F, fig.cap = 'Figure 2. The range of model simulation based on parameter distribution'} -data(APAP) -par(mfrow = c(2,2), mar = c(2,2,1,1), oma = c(2,2,1,1)) -pksim(out, vars = "lnCPL_APAP_mcgL") -text(1, 15, "APAP",cex = 1.2) +```{r, eval=F, fig.cap = '**Figure 1.** Coverage checks of prior PBPK model predictions with calibrated APAP data', 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", legend = F) -text(1, 15, "AG",cex = 1.2) +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", legend = F) -text(1, 15, "AS",cex = 1.2) +pksim(out, vars = "lnCPL_AS_mcgL", xlab = "Time (h)", main = "APAP-S", + ylab = " ", legend = FALSE) 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. +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, eval=F} -check(out) +set.seed(1234) +x <- rfast99(params = params, n = 512, q = q, q.arg = q.arg, replicate = 10) ``` -The `check()` also provides some feasible argument to specify the target output or change the cut-off value. + +### 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, eval=F} -check(out, vars = "lnCPL_APAP_mcgL", SI.cutoff = 0.1, CI.cutoff = 0.1) +out <- solve_mcsim(x, mName = mName, + params = params, + time = times, + vars = vars, + condition = conditions, + rtol = 1e-7, atol = 1e-9) ``` -## Heatmap visualization combined with an index “cut-off” -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). +### 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, eval=F, 'Figure 3: Heatmap of sensitivity index for interaction'} -heat_check(out, order = "interaction") +```{r, eval=F, fig.cap = '**Figure 2.** Time-dependent SI of the plasma APAP concentration estimated from APAP-PBPK model during 12-hr time period post intake.', fig.height=8, fig.width=8} +plot(out, vars = "lnCPL_APAP_mcgL") ``` -```{r, eval=F, 'Figure 4: Heatmap of convergence index for total order'} +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, eval=F, fig.cap = '**Figure 3.** Heat map representation of time-dependent total SI computed for global SA method', fig.height=5} heat_check(out, order = "total order") ``` -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. +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, eval=F, fig.height=9, 'Figure 5: Heatmap of convergence index'} -heat_check(out, index = "CI") +```{r, eval=F, fig.cap = '**Figure 4.** Heat map representation of time-dependent total CI computed for global SA method', fig.height=5} +heat_check(out, index = "CI", order = "total order") ``` + ## References diff --git a/vignettes/pbtk1cpt.Rmd b/vignettes/pbtk1cpt.Rmd index cce540c..b98b38f 100644 --- a/vignettes/pbtk1cpt.Rmd +++ b/vignettes/pbtk1cpt.Rmd @@ -1,5 +1,5 @@ --- -title: "One-compartment PBTK model" +title: "One-compartment PK model" author: "Nan-Hung Hsieh" date: "`r Sys.Date()`" bibliography: references.bib @@ -7,7 +7,7 @@ output: rmarkdown::html_vignette: fig_caption: yes vignette: > - %\VignetteIndexEntry{One-compartment PBTK model} + %\VignetteIndexEntry{One-compartment PK model} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- @@ -20,7 +20,7 @@ knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.height=4, - fig.width=6 + fig.width=7 ) fn = local({ # not used function for fig caption @@ -40,40 +40,54 @@ In this example, We use a simple, one-compartment PK model from `httk` package [ $$\frac{dA_{gutlumen}}{dt} = -k_{gutabs} \cdot A_{gutlumen} + g(t)$$ $$\frac{dC_{rest}}{dt} = \frac{k_{gutabs}}{V_{dist}}-k_{elim} \cdot C_{rest}$$ -where $A_{gutlumen}$ is the state variable that describes the quantity of compound in gut lumen (mol), $k_{gutabs}$ is the absorption rate constant that describes the chemical absorption from the gut lumen into gut tissue through first-order processes (/h), $V_{dist}$ is the volume of distribution (L), and $k_{elim}$ is the elimination rate constant (/h), which is equal to the total clearance divided by the volume of distribution. The time-dependent function $g(t)$ is used to describe the oral dosing schedule. $C_{rest}$ is the chemical concentration in plasma that can be used to compare with observed results in a pharmacokinetic experiment (mol/L). +where $A_{gutlumen}$ is the state variable that describes the quantity of compound in the gut lumen (mg) and $A_{rest}$ is the quantity of compound in rest of body and blood (mg). The parameter $k_{gutabs}$ is the absorption rate constant that describes the chemical absorption from the gut lumen into gut tissue through first-order processes (/h) and $k_{elim}$ is the elimination rate constant (/h), which is equal to the total clearance divided by the volume of distribution. The time-dependent function $g(t)$ is used to describe the oral dosing schedule. -## Model implementations with R deSolve package +The concentration of the chemical in the rest of body and blood ($C_{rest}$, mg/L) can be calculated as + +$$ C_{rest} = A_{rest} / V_{dist} \cdot BW$$ -In the beginning, we need to pre-install GNU MCSim [@JSSv002i09] and related compiler to let us generate .c file and executable file. The GNU MCSim can be installed by following the instruction in GNU MCSim’s manual on https://www.gnu.org/software/mcsim/mcsim.html or using the build-in function `mcsim_install()` in `pksensi`. The GNU compiler is necessary for users that use Linux or MacOS. For Windows users, you should install Rtools on https://cran.r-project.org/bin/windows/Rtools/ and use `Sys.setenv()` to set the working path of compiler. The `Sys.which("gcc")` and `system('g++ -v')` can check whether we can run compiler correctly. +where $V_{dist}$ is the volume of distribution (L/kg BW) and $BW$ is the body weight (kg). The $C_{rest}$ can also be seen as the chemical concentration in plasma that can be further used to compare with observed results in a pharmacokinetic experiment. The bioavailability is assumed to be 100% in this model. + +## Model implementations with R deSolve package -We first implemented this model in R by compiling the file written in C. pksensi allows users to select the preferred method to solve the pharmacokinetic model, either with the `deSolve` package or with GNU MCSim through the compile function. However, running model under GNU MCSim native code can have faster speed to obtain the model outputs. +To start, we implemented the one-compartment PK model in R by compiling the file written under the **GNU MCSim**'s model format. **pksensi** allows users select the preferred method to solve the pharmacokinetic model, either with the **deSolve** [@JSSv033i09] package or with **GNU MCSim** [@JSSv002i09] through the compile function. Running model under **GNU MCSim** native code can have a faster speed to obtain the model outputs. This section mainly focus on how to conduct global SA with **deSolve** package. Then compare the computing time with **GNU MCSim**. -The following R script can download and compile the model description file (`pbtk1cpt.model`) and use `deSolve` package [@JSSv033i09] to solve ordinary differential equations in our model. The example model code of one-compartment PBTK model is available with pksensi package: +The following R script will download the one-compartment PK model from https://github.com/nanhung/pksensi/blob/master/tests/pbtk1cpt.model. The model mainly includes two state variables that are the quantity of compound in the gut lumen (`Agutlument`) and the rest of body (`Acompartmant`). The `Ametabolized` is the quantity of compound transform and metabolize through hepatic clearance. The `AUC` is the calculated area under the curve of the rest body ($mg \cdot h/L$). ```{r, eval=F} pbtk1cpt_model() cat(readLines("pbtk1cpt.model"), sep = "\n") ``` -Then, use `compile_model()` to generate the executable files (`pbtk1cpt.dll` on Windows or `pbtk1cpt.so` on other systems) and R file (`pbtk1cpt_inits.R`) with default input parameters and initial state settings with the definition of `application = "R"`. +The code must first be compiled to run the model. After create the model file, we can use `compile_model` to generate the file that has dynamic-link library (DLL) or share object (SO) extention and can be linked dynamically into an R session (`pbtk1cpt.dll` on Windows or `pbtk1cpt.so` on other systems) and R file (`pbtk1cpt_inits.R`) with default input parameters and initial state settings with the definition of `application = "R"`. -```{r, eval=F} +```{r eval=F, message=FALSE, warning=FALSE} mName <- "pbtk1cpt" compile_model(mName, application = "R") +``` + +The `pbtk1cpt_inits.R` file includes `initParms` and `initStates` functions and `Outputs` variable. The created function have default value of model parameters and initial conditions that can further use to customize in the simulation. + +```{r, eval=F} source(paste0(mName, "_inits.R")) ``` -The parameter values and initial states can be customized to specify the properties and schedule for the given dosing scenario. +The parameter values and initial states can be customized to the specific condition. It can also schedule for the given dosing scenario. In the current setting, we assumed the initial condition of the intake dose to be 1000 mg. We can use `initParms` and `initStates` functions to customize the parameter values and the initial state that will be used in the following modeling and SA. These additional functions are generated when compiling the model file. + +We used the corresponding parameter value of (APAP) in this example. The parameter value can be generated from `parameterize_1comp` function in \pkg{httk} package, which includes physico-chemical properties for over 500 chemicals. The given value of `vdist`, `ke`, and `kgutabs` in \pkg{httk} are 1.1 (L/kg BW), 0.23 (/h), and 2.18 (/h), respectively. The body weight is assumed to be 70 (kg). ```{r, eval=F} parms <- initParms() -parms["vdist"] <- 0.74 -parms["ke"] <- 0.28 +parms["vdist"] <- 1.1 +parms["ke"] <- 0.23 parms["kgutabs"] <- 2.18 +parms["BW"] <- 70 initState <- initStates(parms=parms) initState["Agutlument"] <- 10 ``` +Here shows the given parameter value (`parms`), initial state condition (`initStat`), and the output variable (`Outputs`) that need to specify in model solving. + ```{r, eval=F} parms ``` @@ -86,109 +100,148 @@ initState Outputs ``` -In the current setting, we assumed the initial condition of the intake chemical to be 10 mol. The `initParms` and `initStates` functions were used to customize the parameter values and the initial state that will be used in the `solve_fun` function. These parameter value can be adopted from the `httk` package, which includes physico-chemical and drug biological properties for 553 chemicals. In this case, we used the parameter value of theophylline in this example. The given `vdist`, `ke`, and `kgutabs` are 0.74, 0.28, and 2.18, respectively. - -Through `ode` function in `deSolve` package, we can visualize the pharmacokinetic according to the given parameter conditions such as time points (`times`): +Using the `ode` function in **deSolve** package, we can visualize the pharmacokinetic profile according to the given parameter baseline and the time points (`t`): ```{r, eval=F} -times <- seq(from = 0.01, to = 24.01, by = 1) -y <- deSolve::ode(initState, times, func = "derivs", parms = parms, - dllname = mName, initfunc = "initmod", nout = 1, outnames = Outputs) +t <- seq(from = 0.01, to = 24.01, by = 1) +y <- deSolve::ode(initState, t, func = "derivs", parms = parms, + dllname = mName, initfunc = "initmod", + nout = 1, outnames = Outputs) ``` -```{r, eval=F, fig.cap = 'Figure 1. Simulation results of one-compartment PBTK model.'} +```{r, eval=F, fig.cap = '**Figure 1.** Simulation results of one-compartment PK model.'} plot(y) ``` -To conduct sensitivity analysis for the parameters in one-compartment pharmacokinetic model in this case, we want to quantify the impact of these three parameters on the chemical concentration in plasma during 24-hour time period post intake. We assume a uniform distribution for the estimate for each parameter with the coefficient of uncertainty within 50%. The parameter ranges are assumed to be (0.37, 1.12) for `vdist`, (0.0058, 0.0174) for `ke`, and (0.045, 0.136) for `kgutabs`. The sample number determines the robustness of the result of sensitivity analysis. Higher sample numbers can generate narrow confidence intervals for sensitivity measurements across different replications. However, they might cause heavy computational burden for complex models. Here we use a sample number of 400 with 20 replications: +We conduct global SA for the four parameters in one-compartment PK model to investigate the parameter impact on the plasma concentration during the 24-hr time period post intake. The distribution of parameter is taken to be uniform with bounds corresponding to 50% and 200% of the nominal value. Therefore, the parameter ranges are assumed to be 0.55 and 2.2 L/kg BW for `vdist`. The `ke` are ranged from 0.12 to 0.46 /h, corresponding to half-times of 1.5 and 5.8 hr. The `ka` are ranged 1.09 to 4.36 /h, corresponding to half-times of 0.32 and 0.64 hr. The `BW` is assumed to a normal distribution with mean = 60 kg and sd = 5 kg. The number of sample size determines the robustness of the result of SA. Higher number of sample size lead to narrower confidence intervals for sensitivity measurements across different replications. However, it will take a longer time in computation. Here we use a sample size of 200 with 20 replications. + +```{r, eval=F} +q <- c("qunif", "qunif", "qunif", "qnorm") +q.arg <- list(list(min = parms["vdist"] / 2, max = parms["vdist"] * 2), + list(min = parms["ke"] / 2, max = parms["ke"] * 2), + list(min = parms["kgutabs"] / 2, max = parms["kgutabs"] * 2), + list(mean = parms["BW"], sd = 5)) +params <- c("vdist", "ke", "kgutabs", "BW") +``` + +Through `rfast99` function, a S3 object with class 'rfast99' will be created. The `set.seed` can use to reproduce the same parameter matrix in the random sampling. ```{r, eval=F} -LL <- 0.5 -UL <- 1.5 -q <- "qunif" -q.arg <- list(list(min = parms["vdist"] * LL, max = parms["vdist"] * UL), - list(min = parms["ke"] * LL, max = parms["ke"] * UL), - list(min = parms["kgutabs"] * LL, max = parms["kgutabs"] * UL)) set.seed(1234) -params <- c("vdist", "ke", "kgutabs") x <- rfast99(params, n = 200, q = q, q.arg = q.arg, replicate = 20) ``` -Because the pharmacokinetic model is being used to describe a continuous process for the chemical concentration over time, the sensitivity measurements can also show the time-dependent relationships for each model parameter. Here we define the output time points to examine the change of the parameter sensitivity over time. To solve the pharmacokinetic model through deSolve, we need to provide the details of the argument: +The generated parameters are stored as a 3-D array under the named `a`, with the dimension of sample size, number of replication, and number of parameter, respectively. ```{r, eval=F} -out <- solve_fun(x, times, initState = initState, outnames = Outputs, dllname = mName) +dim(x$a) +``` + +The sample number is 200 with 4 model parameters, which generates 800 model evaluations. The replication is set to 20. Therefore, the totol of 16,000 parameter sequence will be used to compute the corresponding outputs. + +```{r, eval=F, fig.cap = '**Figure 2.** The parameter search curve generated by eFAST method. The x-axis is the sample number and y-axis is the parameter value for `vdist`, \`ke`, `kgutabs`, and \`BW` (top to down), respectively. The bar plot shows the distribution of sampling parameter.', fig.height=5} +par(mfrow=c(4,4),mar=c(0.8,0.8,0.8,0),oma=c(4,4,2,1), pch =".") +for (j in c("vdist", "ke", "kgutabs", "BW")) { + if ( j == "BW") { + plot(x$a[,1,j], ylab = "BW") + } else plot(x$a[,1,j], xaxt="n", ylab = "") + for (i in 2:3) { + if ( j == "BW") { + plot(x$a[,i,j], ylab = "", yaxt="n") + } else plot(x$a[,i,j], xaxt="n", yaxt="n", ylab = "") + } + hist <- hist(x$a[,,j], plot=FALSE, + breaks=seq(from=min(x$a[,,j]), to=max(x$a[,,j]), length.out=20)) + barplot(hist$density, axes=FALSE, space=0, horiz = T, main = j) +} +mtext("Model evaluation", SOUTH<-1, line=2, outer=TRUE) +``` + +Because the PK model is being used to describe a continuous process for the chemical concentration over time, the sensitivity measurements, therefore, have the time-dependent relationships for each model parameter. Here we define the output time points (`t`) to examine the change of the parameter sensitivity over time. To solve the model through **deSolve**, we need to provide the details of the argument, which include time (`t`), initial conditions of state variable (`initState`), output variables (`outnames`), and name of the shared library (`mName`, without extension). To create the time-dependent sensitivity measurement, we set the time duration from 0.01 to 24.01 hours with the time segment of 1 hour as the above definition in `ode` function in this example. The initial time point should avoid 0 to prevent computational error in misconduct. The `outnames`, `dllname`, are based on the arguments from the `ode` function in **deSolve** package. The details of the model structure and these arguments are defined in `pbtk1comp.c`. and `pbtk1comp_inits.R` files. + +```{r, eval=F} +outputs <- c("Ccompartment", "Ametabolized") +out <- solve_fun(x, time = t, initState = initState, + outnames = outputs, dllname = mName) ``` -To create the time-dependent sensitivity measurement, we set the time duration from 0.01 to 24.01 hours in the example. The `initParmsfun` is used to generate the sampling value for each parameter. The `outnames`, `dllname`, `func`, `initfunc` are based on the arguments from the ode function in `deSolve` package. The details of the model structure and these arguments are defined in `pbtk1comp.c`. and `pbtk1comp_inits.R`. Finally, the `tell2` function is used to integrate the parameter values and the output results of numerical analysis that were generated and stored in variables x and y. The result of object x is an object of rfast99, which has specific `print`, `plot`, and `check` method. The print function gives the sensitivity and convergence indices for main, interaction, and total order at each time point. In addition to print out the result of sensitivity analysis, the more efficient way to distinguish the influence of model parameter is to visualize them. The time-dependent sensitivity indices are shown in Figure 2. +The output result `out` is an S3 object of `rfast99` as well, which can link with `print`, `plot`, and `check` method to examine the sensitivity measurements. The `print` function gives the sensitivity and convergence indices for main, interaction, and total order at each time point. In addition to print out the result of SA, the more efficient way to distinguish the influence of model parameter is to visualize these indices. -```{r, eval=F, fig.cap = 'Figure 2. Time-dependent sensitivity indices of the plasma concentration estimated from one-compartment PBTK model during 24 hour time period intake.'} +```{r, eval=F, fig.cap = '**Figure 3.** Time-dependent SI of the plasma concentration estimated from one-compartment PK model during 24-hr time period intake. The solid line represent the total (black) and first (red) order SI with 95\\% confidence interval (polygon). The dash line is the cut-off with default value of 0.05.'} plot(out) ``` -Here, we can find that `vdist` and `ke` are dominating the plasma concentration in the before and after 5-hour post chemical intake, respectively, representing that the elimination is a key parameter to dominate the plasma concentration. Besides, the `kgutabs` only plays a crucial role to determine the plasma concentration in the first hour. The relationship between concentration and the parameters can be plotted as follow (Figure 3): +The SI has computed range from 0 (no impact) to 1 (high impact) and represent the contribution percentage of output variance under the given parameter distributions. -```{r, eval=F, fig.cap = 'Figure 3. The relationship between model parameter and estimated concentration under the time-point of 0.01, 2.01, and 24.01 hr'} -par(mfrow = c(3,3), mar = c(2,2,2,2), oma = c(2,2,1,1)) -plot(x$a[,1,"vdist"], out$y[,1,"0.01",], main = "vdist") -text(1, .7, "t=0.01",cex = 1.2) -plot(x$a[,1,"ke"], out$y[,1,"0.01",], main = "ke") -plot(x$a[,1,"kgutabs"], out$y[,1,"0.01",], main = "kgutabs") -plot(x$a[,1,"vdist"], out$y[,1,"2.01",]) -text(1, 18, "t=2.01",cex = 1.2) -plot(x$a[,1,"ke"], out$y[,1,"2.01",]) -plot(x$a[,1,"kgutabs"], out$y[,1,"2.01",]) -plot(x$a[,1,"vdist"], out$y[,1,"24.01",]) -text(1, .7, "t=24.01",cex = 1.2) -plot(x$a[,1,"ke"], out$y[,1,"24.01",]) -plot(x$a[,1,"kgutabs"], out$y[,1,"24.01",]) -mtext("parameter", SOUTH<-1, line=0.4, outer=TRUE) -mtext("Ccompartment", WEST<-2, line=0.4, outer=TRUE) +Here, we can see that `vdist` and `ke` are dominating the plasma concentration before and after 5-hr post chemical intake, respectively. Besides, the `kgutabs` only plays a crucial role to determine the plasma concentration in the first hour. However, the current result is only based on the setting distribution of model parameters for APAP. The different given input conditions (e.g., range of parameter uncertainty, chemical dependent parameter value) can change the result (result not shown). + +The default output in the plotting is setting at the first variable. To exam the time-sependent SI of other variables, such as `Ametabolized` in this case, we need to assign the variable name `vars = "Ametabolized"` in `plot` function. + +```{r, eval=F, fig.cap = '**Figure 4.** Time-dependent SI of the plasma concentration estimated from one-compartment PK model during 24-hr time period intake.'} +plot(out, vars = "Ametabolized") ``` -The x is a list of class "rfast99", containing all the input arguments detailed before and the calculated sensitivity indices of first order (`mSI`), interaction (`iSI`), and total order (`tSI`). The convergence indices are also stored in the list named `mCI`, `iCI`, and `tCI`. The parameter values are stored in an array `x$a` with c(model evaluation, replication, parameters). +The amount of metabolized is also determined by parameter `ke`. Same as `Ccompartment`, the `kgutabs` contribute about 30 - 40% variation of model output in the first hour. The `BW` is the least important parameter in the current analysis, and therefore can be fixed in the model fitting to data and additional applications. -```{r, eval=F} -dim(x$a) +In addition to use time-SI profile to investigate the parameter impact on model output, we can construct the relationship between parameter and the corresponding value of model output in the examination. + +```{r, eval=F, fig.cap = '**Figure 5.** The relationship between model parameter and estimated concentration at times 0.01, 2.01, 6.01, and 24.01 hr (top to bottom), respectively.'} +par(mfcol=c(4,4),mar=c(0.8,0.8,0,0),oma=c(4,4,2,1), pch = ".") +plot(x$a[,1,"vdist"], out$y[,1,"0.01",1], xaxt="n", main = "\nvdist") +plot(x$a[,1,"vdist"], out$y[,1,"2.01",1], xaxt="n") +plot(x$a[,1,"vdist"], out$y[,1,"6.01",1], xaxt="n") +plot(x$a[,1,"vdist"], out$y[,1,"24.01",1]) +for (j in c("ke", "kgutabs", "BW")){ + for (k in c("0.01", "2.01", "6.01", "24.01")){ + if (k == "0.01") { + plot(x$a[,1,j], out$y[,1,k,1], yaxt = "n", xaxt="n", main = paste0("\n", j)) + } else if (k == "24.01") { + plot(x$a[,1,j], out$y[,1,k,1], yaxt = "n") + } else plot(x$a[,1,j], out$y[,1,k,1], xaxt = "n", yaxt = "n") + } +} +mtext("Parameter", SOUTH<-1, line=2, outer=TRUE) +mtext("Ccompartment", WEST<-2, line=2, outer=TRUE) ``` -In addition, the output are also formated with c(model evaluation, replication, time, variable). +The output variable `out` containing all the input arguments detailed before and the calculated SI of first order (`mSI`), interaction (`iSI`), and total order (`tSI`). The convergence indices are also stored in the list named `mCI`, `iCI`, and `tCI`. The output are formated as 4-D array in `y` with the dimension name of model evaluation, number of replications, number of time points, and number of output variables, respectively. ```{r, eval=F} dim(out$y) ``` -The `check()` is a useful function to determine which parameters have relative lower sensitivity measurement across the given time interval, and therefore can be applied parameter fixing in model calibration. The argument of `SI.cutoff` is setting at 0.5 to detect the relative non-influential parameters in this case. +Some functions in **pksensi** provide a efficient way to check the result from global SA. The `check` can determine which parameters have relatively lower sensitivity measurement across the given time points and model outputs, and therefore can be applied parameter fixing in model calibration. The `check` also provides some feasible argument to specify the target output or change the cut-off value. The argument of `SI.cutoff` is setting at 0.05 to detect the relative non-influential parameters as default. ```{r, eval=F} -check(out, SI.cutoff = 0.5) +check(out, SI.cutoff = 0.05) ``` -Based on the sensitivity measurement of the total order, the result shows that `kgutabs` has relative lower measurement of sensitivity index. +Based on the sensitivity measurement of the total order, the result shows that `BW` has a relative lower measurement of SI. However, all parameters do not converge to the setting cut-off, which means the larger sample size is required in further sensitivity testing. Same as `plot` function that can assign specific output variable in the examination, the `check` function can also use the assignment (`vars`) to exam the given output. -## Model implementations with GNU MCSim +### Model implementation with GNU MCSim -Alternatively, to solve ODE by using GNU MCSim, we need to change the argument to `application = mcsim` in `compile_model()`. Rather than apply R `deSolve` to solve differential equations, the GNU MCSim can provide higher computational speed in global sensitivity analysis. +In addtion to use **deSolve** to solve differential equations in PK model, the **GNU MCSim** can provide better computational efficiency. To solve ODE through **GNU MCSim**, we need to change the argument to `application = mcsim` in `compile_model`. The computing time of using `solve_fun` in SA is estimated as, ```{r, eval=F} -system.time(y<-solve_fun(x, times, initState = initState, outnames = Outputs, dllname = mName)) +system.time(out<-solve_fun(x, t, initState = initState, + outnames = Outputs, dllname = mName)) ``` -```{r, eval=F} +Then, before we conduct the SA through **GNU MCSim**, The following code is used to compile the **GNU MCSim** model code to the exexutable program. + +```{r, eval=F, message=F, warning=F} compile_model(mName, application = "mcsim") ``` -Sililiar to `solve_fun` that can define the initial parameter and state values through input function, the `solve_mcsim` has a `condition` argument that is used to givien the specific input value such as oral dose or fixing parameter value or initial state variable. +Similar to `solve_fun` function that can define the initial value of parameter and state variable through generated functions, the `solve_mcsim` also has a `condition` argument that can be used to give the specific input value such as exposure dose, fixed parameter value or initial condition of state variable. ```{r, eval=F} -conditions <- c("Agutlument = 10") # Set the initial state of Agutlument = 10 -system.time(y<-solve_mcsim(x, mName = mName, - params = params, - vars = Outputs, - time = times, - condition = conditions)) +conditions <- c("Agutlument = 10") +system.time(out <- solve_mcsim(x, mName = mName, params = params, + vars = Outputs, time = t, + condition = conditions)) ``` -Under the same given condition, it takes 5-6 (`deSolve`) and 1-2 (GNU MCSim) seconds to solve model. The `solve_mcsim()` shows the better computational performance than `solve_fun()` in `pksensi`. +After solving the equations under the same given condition, we can find that **GNU MCSim** has about 4 - 5 times faster in computing performance than using **deSolve**. In this case, we only focus on performing the global SA alone for generic PK model without additonal comparison with experiment data. The next example will display and reproduce our previous published result [@fphar201800588] with full global SA workflow. ## References