In this example, We use a simple, one-compartment PK model from `httk`

package (Pearce et al. 2017) to demonstrate how `pksensi`

can be applied to pharmacokinetic studies. The differential equations for the one-compartment pharmacokinetic model can be written as:

\[\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).

In the beginning, we need to pre-install GNU MCSim (Bois and Maszle 1997) 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.

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.

The following R script can download and compile the model description file (`pbtk1cpt.model`

) and use `deSolve`

package (Soetaert, Petzoldt, and Setzer 2010) to solve ordinary differential equations in our model. The example model code of one-compartment PBTK model is available with pksensi package:

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 parameter values and initial states can be customized to specify the properties and schedule for the given dosing scenario.

```
parms <- initParms()
parms["vdist"] <- 0.74
parms["ke"] <- 0.28
parms["kgutabs"] <- 2.18
initState <- initStates(parms=parms)
initState["Agutlument"] <- 10
```

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`

):

```
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)
```

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:

```
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:

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.

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):

```
par(mfrow = c(3,3), mar = c(2,2,2,2), oma = c(2,2,1,1))
plot(x$a[,1,"vdist"], y[,1,"0.01",], main = "vdist")
text(1, .7, "t=0.01",cex = 1.2)
plot(x$a[,1,"ke"], y[,1,"0.01",], main = "ke")
plot(x$a[,1,"kgutabs"], y[,1,"0.01",], main = "kgutabs")
plot(x$a[,1,"vdist"], y[,1,"2.01",])
text(1, 18, "t=2.01",cex = 1.2)
plot(x$a[,1,"ke"], y[,1,"2.01",])
plot(x$a[,1,"kgutabs"], y[,1,"2.01",])
plot(x$a[,1,"vdist"], y[,1,"24.01",])
text(1, .7, "t=24.01",cex = 1.2)
plot(x$a[,1,"ke"], y[,1,"24.01",])
plot(x$a[,1,"kgutabs"], y[,1,"24.01",])
mtext("parameter", SOUTH<-1, line=0.4, outer=TRUE)
mtext("Ccompartment", WEST<-2, line=0.4, outer=TRUE)
```

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

In addition, the output are also formated with c(model evaluation, replication, time, variable).

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.

Based on the sensitivity measurement of the total order, the result shows that `kgutabs`

has relative lower measurement of sensitivity index.

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.

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.

```
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))
```

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`

.

Bois, Frederic, and Don Maszle. 1997. “MCSim: A Monte Carlo Simulation Program.” *Journal of Statistical Software, Articles* 2 (9): 1–60. https://doi.org/10.18637/jss.v002.i09.

Pearce, Robert, R. Woodrow Setzer, Cory Strope, Nisha Sipes, and John Wambaugh. 2017. “httk: R Package for High-Throughput Toxicokinetics.” *Journal of Statistical Software, Articles* 79 (4): 1–26. https://doi.org/10.18637/jss.v079.i04.

Soetaert, Karline, Thomas Petzoldt, and R. Woodrow Setzer. 2010. “Solving Differential Equations in R: Package deSolve.” *Journal of Statistical Software, Articles* 33 (9): 1–25. https://doi.org/10.18637/jss.v033.i09.