Skip to contents

Bisphenol A / S / F / AF PBPK (Karrer 2018)

Source: vignettes/articles/Karrer_2018_bisphenols_pbpk.Rmd
Karrer_2018_bisphenols_pbpk.Rmd

Model and source

  • Citation: Karrer C, Roiss T, von Goetz N, Gramec Skledar D, Peterlin Masic L, Hungerbuhler K. (2018). Physiologically Based Pharmacokinetic (PBPK) Modeling of the Bisphenols BPA, BPS, BPF, and BPAF with New Experimental Metabolic Parameters: Comparing the Pharmacokinetic Behavior of BPA with Its Substitutes. Environmental Health Perspectives 126(7):077002.
  • Article: https://doi.org/10.1289/EHP2739 (PMCID PMC6108829)
  • ETH Zurich. Original implementation in R with the deSolve package; the complete model code and its input tables are deposited in the paper’s Supplemental Material.

Bisphenol A (BPA) is a high-production-volume plasticiser and a known endocrine disruptor; BPS, BPF, and BPAF are the structural analogs increasingly used to replace it. This paper builds a whole-body PBPK model for each bisphenol and its glucuronide and sulfate conjugates, after both peroral and dermal exposure, to compare the internal exposures that would result from replacing BPA with each analog. The BPA model is a recalibration of Yang et al. (2015) against the Thayer et al. (2015) biomonitoring data (the authors found and fixed a body-weight scaling error in the published gut-glucuronidation vmax); BPS was calibrated against Oh et al. (2018); BPF and BPAF could not be calibrated (no in-vivo data exist) and were parameterised by analogy with new in-vitro glucuronidation kinetics measured in human microsomes.

Four models are packaged, one per bisphenol, matching the four the authors built:

mod_bpa <- nlmixr2lib::readModelDb("Karrer_2018_bisphenolA_pbpk")
mod_bps <- nlmixr2lib::readModelDb("Karrer_2018_bisphenolS_pbpk")
mod_bpf <- nlmixr2lib::readModelDb("Karrer_2018_bisphenolF_pbpk")
mod_bpaf <- nlmixr2lib::readModelDb("Karrer_2018_bisphenolAF_pbpk")

They share one ODE structure and the chemical-independent physiology; they differ in molecular weight, tissue:serum partition coefficients, glucuronidation kinetics (and its functional form), the fraction of glucuronide subject to enterohepatic recirculation (EHR), and the oral-uptake / urinary-clearance constants inherited from the calibrated BPA or BPS model. All four are deterministic typical-value simulators: the paper reports single point estimates with no between-subject variance and no residual-error model, so none of the files carries eta* or residual-error terms. Inventing variances would be fabrication; see Assumptions, deviations and errata.

Population and biological context 

data.frame(
  Field = c("Species", "Reference individual", "Routes", "Doses compared", "Focus"),
  Value = c(
    "Human (physiology Edginton et al. 2006 / ICRP 2002)",
    "Adult female, 60 kg, 30 years (Table 4); structure also covers infants, toddlers, children, adolescents and adult males via Table S1",
    "Peroral (gastric emptying then gut + hepatic first pass) and dermal (two parallel skin-surface depots: thermal paper and personal care products)",
    "Single 500 ng/kg BW peroral or dermal (Fig 4); age/route-specific environmental estimates (Fig 5)",
    "Females of childbearing age; serum and gonadal concentrations of the unconjugated bisphenols"
  )
) |> knitr::kable()
Field Value
Species Human (physiology Edginton et al. 2006 / ICRP 2002)
Reference individual Adult female, 60 kg, 30 years (Table 4); structure also covers infants, toddlers, children, adolescents and adult males via Table S1
Routes Peroral (gastric emptying then gut + hepatic first pass) and dermal (two parallel skin-surface depots: thermal paper and personal care products)
Doses compared Single 500 ng/kg BW peroral or dermal (Fig 4); age/route-specific environmental estimates (Fig 5)
Focus Females of childbearing age; serum and gonadal concentrations of the unconjugated bisphenols

The reference physiology (cardiac output qcc, the fv_* tissue-volume fractions and the fq_* blood-flow fractions) is held at the adult-female values. Cardiac output is an absolute per-age-group quantity and is not body-weight scaled; the tissue volumes, the saturable-metabolism vmax terms, the first-order rate constants and the urinary clearances all scale with body weight WT (linearly, by liver mass, by bw^-0.25, or by bw^0.75 as the deposited code dictates). To simulate another age group, override WT, qcc and the fv_* / fq_* fractions with that group’s row from Table S1.

Model structure

Twenty-one ODE states per bisphenol: eleven for the parent, four each for the glucuronide and sulfate, and two dermal skin-surface depots.

data.frame(
  Layer = c(
    rep("Parent", 11), rep("Glucuronide", 4), rep("Sulfate", 4), rep("Dermal depots", 2)
  ),
  State = c(
    "stomach", "a_small_intestine", "a_liver", "a_fat", "a_gonads", "a_skin",
    "a_brain", "a_rapidly_perfused", "a_slowly_perfused", "plasma", "urine",
    "a_small_intestine_gluc", "gut_lumen_gluc", "serum_gluc", "urine_gluc",
    "a_small_intestine_sulf", "gut_lumen_sulf", "serum_sulf", "urine_sulf",
    "depot_td_thermalpaper", "depot_td_personalcare"
  ),
  Role = c(
    "Oral dose site; empties to the small intestine at rate ge",
    "Enterocyte lumen; gut glucuronidation vs absorption to liver compete",
    "Liver; gut-absorbed parent + first-pass glucuronidation and sulfation",
    "Fat (flow-limited)", "Gonads (flow-limited; endocrine-disruption endpoint)",
    "Skin (flow-limited; receives dermal penetration)", "Brain (flow-limited)",
    "Richly perfused lump (flow-limited)", "Slowly perfused lump (flow-limited)",
    "Well-mixed serum pool (arterial supply + venous return - renal loss)",
    "Cumulative urinary parent",
    "Glucuronide formed in enterocytes", "Glucuronide in the gut EHR pool",
    "Glucuronide in serum submodule", "Cumulative urinary glucuronide",
    "Sulfate formed in enterocytes (inert: gut sulfation off)",
    "Sulfate in the gut EHR pool (inert: no sulfate EHR)",
    "Sulfate in serum submodule", "Cumulative urinary sulfate",
    "Skin-surface depot, thermal paper (20% absorbed, 6 h penetration half-life)",
    "Skin-surface depot, personal care products (60% absorbed, 10 min half-life)"
  )
) |> knitr::kable()
Layer State Role
Parent stomach Oral dose site; empties to the small intestine at rate ge
Parent a_small_intestine Enterocyte lumen; gut glucuronidation vs absorption to liver compete
Parent a_liver Liver; gut-absorbed parent + first-pass glucuronidation and sulfation
Parent a_fat Fat (flow-limited)
Parent a_gonads Gonads (flow-limited; endocrine-disruption endpoint)
Parent a_skin Skin (flow-limited; receives dermal penetration)
Parent a_brain Brain (flow-limited)
Parent a_rapidly_perfused Richly perfused lump (flow-limited)
Parent a_slowly_perfused Slowly perfused lump (flow-limited)
Parent plasma Well-mixed serum pool (arterial supply + venous return - renal loss)
Parent urine Cumulative urinary parent
Glucuronide a_small_intestine_gluc Glucuronide formed in enterocytes
Glucuronide gut_lumen_gluc Glucuronide in the gut EHR pool
Glucuronide serum_gluc Glucuronide in serum submodule
Glucuronide urine_gluc Cumulative urinary glucuronide
Sulfate a_small_intestine_sulf Sulfate formed in enterocytes (inert: gut sulfation off)
Sulfate gut_lumen_sulf Sulfate in the gut EHR pool (inert: no sulfate EHR)
Sulfate serum_sulf Sulfate in serum submodule
Sulfate urine_sulf Cumulative urinary sulfate
Dermal depots depot_td_thermalpaper Skin-surface depot, thermal paper (20% absorbed, 6 h penetration half-life)
Dermal depots depot_td_personalcare Skin-surface depot, personal care products (60% absorbed, 10 min half-life)

Source trace

Every parameter in each model file carries an in-file comment pointing to its source location. The chemical-specific parameters come from the following tables.

data.frame(
  Quantity = c(
    "Molecular weight (g/mol)",
    "Tissue:serum partition coefficients",
    "Hepatic glucuronidation Km / vmax / Ksi",
    "Enterocyte glucuronidation Km / vmax / Ksi",
    "Hepatic sulfation Km / vmax (all = BPA)",
    "Gastric emptying, enterocyte volume",
    "Oral uptake small intestine -> liver",
    "Urinary clearance (parent)",
    "EHR fraction (1 - met1g)",
    "Dermal extent of absorption / half-life / period",
    "Peroral uptake period (3 min)",
    "Reference physiology (qcc, fv_*, fq_*)"
  ),
  Source = c(
    "Discussion (glucuronide MW minus glucuronyl increment 176.12); BPA from deposited code",
    "Table 6 (Doerge et al. 2011 for BPA; Zhang and Zhang 2006 QSAR otherwise)",
    "Table 7 ('used in basic model'; Ksi only for BPAF liver)",
    "Table 7 ('used in basic model'; Ksi for BPS, BPF, BPAF gut)",
    "Table 7 (Kurebayashi et al. 2010; approximated for the analogs)",
    "Table 8 / deposited code",
    "Table 8 (BPA 2.1; BPS/BPAF 5.0 calibrated; BPF = BPA)",
    "Table 8 (BPA/BPF 0.06; BPS/BPAF 0.30)",
    "Table 8 (BPA 10%, BPS 67%, BPF 0%, BPAF 70%)",
    "Table 5",
    "deposited code (period.O = 3/60 h), not Table 5's 0.25 h",
    "Table S1 / Table C3 (adult female)"
  )
) |> knitr::kable()
Quantity Source
Molecular weight (g/mol) Discussion (glucuronide MW minus glucuronyl increment 176.12); BPA from deposited code
Tissue:serum partition coefficients Table 6 (Doerge et al. 2011 for BPA; Zhang and Zhang 2006 QSAR otherwise)
Hepatic glucuronidation Km / vmax / Ksi Table 7 (‘used in basic model’; Ksi only for BPAF liver)
Enterocyte glucuronidation Km / vmax / Ksi Table 7 (‘used in basic model’; Ksi for BPS, BPF, BPAF gut)
Hepatic sulfation Km / vmax (all = BPA) Table 7 (Kurebayashi et al. 2010; approximated for the analogs)
Gastric emptying, enterocyte volume Table 8 / deposited code
Oral uptake small intestine -> liver Table 8 (BPA 2.1; BPS/BPAF 5.0 calibrated; BPF = BPA)
Urinary clearance (parent) Table 8 (BPA/BPF 0.06; BPS/BPAF 0.30)
EHR fraction (1 - met1g) Table 8 (BPA 10%, BPS 67%, BPF 0%, BPAF 70%)
Dermal extent of absorption / half-life / period Table 5
Peroral uptake period (3 min) deposited code (period.O = 3/60 h), not Table 5’s 0.25 h
Reference physiology (qcc, fv_, fq_) Table S1 / Table C3 (adult female)

Validation

The Supplemental Material deposits the authors’ complete deSolve model code, so the strongest validation is that the packaged rxode2 model reproduces that executable source. The two headline predictions the paper states in the text are the serum Cmax of unconjugated BPAF in a 60 kg adult female after a single 500 ng/kg BW dose: 1.79 pM peroral and 7.66 pM dermal (Results, “Comparison of Model Predictions”; Figure 4A).

# Molecular weights (g/mol); see the source trace above.
mw <- c(BPA = 228.28, BPS = 249.88, BPF = 199.88, BPAF = 335.88)
mods <- list(BPA = mod_bpa, BPS = mod_bps, BPF = mod_bpf, BPAF = mod_bpaf)
WT <- 60 # adult female reference (Table 4)

# Single-dose scenario: amt (nmol) = dose (ng/kg) * WT / MW, dosed into the
# named compartment with a modeled zero-order duration (rate = -2). Peroral
# EoA = 100% (applied as the dose); dermal EoA enters via the model's f().
solve_single <- function(mod, bp, route, dose_ng_per_kg = 500,
                         tmax = 120, dt = 0.05) {
  dose_nmol <- dose_ng_per_kg * WT / mw[[bp]]
  cmt <- if (route == "oral") "stomach" else "depot_td_thermalpaper"
  ev <- rxode2::et(amt = dose_nmol, cmt = cmt, rate = -2) |>
    rxode2::et(seq(0, tmax, by = dt))
  rxode2::rxSolve(mod, ev, c(WT = WT), returnType = "data.frame")
}

BPAF single-dose serum Cmax (reproducing the deposited code) 

bpaf_oral <- solve_single(mod_bpaf, "BPAF", "oral", tmax = 48, dt = 0.02)
bpaf_derm <- solve_single(mod_bpaf, "BPAF", "dermal", tmax = 120, dt = 0.05)

bpaf_cmp <- data.frame(
  Route = c("Peroral", "Dermal (thermal paper)"),
  `Paper Fig 4 (pM)` = c(1.79, 7.66),
  `Deposited code (pM)` = c(1.566, 7.745),
  `This model (pM)` = c(max(bpaf_oral$Cc) * 1000, max(bpaf_derm$Cc) * 1000),
  check.names = FALSE
)
bpaf_cmp |>
  dplyr::mutate(dplyr::across(dplyr::where(is.numeric), \(x) round(x, 3))) |>
  knitr::kable()
Route Paper Fig 4 (pM) Deposited code (pM) This model (pM)
Peroral 1.79 1.566 1.566
Dermal (thermal paper) 7.66 7.745 7.745

The packaged rxode2 model reproduces the deposited deSolve code to better than 0.1% on both routes. The dermal prediction is within 1.1% of the value the paper prints in Figure 4 (7.745 vs 7.66 pM); the peroral prediction (1.566 pM) is 12% below the printed 1.79 pM. That gap is a paper-versus-deposited-code discrepancy, not a transcription error: running the authors’ own deposited single-dose code with the printed Table 6-8 parameters gives 1.566 pM, and the rxode2 model matches it. The paper itself notes that its deposited input tables carry roundings; see Assumptions, deviations and errata. The value is not tuned to the printed figure (tuning a parameter to hit a target is prohibited).

stopifnot(
  # rxode2 reproduces the deposited-code BPAF serum Cmax to <0.5%.
  abs(max(bpaf_oral$Cc) * 1000 - 1.566) < 0.02,
  abs(max(bpaf_derm$Cc) * 1000 - 7.745) < 0.05,
  # BPAF is absorbed with heavy gut + hepatic first pass, so peroral serum is
  # far lower than dermal -- the paper's qualitative finding for BPAF.
  max(bpaf_oral$Cc) < max(bpaf_derm$Cc)
)

Cross-compound ordering (Figure 4A)

The paper’s central comparative finding is that, across both routes, BPS gives the highest serum Cmax and BPAF the lowest, with BPA and BPF similar. The packaged models reproduce that ordering.

cc <- lapply(names(mods), function(bp) {
  so <- solve_single(mods[[bp]], bp, "oral", tmax = 48, dt = 0.02)
  sd <- solve_single(mods[[bp]], bp, "dermal", tmax = 120, dt = 0.05)
  data.frame(
    Bisphenol = bp,
    `Peroral Cmax (pM)` = max(so$Cc) * 1000,
    `Dermal Cmax (pM)` = max(sd$Cc) * 1000,
    check.names = FALSE
  )
}) |> dplyr::bind_rows()
cc |>
  dplyr::mutate(dplyr::across(dplyr::where(is.numeric), \(x) signif(x, 3))) |>
  knitr::kable()
Bisphenol Peroral Cmax (pM) Dermal Cmax (pM)
BPA 29.00 14.60
BPS 569.00 28.00
BPF 35.40 17.20
BPAF 1.57 7.75 
oral <- setNames(cc$`Peroral Cmax (pM)`, cc$Bisphenol)
derm <- setNames(cc$`Dermal Cmax (pM)`, cc$Bisphenol)
stopifnot(
  # BPS highest, BPAF lowest, for BOTH routes (Fig 4A / Results text).
  which.max(oral) == which(names(oral) == "BPS"),
  which.min(oral) == which(names(oral) == "BPAF"),
  which.max(derm) == which(names(derm) == "BPS"),
  which.min(derm) == which(names(derm) == "BPAF"),
  # BPA and BPF are "rather similar" (within a factor of 2).
  abs(log(oral[["BPF"]] / oral[["BPA"]])) < log(2)
)

Serum concentration-time profiles 

prof <- lapply(names(mods), function(bp) {
  so <- solve_single(mods[[bp]], bp, "oral", tmax = 96, dt = 0.1)
  sd <- solve_single(mods[[bp]], bp, "dermal", tmax = 96, dt = 0.1)
  dplyr::bind_rows(
    data.frame(Bisphenol = bp, Route = "Peroral", time = so$time, Cc_pM = so$Cc * 1000),
    data.frame(Bisphenol = bp, Route = "Dermal", time = sd$time, Cc_pM = sd$Cc * 1000)
  )
}) |> dplyr::bind_rows()

ggplot(dplyr::filter(prof, time > 0, Cc_pM > 0), aes(time, Cc_pM, colour = Bisphenol)) +
  geom_line() +
  facet_wrap(~Route, scales = "free_y") +
  scale_y_log10() +
  labs(x = "Time (h)", y = "Serum concentration (pM)") +
  theme_bw()
Unconjugated-bisphenol serum profiles in a 60 kg adult female after a single 500 ng/kg BW peroral or dermal dose (reproducing Figure 4A).

Unconjugated-bisphenol serum profiles in a 60 kg adult female after a single 500 ng/kg BW peroral or dermal dose (reproducing Figure 4A).

Mass balance

A single peroral dose must be fully accounted for at long times as parent in urine plus the two conjugates in urine plus whatever remains in the body. With the renal and metabolic sinks the dose is essentially completely excreted.

mb <- solve_single(mod_bpaf, "BPAF", "oral", tmax = 240, dt = 0.5)
dose_nmol <- 500 * WT / mw[["BPAF"]]
body_states <- c(
  "stomach", "a_small_intestine", "a_liver", "a_fat", "a_gonads", "a_skin",
  "a_brain", "a_rapidly_perfused", "a_slowly_perfused", "plasma",
  "a_small_intestine_gluc", "gut_lumen_gluc", "serum_gluc",
  "a_small_intestine_sulf", "gut_lumen_sulf", "serum_sulf"
)
last <- mb[nrow(mb), ]
recovered <- last$urine + last$urine_gluc + last$urine_sulf +
  sum(as.numeric(last[body_states]))
data.frame(
  Quantity = c("Dose (nmol)", "Urine parent", "Urine glucuronide",
               "Urine sulfate", "Remaining in body", "Total recovered / dose"),
  Value = c(
    round(dose_nmol, 3), round(last$urine, 3), round(last$urine_gluc, 3),
    round(last$urine_sulf, 3), round(sum(as.numeric(last[body_states])), 3),
    round(recovered / dose_nmol, 4)
  )
) |> knitr::kable()
Quantity Value
Dose (nmol) 89.318
Urine parent 0.170
Urine glucuronide 88.574
Urine sulfate 0.573
Remaining in body 0.000
Total recovered / dose 1.000 

stopifnot(abs(recovered / dose_nmol - 1) < 0.01)

Assumptions, deviations and errata

  • Deposited code is authoritative; the packaged model reproduces it exactly. The Supplemental Material deposits the authors’ deSolve code and its input tables. The packaged rxode2 model reproduces that code to machine precision on both exposure routes for BPAF (1.566 pM peroral, 7.745 pM dermal). No parameter was tuned to any figure.
  • Peroral BPAF, code vs printed figure. The deposited code gives 1.566 pM for the BPAF single-dose peroral serum Cmax, 12% below the 1.79 pM the paper prints in Figure 4. The deposited input tables carry 2-3 significant-figure roundings of the full parameter set (the paper notes, e.g., that input table C2 mis-shifts one age group’s cardiac output), which is the most likely origin of the small difference; the figure’s exact simulation inputs are not recoverable from the published material. The model follows the deposited code and the printed Table 6-8 parameters.
  • Oral uptake period. Table 5 prints a peroral uptake period of 0.25 h, but the deposited code uses 3 min (period.O <- 3/60). The code value is used: 3 min reproduces the published peroral tmax, 0.25 h does not.
  • Conjugate volume of distribution. Table 8 lists the conjugate volume of distribution as 0.0435 L/kg (ICRP), but the deposited code sets it equal to the age-group plasma volume (VbodygC <- VplasmaC). The code is followed.
  • Single dermal source. The age-comparison figures (Fig 4/5) use one dermal source at a time (thermal paper, 20% absorbed, 6 h half-life). The second depot (depot_td_personalcare, 60% / 10 min) is present and is used when both sources act in parallel (the environmental-exposure scenarios); the validation here doses only the thermal-paper depot, matching Figure 4.
  • Inert sulfate gut pools. Gut sulfation is switched off in the deposited code (vmaxgutsC ~ 0, fguts = 0) and no sulfate undergoes EHR (met1s = 1), so a_small_intestine_sulf and gut_lumen_sulf stay empty. They are carried for structural symmetry with the deposited code and so that a future sulfation scenario can switch them on.
  • Molecular weights. BPA (228.28) is from the deposited code. The analog parent weights are derived from the Discussion’s glucuronide weights (BPS-g 426, BPF-g 376, BPAF-g 512 g/mol) minus the glucuronyl increment (176.12), giving 249.88, 199.88 and 335.88 g/mol. The BPAF value so derived reproduces the published serum Cmax, which independently checks it.
  • No IIV or residual error. The paper’s basic models are deterministic typical-value simulators; the 2D Monte-Carlo uncertainty/variability analysis (Tables S6-S7) is a separate layer not reproduced here. Inventing variances would be fabrication.
  • BPF / BPAF are uncalibrated. No in-vivo BPF or BPAF data exist; their uptake and clearance constants are inherited from the calibrated BPA (BPF) and BPS (BPAF) models, their EHR fractions are educated guesses from glucuronide molecular weight, and the paper states their predictions “should be considered with caution”.