# Sensitivity Analysis for PBPK Model

The post was updated on 2019-12-31 under Windows 10 x64.

## Intro

In my previous post, I used python SALib to conduct the Sobol sensitivity analysis for the simple example (Ishigami function). Now, I want to take the more complicated case to understand how to apply sensitivity analysis in a dynamic model.

There is an excellent example of the dynamic model, the physiological-based pharmacokinetic (PBPK) model provided from the book In silico Methods for Predicting Drug Toxicity (chapter 3). This model was used to simulate the time-course of the concentration for the 1, 3-Butadiene in human body fluids and organs. The authors also provided the R source code to conduct the reproducible PBPK modeling.

It is an opportunity to apply the sensitivity analysis for this model. Therefore, I choose the elementary effect-based Morris screening method (Morris, 1991), which is also the advance local sensitivity analysis that change one-factor-at-a-time (OAT) to investigate the parameter impact on the model output.

The Morris method is an extension of the traditional OAT process that provide global converge over multiple points in the multi-dimensional parameter spaces. For each parameter, the Morris approach can calculate the mean OAT sensitivities $$\mu$$ and its standard deviation $$\sigma$$, which represent the overall influence and the interacts/non-linear effects on the model output, respectively.

An improved elementary effect method was developed by Campolongo et al. (2007). They improved both the reliability of the sensitivity measures (adding the mean absolute value of the OAT sensitives, denoted $$\mu^*$$) as well as the efficiency of the sampling strategy (generating a large number of proposed sampling trajectories and selecting a subset with the highest diversion across parameter space).

Compare with variance-based sensitivity analysis, the Morris method has high computing stable (lower computational cost) and only need a few sample number (lower computational resource) to reach the convergence. However, there are some disadvantages of the local method, e.g., the percentage of the output variations cannot be quantified.

Due to this PBPK case study is running in the pure R with the deSolve package, it might take a longer time to solve the ODEs. To have better computational performance, the Morris method is the best candidate.

# PBPK modeling

First, we need to define and initialize the state variables. The used PBPK model was constructed from three tissue compartments (fat, well-perfused, and poorly-perfused tissues).

y <- c("Q_fat" = 0,  # Quantity of butadiene in fat (mg)
"Q_wp" = 0,   # ~ in well-perfused (mg)
"Q_pp" = 0,   # ~ in poorly-perfused (mg)
"Q_met" = 0)  # ~ metabolized (mg)

## Define the model parameters

In this PBPK model, the unit of volumes is liter, the unit of time is minute, and the unit of flows is liter per minute.

parameters <- c("BDM" = 73,    # Body mass (kg)
"Height" = 1.6, # Body height (m)
"Age" = 40,     # in years
"Sex" = 1,      # code 1 is male, 2 is female
"Flow_pul" = 5, # Pulmonary ventilation rate (L/min)
"Vent_Perf" = 1.14,    # Ventilation over perfusion ratio
"Pct_LBDM_wp" = 0.2,   # wp tissue as fraction of lean mass
"Pct_Flow_fat" = 0.1,  # Fraction of cardiac output to fat
"Pct_Flow_pp" = 0.35,  # ~ to pp
"PC_art" = 2,          # Blood/air partition coefficient
"PC_fat" = 22,         # Fat/blood ~
"PC_wp" = 0.8,         # wp/blood ~
"PC_pp" = 0.8,         # pp/blood ~
"Kmetwp" = 0.25)       # Rate constant for metabolism

## Define the exposure function.

The exposure is set to the periodic exposure (8-hour exposure per day) in a week.

min_per_week <- 60 * 24 * 7

x1 <- seq(1, min_per_week, 1440) # start exposure
x2 <- seq(480, min_per_week, 1440) # stop exposure
duration <- as.vector(matrix(c(x1, x2), nrow = 2, byrow = TRUE))
exposure <- rep(c(10,0.1),length(duration)/2)
C_inh <- approxfun(x = duration, y = exposure, method="constant", f=0, rule=2)

# Check the input concentration profile just defined
plot(C_inh(1:min_per_week), xlab = "Time (min)",
ylab = "Butadiene air concentration (ppm)", type = "l")

## Define the PBPK model

# Define the model equations
bd.model = function(t, y, parameters) {
with (as.list(y), {
with (as.list(parameters), {

# Define the fixed parameters
Height = 1.6
Age = 40
Sex = 1

# Define some useful constants
MW_bu = 54.0914 # butadiene molecular weight (in grams)
ppm_per_mM = 24450 # ppm to mM under normal conditions
# Conversions from/to ppm
ppm_per_mg_per_l = ppm_per_mM / MW_bu
mg_per_l_per_ppm = 1 / ppm_per_mg_per_l
# Calculate Flow_alv from total pulmonary flow
Flow_alv = Flow_pul * (1 - Pct_Deadspace)
# Calculate total blood flow from Flow_alv and the V/P ratio
Flow_tot = Flow_alv / Vent_Perf
# Calculate fraction of body fat
Pct_BDM_fat = (1.2 * BDM / (Height * Height) - 10.8
*(2 - Sex) +
0.23 * Age - 5.4) * 0.01
# Actual volumes, 10% of body mass (bones…) get no butadiene
Eff_V_fat = Pct_BDM_fat * BDM
Eff_V_wp = Pct_LBDM_wp * BDM * (1 - Pct_BDM_fat)
Eff_V_pp = 0.9 * BDM - Eff_V_fat - Eff_V_wp
# Calculate actual blood flows from total flow and percent flows
Flow_fat = Pct_Flow_fat * Flow_tot
Flow_pp = Pct_Flow_pp * Flow_tot
Flow_wp = Flow_tot * (1 - Pct_Flow_pp - Pct_Flow_fat)
# Calculate the concentrations
C_fat = Q_fat / Eff_V_fat
C_wp = Q_wp / Eff_V_wp
C_pp = Q_pp / Eff_V_pp
# Venous blood concentrations at the organ exit
Cout_fat = C_fat / PC_fat
Cout_wp = C_wp / PC_wp
Cout_pp = C_pp / PC_pp
# Sum of Flow * Concentration for all compartments
dQ_ven = Flow_fat * Cout_fat + Flow_wp * Cout_wp + Flow_pp * Cout_pp
C_inh.current = C_inh(t) # to avoid calling C_inh() twice
# Arterial blood concentration
# Convert input given in ppm to mg/l to match other units
C_art = (Flow_alv * C_inh.current * mg_per_l_per_ppm + dQ_ven) / (Flow_tot + Flow_alv / PC_art)
# Venous blood concentration (mg/L)
C_ven = dQ_ven / Flow_tot
# Alveolar air concentration (mg/L)
C_alv = C_art / PC_art
# Exhaled air concentration (ppm!)
if (C_alv <= 0) {
C_exh = 10E-30 # avoid round off errors
} else {
C_exh = (1 - Pct_Deadspace) * C_alv * ppm_per_mg_per_l +
}
# Quantity metabolized in liver (included in well-perfused)
dQmet_wp = Kmetwp * Q_wp
# Differentials for quantities
dQ_fat = Flow_fat * (C_art - Cout_fat)
dQ_wp = Flow_wp * (C_art - Cout_wp) - dQmet_wp
dQ_pp = Flow_pp * (C_art - Cout_pp)
dQ_met = dQmet_wp
# The function bd.model must return at least the derivatives
list(c(dQ_fat, dQ_wp, dQ_pp, dQ_met), # derivatives
c("C_ven" = C_ven, "C_art" = C_art)) # extra outputs
}) # end with parameters
}) # end with y
} # end bd.model

Define the computation output times and call the ODE solver.

library(deSolve)
times <- seq(from = 0, to = min_per_week, by = 10) # A week
results <- ode(times = times, func = bd.model, y = y, parms = parameters)

Plot the model simulation results.

plot(results, col="red")

# Sensitivity analysis

Firstly, we use the Morris method which was provided from sensitivity package.

library(sensitivity)
## Registered S3 method overwritten by 'sensitivity':
##   method    from
##   print.src dplyr

## Define the testing parameters

Twelve parameters are used in this test.

factors <- c("BDM", "Flow_pul", "Pct_Deadspace", "Vent_Perf",
"Pct_LBDM_wp", "Pct_Flow_fat", "Pct_Flow_pp",
"PC_art", "PC_fat", "PC_wp", "PC_pp",
"Kmetwp")

## Parameter uncertainty

Then, we need to set up the uncertainty for each parameter. Here we set the 10% variation in the testing parameters.

LL <- 0.9 # 10% lower limit
UL <- 1.1 # 10% upper limit

# Define the lower and upper limits that will be the input to the morris function
binf <- c(parameters["BDM"]*LL,
parameters["Flow_pul"]*LL,
parameters["Vent_Perf"]*LL,
parameters["Pct_LBDM_wp"]*LL,
parameters["Pct_Flow_fat"]*LL,
parameters["Pct_Flow_pp"]*LL,
parameters["PC_art"]*LL,
parameters["PC_fat"]*LL,
parameters["PC_wp"]*LL,
parameters["PC_pp"]*LL,
parameters["Kmetwp"]*LL)
bsup <- c(parameters["BDM"]*UL,
parameters["Flow_pul"]*UL,
parameters["Vent_Perf"]*UL,
parameters["Pct_LBDM_wp"]*UL,
parameters["Pct_Flow_fat"]*UL,
parameters["Pct_Flow_pp"]*UL,
parameters["PC_art"]*UL,
parameters["PC_fat"]*UL,
parameters["PC_wp"]*UL,
parameters["PC_pp"]*UL,
parameters["Kmetwp"]*UL)

To check the convergence of the sensitivity index, we can set up the sequence of the sampling number. After that, we can apply the for-loop to do the morris sensitivity analysis and link with numerical analysis from deSolve::ode function. In this case, I examine the parameter sensitivity under the exposure time at 8-hour (stop point of the first exposure) and investigate the parameter sensitivity. Also, we only focused on the parameter effect on venous blood concentration.

sample <- seq(from = 20, to = 140, by = 20)

for (i in 1:length(sample)) {
set.seed(12345)
x <- morris(model = NULL, factors = factors, r = sample[i],
design = list(type = "oat", levels = 6, grid.jump = 3),
binf = binf, bsup = bsup, scale = TRUE)

for (iteration in 1:nrow(x$X)) { parameters["BDM"] = x$X[iteration,"BDM"]
parameters["Flow_pul"] = x$X[iteration,"Flow_pul"] parameters["Pct_Deadspace"] = x$X[iteration,"Pct_Deadspace"]
parameters["Vent_Perf"] = x$X[iteration,"Vent_Perf"] parameters["Pct_LBDM_wp"] = x$X[iteration,"Pct_LBDM_wp"]
parameters["Pct_Flow_fat"] = x$X[iteration,"Pct_Flow_fat"] parameters["Pct_Flow_pp"] = x$X[iteration,"Pct_Flow_pp"]
parameters["PC_art"] = x$X[iteration,"PC_art"] parameters["PC_fat"] = x$X[iteration,"PC_fat"]
parameters["PC_wp"] = x$X[iteration,"PC_wp"] parameters["PC_pp"] = x$X[iteration,"PC_pp"]
parameters["Kmetwp"] = x$X[iteration,"Kmetwp"] # We focus on time at 8 hour (stop exposure), times = c(0, 480) # Integrate tmp = ode(times = times, func = bd.model, y = y, parms = parameters) if (iteration == 1) { # initialize results = tmp[2,-1] sampled.parms = c(parameters["BDM"], parameters["Flow_pul"], parameters["Pct_Deadspace"], parameters["Vent_Perf"], parameters["Pct_LBDM_wp"], parameters["Pct_Flow_fat"], parameters["Pct_Flow_pp"], parameters["PC_art"], parameters["PC_fat"], parameters["PC_wp"], parameters["PC_pp"], parameters["Kmetwp"]) } else { # accumulate results = rbind(results, tmp[2,-1]) sampled.parms = rbind(sampled.parms, c(parameters["BDM"], parameters["Flow_pul"], parameters["Pct_Deadspace"], parameters["Vent_Perf"], parameters["Pct_LBDM_wp"], parameters["Pct_Flow_fat"], parameters["Pct_Flow_pp"], parameters["PC_art"], parameters["PC_fat"], parameters["PC_wp"], parameters["PC_pp"], parameters["Kmetwp"])) } } # tell(x, results[,"C_ven"]) # We focus on parameter effect on venous blood concentration if (i == 1){ X <- apply(abs(x$ee), 2, mean)
} else {
X <- rbind(X, apply(abs(x$ee), 2, mean)) } } ## Number of repetitions and parameter sensitivity Through this plot, we can examine the variation and converge of the parameter. It seems that each sensitivity index of the parameter was stable across all sampling number. Also, we can find the fraction of pulmonary dead space Pct_Deadspace is the most influential parameter. # Manipulate the output result for plotting row.names(X) <- sample meltX <- reshape::melt(X) ee_lim <- c(min(abs(x$ee)), max(abs(x$ee))) library(viridis) # use viridis to create distinct colors  ## Loading required package: viridisLite par(mar=c(4,4,2,1)) plot(sample, subset(meltX, X2 == factors[1])[,3], type="b", col=viridis_pal()(12)[1], frame.plot = FALSE, xlab="number of repetitions", ylab=expression(paste(mu,"*")), ylim=ee_lim, xlim = c(20,160)) for( i in 2:12){ lines(sample, subset(meltX, X2 == factors[i])[,3], type="b", col=viridis_pal()(12)[i]) } text(140, subset(meltX, X1 == 140)[,3], factors, pos = 4) ## Mean OAT sensitivities and standard deviation This visualization provides the information to examine the parameter sensitivity and the relationships between model parameters and its related output. According to the output result, we can find the relationship between the parameter values and model outputs are between linear and monotonic. plot(x, xlim=ee_lim, main ="number of repetitions = 140") abline(0,1) # non-linear and/or non-monotonic abline(0,0.5, lty = 2) # monotonic abline(0,0.1, lty = 3) # almost linear legend("bottomright", legend = c("non-linear and/or non-monotonic", "monotonic", "linear"), lty = c(1:3)) ## Parameter variation and model output To investigate the parameter effect on the PBPK model outputs, we used the following visualization methods for all parameters and the most influential parameter Pct_Deadspace and less influential parameter PC_pp. For all parameters par(mfrow=c(3,4)) for (i in 1:12){ plot(x$X[,i], x$y, main = x$factors[i], xlab = "", ylab = "predict conc.")
}

For most influential parameter Pct_Deadspace

library(ggplot2)
library(plyr)

sensiX <- data.frame(x$X[,"Pct_Deadspace"], x$y)
musensi <- ddply(sensiX, "Pct_Deadspace", summarize, mean.conc=mean(Predict.conc))

linetype="dashed") +
theme_minimal()

For less influential parameter PC_pp

insensiX <- data.frame(x$X[,"PC_pp"], x$y)
names(insensiX) <- c("PC_pp", "Predict.conc")
muinsensi <- ddply(insensiX, "PC_pp", summarize, mean.conc=mean(Predict.conc))

ggplot(insensiX, aes(x=Predict.conc, color=as.factor(PC_pp))) +
geom_density() + scale_color_discrete(name="PC_pp") +
geom_vline(data=muinsensi, aes(xintercept=mean.conc, color=as.factor(PC_pp)),
linetype="dashed") +
theme_minimal()

Through this method, we can quickly check the parameter sensitivity in the PBPK model and compare the parameter effects on model output. However, this is just a simple test case. To do the “complete” sensitivity analysis, we still need to find the reasonable parameter uncertainty and its range (the sensitivity may change if we revise the range of parameter uncertainty). However, this step-by-step method can provide essential information for beginners to apply the sensitivity in the dynamic model. Also, the parameter sensitivity might change with time. It would be better to investigate the time-varying sensitivity index in the factor prioritization.

# Session info

sessionInfo()
## R version 3.6.2 (2019-12-12)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17763)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base
##
## other attached packages:
## [1] plyr_1.8.5         ggplot2_3.2.1      viridis_0.5.1      viridisLite_0.3.0
## [5] sensitivity_1.16.3 deSolve_1.25
##
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.3        compiler_3.6.2    pillar_1.4.2      tools_3.6.2
##  [5] boot_1.3-23       digest_0.6.23     evaluate_0.14     lifecycle_0.1.0
##  [9] tibble_2.1.3      gtable_0.3.0      pkgconfig_2.0.3   rlang_0.4.2
## [13] cli_2.0.0         rstudioapi_0.10   yaml_2.2.0        blogdown_0.17
## [17] xfun_0.11         gridExtra_2.3     dplyr_0.8.3       withr_2.1.2
## [21] stringr_1.4.0     knitr_1.26        tidyselect_0.2.5  grid_3.6.2
## [25] reshape_0.8.8     glue_1.3.1        R6_2.4.1          fansi_0.4.0
## [29] rmarkdown_2.0     bookdown_0.16     sessioninfo_1.1.1 farver_2.0.1
## [33] purrr_0.3.3       magrittr_1.5      scales_1.1.0      htmltools_0.4.0
## [37] assertthat_0.2.1  colorspace_1.4-1  labeling_0.3      stringi_1.4.3
## [41] lazyeval_0.2.2    munsell_0.5.0     crayon_1.3.4