TASK

Based on the OSHA permissible exposure limit (PEL), the time-weighted average is 1 ppm. If a worker is under the exposure for a long time, what is the estimated blood concentration.

SOLUTION

Packages

library(deSolve)

Parameters, Variables, and Model

Parameters

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

Variable

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)

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

Exposure setting

t_exposure <- 43200 # 30 days
C_inh <- approxfun(x = c(0, t_exposure), y=c(10,0), method = "constant", f = 0, rule = 2) 
plot(C_inh(1:t_exposure), type="l", xlab = "Time (min)", ylab = "Butadiene inhaled concentration (ppm)")

Simulation

t <- seq(from=0, to=t_exposure, by=100)
out <- ode(times=t, func=bd.model, y=y, parms=parameters)
tail(out, 10)
##         time    Q_fat       Q_wp      Q_pp    Q_met      C_ven       C_art
## [424,] 42300 11.24631 0.05555652 0.7245784 582.4532 0.01517190 0.025727662
## [425,] 42400 11.24632 0.05555653 0.7245785 583.8421 0.01517190 0.025727663
## [426,] 42500 11.24632 0.05555653 0.7245785 585.2310 0.01517190 0.025727664
## [427,] 42600 11.24632 0.05555653 0.7245785 586.6199 0.01517190 0.025727664
## [428,] 42700 11.24633 0.05555653 0.7245785 588.0088 0.01517190 0.025727665
## [429,] 42800 11.24633 0.05555653 0.7245785 589.3977 0.01517190 0.025727665
## [430,] 42900 11.24633 0.05555653 0.7245785 590.7866 0.01517191 0.025727666
## [431,] 43000 11.24633 0.05555653 0.7245785 592.1756 0.01517191 0.025727666
## [432,] 43100 11.24634 0.05555653 0.7245786 593.5645 0.01517191 0.025727667
## [433,] 43200 11.24634 0.05555619 0.7245784 594.9534 0.01517188 0.009663619
plot(out)