One-compartment PK model

Nan-Hung Hsieh



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

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\]

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

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 (Soetaert, Petzoldt, and Setzer 2010) package or with GNU MCSim (Bois and Maszle 1997) 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 will download the one-compartment PK model from 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\)).

cat(readLines("pbtk1cpt.model"), sep = "\n")

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 on other systems) and R file (pbtk1cpt_inits.R) with default input parameters and initial state settings with the definition of application = "R".

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.

source(paste0(mName, "_inits.R"))

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 package, which includes physico-chemical properties for over 500 chemicals. The given value of vdist, ke, and kgutabs in are 1.1 (L/kg BW), 0.23 (/h), and 2.18 (/h), respectively. The body weight is assumed to be 70 (kg).

parms <- initParms()
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.


Using the ode function in deSolve package, we can visualize the pharmacokinetic profile according to the given parameter baseline and the time points (t):

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)

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.

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.

x <- rfast99(params, n = 200, q = q, q.arg = q.arg, replicate = 20)

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.


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.

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.

outputs <- c("Ccompartment", "Ametabolized")
out <- solve_fun(x, time = t, initState = initState, 
                 outnames = outputs, dllname = mName)

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.


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.

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.

plot(out, vars = "Ametabolized")

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.

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.

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)

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.


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.

check(out, SI.cutoff = 0.05)

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 implementation with GNU MCSim

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,

system.time(out<-solve_fun(x, t, initState = initState, 
                           outnames = Outputs, dllname = mName))

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.

compile_model(mName, application = "mcsim")

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.

conditions <- c("Agutlument = 10") 
system.time(out <- solve_mcsim(x, mName = mName, params = params, 
                               vars = Outputs, time = t, 
                               condition = conditions))

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 (Hsieh et al. 2018) with full global SA workflow.


Bois, Frederic, and Don Maszle. 1997. “MCSim: A Monte Carlo Simulation Program.” Journal of Statistical Software, Articles 2 (9): 1–60.

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.

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.

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.