Exercise 1.4 Developing the PBPK model

TASK

• Instead of simple PK model, in this exercise we want to apply the well developed PBPK model for 1,3 butadiene.
• First, reproduce the simulation result from the published paper*.
• Exposure: Continuous exposure of the 10 ppm 1,3-butadiene for 2 hour.

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

SOLUTION

Packages

library(deSolve)

Model

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

      # Calculate flow and volumes
      
      # Calculate the tissue, blood, and air

      # Differentials for quantities

      # 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

Constants

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

Air and blood flow

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

Volumes

# 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

Concentrations - tissues & blood

# 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

Concentrations - air

# 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
}

Differentials for quantities

# 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

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

Inputs

C_inh <- approxfun(x = c(0, 120), y=c(10,0), method = "constant", f = 0, rule = 2) 

plot(C_inh(1:300), type="l", xlab = "Time (min)", ylab = "Butadiene inhaled concentration (ppm)")

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 metabolism

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)

Simulation

# Define the computation output times
t <- seq(from=0, to=1440, by=10)
t

##   [1]    0   10   20   30   40   50   60   70   80   90  100  110  120  130
##  [15]  140  150  160  170  180  190  200  210  220  230  240  250  260  270
##  [29]  280  290  300  310  320  330  340  350  360  370  380  390  400  410
##  [43]  420  430  440  450  460  470  480  490  500  510  520  530  540  550
##  [57]  560  570  580  590  600  610  620  630  640  650  660  670  680  690
##  [71]  700  710  720  730  740  750  760  770  780  790  800  810  820  830
##  [85]  840  850  860  870  880  890  900  910  920  930  940  950  960  970
##  [99]  980  990 1000 1010 1020 1030 1040 1050 1060 1070 1080 1090 1100 1110
## [113] 1120 1130 1140 1150 1160 1170 1180 1190 1200 1210 1220 1230 1240 1250
## [127] 1260 1270 1280 1290 1300 1310 1320 1330 1340 1350 1360 1370 1380 1390
## [141] 1400 1410 1420 1430 1440

# Solve ODE
out <- ode(times=t, func=bd.model, y=y, parms=parameters)
head(out)

##      time      Q_fat       Q_wp       Q_pp      Q_met       C_ven
## [1,]    0 0.00000000 0.00000000 0.00000000 0.00000000 0.000000000
## [2,]   10 0.02293618 0.03724892 0.07427798 0.06645654 0.003338318
## [3,]   20 0.04722954 0.04026245 0.14189019 0.16431358 0.004379098
## [4,]   30 0.07221315 0.04152080 0.20176415 0.26661643 0.005210310
## [5,]   40 0.09777386 0.04256838 0.25471138 0.37175647 0.005941936
## [6,]   50 0.12383371 0.04349410 0.30153555 0.47935418 0.006589697
##           C_art
## [1,] 0.01606403
## [2,] 0.01819035
## [3,] 0.01885327
## [4,] 0.01938270
## [5,] 0.01984870
## [6,] 0.02026129

Plot

plot(out)