Exercise 2.3: Monte Carlo Simulation for PBPK model
Nan-Hung Hsieh
2019/12/09 (update: 2019-12-07)
TASK
- Reproduce the published Monte Carlo analysis result in Bois and Brochot (2016)
- Examine the butadiene concentration in fat after the exposure of 1 day.
- Testing parameter include body mass (
BDM), pulmonary flow (Flow_pul), partition coefficient of arterial blood (PC_art) and metabolism rate (Kmetwp)
- Construct the relationship between body mass and quantity in fat
SOLUTION
Package
library(deSolve)Inputs
C_inh <- approxfun(x = c(0, 120), y=c(10,0), method = "constant", f = 0, rule = 2) Parameters and outputs
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)
"Pct_Deadspace" = 0.7, # Fraction of pulmonary deadspace
"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 metabolismy <- 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)Model
# Define the model equations
bd.model = function(t, y, parameters) {
with (as.list(y), {
with (as.list(parameters), {
# Define some known constants
Height = 1.6 # use to calculate fraction of body fat
Age = 40 # use to calculate fraction of body fat
Sex = 1 # use to calculate fraction of body fat
MW_bu = 54.0914 # butadiene molecular weight (in grams)
# Conversions from/to ppm
ppm_per_mM = 24450 # ppm to mM under normal conditions
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 + Pct_Deadspace * C_inh.current
}
# 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.modelMonte Carlo simulation
Sample randonly some parameters
start_time <- Sys.time()
set.seed(1)
for (iteration in 1:1000){
parameters["BDM"] <- rnorm(1, 73, 7.3)
parameters["Flow_pul"] <- rnorm(1, 5, 0.5)
parameters["PC_art"] <- rnorm(1, 2, 0.2)
parameters["Kmetwp"] <- rnorm(1, 0.25, 0.025)
# We focus on time at 24 hour
times <- c(0, 1440)
# 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["PC_art"],
parameters["Kmetwp"])
} else { # accumulate
results <- rbind(results, tmp[2,-1])
sampled.parms <- rbind(sampled.parms,
c(parameters["BDM"],
parameters["Flow_pul"],
parameters["PC_art"],
parameters["Kmetwp"]))
}
}
end_time <- Sys.time()Output analysis
end_time - start_time## Time difference of 30.19761 secshead(sampled.parms)## BDM Flow_pul PC_art Kmetwp
## sampled.parms 68.42689 5.091822 1.832874 0.2898820
## 75.40541 4.589766 2.097486 0.2684581
## 77.20320 4.847306 2.302356 0.2597461
## 68.46494 3.892650 2.224986 0.2488767
## 72.88181 5.471918 2.164244 0.2648475
## 79.70853 5.391068 2.014913 0.2002662head(results)## Q_fat Q_wp Q_pp Q_met C_ven
## results 0.2284525 1.035820e-04 0.001486063 1.541566 8.649140e-05
## 0.2487108 9.556971e-05 0.001439920 1.499357 7.662418e-05
## 0.2774336 1.144274e-04 0.001620684 1.644112 8.284359e-05
## 0.2131533 9.661675e-05 0.001516974 1.318391 8.228833e-05
## 0.2808005 1.388960e-04 0.001795947 1.775749 9.580122e-05
## 0.2951133 1.440454e-04 0.001571253 1.611567 8.321394e-05
## C_art
## results 5.332478e-05
## 4.964288e-05
## 5.540840e-05
## 5.441045e-05
## 6.274876e-05
## 5.314532e-05BDM <- sampled.parms[,1]
hist(BDM)
Q_fat <- results[,1]
hist(Q_fat)
plot(BDM, Q_fat, xlab="Body mass (kg)", ylab="Quantity in Fat (mg)")
Reference
Bois F.Y., Brochot C. (2016) Modeling Pharmacokinetics. In: Benfenati E. (eds) In Silico Methods for Predicting Drug Toxicity. Methods in Molecular Biology, vol 1425. Humana Press, New York, NY