One-compartment PBTK model

Equations

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

Model implementations with R deSolve package

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:

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

mName <- "pbtk1cpt"
compile_model(mName, application = "R")
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.

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

parms

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

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)

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:

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:

y <- solve_fun(x, times, initState = initState, outnames = Outputs, dllname = mName)
tell2(x,y)

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.

plot(x)

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

dim(x$a)

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

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

check(x, SI.cutoff = 0.5)

Based on the sensitivity measurement of the total order, the result shows that kgutabs has relative lower measurement of sensitivity index.

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

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

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.

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.