The biokinetics of bisphenol A following single oral doses were described using a PBK model previously developed for the plasticizer DPHP (McNally et al., 2021), with adaptions to model structure made, as necessary. The final model described entry of BPA through ingestion with absorption of BPA from the stomach and intestines and a simple model of the lymphatic system describing uptake of BPA via the lacteals in the small intestine and entering venous blood after bypassing the liver. The dose that entered the lymphatic system was coded as a fraction of the administered dose; a fraction of administered dose was coded as entering the liver via the portal vein; a fraction of the administered dose of BPA passed through the intestine without being absorbed and was coded as administered dose minus the lymphatic and hepatic dose components. The model described the metabolism of BPA-to-BPA glucuronide (BPAG) and BPA sulphate (BPAS) in both liver and gut. Sub-models were included to describe the kinetics of BPAG and BPAS, with the models for BPA and the two metabolites connected via gut and liver. Binding of BPA, BPAG and BPAS was coded from arterial blood, with the consequence that only the unbound fraction in blood was available for distribution to organs and tissues, metabolism, and elimination. However, the bound and unbound fractions of BPA in blood are in equilibrium such that the bound fraction is gradually removed to become unbound as the existing unbound fraction is excreted.

The model structures for BPA and the two metabolites differed only in the coding of uptake required for BPA (including the simplified description of a lymphatic compartment). BPA and metabolite models had a stomach and intestine draining into the liver. Adipose, blood, kidney, and slowly and rapidly perfused compartments were included. Elimination of BPA, BPAG and BPAS was coded through the kidney compartment, with first order elimination rates, proportional to kidney tissue concentration, coded in each case. Both BPA and metabolite models included the transport process of enterohepatic recirculation. Uptake of BPA and metabolites from the liver into bile was modelled as a first order uptake process with a delay of 4 hours (to represent transport in bile) before BPA (and metabolites) appeared in the small intestine and were available for reabsorption. First order elimination rates for each substance were coded to account for fractions of recirculated BPA and metabolites that were eliminated in faeces rather than reabsorbed from the small intestine. As a consequence of coding enterohepatic recirculation, the PBPK model was solved as a system of delay differential equations (DDEs).

The final structure of the model described above followed the iterative model development process (incorporating uncertainty and sensitivity analysis) documented in McNally et al. (2021). Table 1 lists a glossary of parameters name and abbreviations. A schematic representation of the model is shown in Figure 1.

The aim of the QIVIVE workflow was to estimate a distribution of ingested concentrations of BPA that was consistent with target tissue concentrations, for comparison against in vitro data. In vivo hepatic tissue BPA concentrations (CVli) and in vivo kidney tissue BPA concentrations (CVki) were selected as suitable target tissue concentrations because HepG2 cells are derived from the human liver and HEK293T cells from the human kidney, respectively and are considered in vitro surrogates for the liver and kidney in vivo, respectively.

The key to the approach described below is the recognition that the PBK model is an imperfect approximation to reality. Exact matching of the chosen PBK model response to an in vitro concentration suggests a higher degree of belief in the model than is warranted and is thus not desirable. By accepting a discrepancy between the two, within a specified threshold, model uncertainty is thus accommodated, and an error term is created that can be exploited by efficient sampling techniques. The modelling framework was adapted from Loizou et al., 2021, with the addition of the application of the VCBA for the estimation of the in vitro free concentration of BPA. Figure 2 shows a schematic of the workflow. The steps were:

1) Probability distributions for tissue volumes, expressed as a fraction of body weight, and tissue blood flows, expressed as a fraction of cardiac output were estimated using a virtual population generated using PopGen (McNally et al., 2014). Specifically, a US population of 10,000 individuals comprised of 2500 Caucasian males; 2500 Caucasian females, 2500 African American males and 2500 African American females, aged 25–45 and with BMI ranging from 19–35 was generated (which captured the characteristics of the human volunteer study population). Probability distributions for tissue volumes and fractional blood flows were based upon statistical analysis of this model output. For other parameters, such as partition coefficients and rate parameters, uniform distributions were ascribed based upon author’s judgement to represent conservative yet credible bounds and refined through the model development process. The probability distributions used in the reported uncertainty and sensitivity analyses and parameter calibrations are given in Supplementary Table SA1 Supplementary Material.

Finally, a PoD, the BMDL₁₀ lower bound in the in vivo dose response relationship was estimated (see section Calculation of in vivo benchmark dose).

The dose-response curves were predicted by calculating the AUC of the liver (CVli) or kidney (CVki) tissue concentrations versus fold induction or percentage activity. The first 3 h of the 24-h simulation period captured the maximum AUC for both tissues (data not shown). The mean, 2.5 and 97.5% of the credible interval values were calculated for the most sensitive parameters identified by GSA that determined CVli and CVki variability whilst also estimating the oral exposure concentration (PORALDOSE). PORALDOSE was estimated in µg/kg BW/day for direct comparison with the temporary Tolerable Daily Intake (t_TDI) of 4 μg/kg bw per day for BPA used in the risk assessment conducted by EFSA.

The BMD, BMDL and BMDU were estimated for the 10% BMRs compared to controls from the 90% confidence interval around the mean in vivo concentrations and corresponding fold inductions or percentage activities. The 10% BMR was calculated for direct comparison with the value used by EFSA. BMD values were calculated by model averaging from four fitted models for continuous (value for each individual) response data.

The model was coded in the GNU MCSim language (version 6.1.0.) ² and run under Windows 10 Pro using RStudio (RStudio Team 2016). Files for running MCSim under windows, tools and instructions for installation are available from Github ³ . All plots were created using R version 4.0.2 and ggplot2 (R Development Core Team 2008; Wickham 2016).

In order to perform probabilistic simulations the model code was further modified to ensure that logical constraints on mass balance and blood flow to the tissues were met by adopting the re-parameterizations described by (Gelman et al., 1996).

The PBK model was evaluated using RVis, an open access PBK modelling platform ⁴ which provides an intuitive user-friendly interface with which to interact with MCSim and the R platform ⁵ . The model equations were solved using MCSim which writes an output file in Tab Separated Values (TSV) format which is then input into the R environment and read by the R packages required for the various analyses. GSA of model outputs (Morris screening test and extended Fourier Amplitude Sensitivity Test (eFAST)) were conducted using the Sensitivity package of R. Reshaping of data and plotting was done using the reshape and ggplot2 packages respectively (Wickham 2007; Pouillot and Delignette-Muller 2010; Soetaert et al., 2010; Pujol et al., 2015). The main effects and total effects (McNally et al., 2011) were computed at each time point and parameter sensitivities were studied over this period using Lowry plots generated as described in McNally et al. (2011).

Benchmark dose values (BMDs) were calculated using PROAST version 69.0 hosted on the EFSA Open Analytics web site ⁶ .

The PBPK model is provided in Supplementary Material. R scripts for interacting with the model and replicating various stages of the analysis are available from the authors on request.

A Dell Precision M4800 with an Intel(R) Core™ i7-4800MQ CPU@2.70 GHz with 32.0 GB RAM running Windows 10 Pro was used for this study.

The parameter ranges for global and individual-specific parameters were refined through calibration of the model using the plasma BPA, BPAG and BPAS data of Thayer et al. (2015). Detailed results from calibration are provided in Supplementary Material; however, a summary of the key results is presented below. The simulations for three individuals only are shown in Figure 3: this subset of results demonstrates the range of behaviour that the calibrated model was able to fit and is used to emphasise key points from the results. Good fits to the plasma data of Thayer et al. (2015) were obtained for all 14 volunteers and the fits to the unique trends of BPA, BPAG and BPAS in plasma from each volunteer are shown in the three panels of Supplementary Figures S2–S12 of Supplementary Material. Summary statistics, that is, prior and refined (posterior) parameter distributions for global (i.e., common to all 14 individuals) and individual-specific parameters, are listed in Table 3 and Table 4.

Three individual-specific parameters governed variability in the uptake of BPA: the fraction of administered BPA that was available for hepatic uptake; the fraction taken up into the lymphatic system; and a lag term (Lymphlag) which determined the time duration before BPA entered venous blood (at the thoracic duct). Results from calibration (Table 4) showed the majority fraction of BPA, between 70 and 90% of administered BPA, was taken up in the gut and entered hepatic circulation with substantial differences in central estimates and credible intervals between volunteers. This majority fraction was subject to first pass metabolism in gut and liver, with simulations indicating that metabolism was rapid, with little BPA entering systemic circulation from this route. In contrast the fraction entering the lymphatic system was small – between 2 and 7% of administered BPA (with differences in central estimates and credible intervals between volunteers) - however this small fraction was not subject to first pass metabolism. Whilst only a small fraction, consideration of lymphatic uptake was important for successfully fitting 1) the peak of BPA in plasma occurring at later time point compared to metabolites BPAG and BPAS; 2) the double peak of BPA seen in individuals one and five (discussed below); 3) accounting for a slower clearance of BPA from plasma. The influence of the lymphatic component is evident in the plasma BPA data from individuals 1 and 5 (Figure 3) where the double peak fitted for these individuals is a consequence of the fraction entering the lymphatic system. For some participants (Figure 3, participant 3) there was no double peak – this behaviour corresponded to small values of Lymphlag. It is interesting to note that the double peak (characterised by the very rapid decrease and rapid increase at approximately 1.5 h) in individual one occurs during the uptake phase whereas it is in the distribution phase (at approximately 2.5 h) in individual five. This is consistent with a more rapid removal of BPA through metabolism (to BPAG and BPAS) due to the much higher estimated MPY in individual 1 compared to individual 5.

Comparison of the two sets of summary statistics indicates a broad consistency, with modest changes to the medians following calibration, but with a consistent narrowing of credible intervals. This was not true, however, for the following partition coefficients; Pgub, Pstb, PliG, PliS and PguS. For these parameters uniform priors over wide ranges were specified, reflecting the a priori uncertainty in partition coefficients predicted using an in-silico algorithm. Following calibration much narrower marginal posterior distributions were associated with these parameters, reflecting that good fits to HBM data could only be achieved within a much narrower range of parameter space. Similarly, there was a very substantial narrowing of the posteriors for KEMAX, BELLYPERM and GIPERM compared with priors.

The AUC tissue concentrations for CVli and CVki were studied with sensitivity analysis. The parameter ranges ascribed to model parameters were the refined limits following parameter calibration – for participant specific parameters, the minimum and maximum values were taken over the summary statistics for the 14 individuals. Due to the stochastic nature of the Morris test, parameter rankings were derived by identifying the mode for each parameter over six simulations. The ten parameters to which each of the outputs were most sensitive were selected from the entire set of 74 model parameters for the second phase of sensitivity analysis using the variance-based eFAST method. The outputs were judged to be insensitive to the other parameters; these were held fixed at default values in the second phase of sensitivity analysis. Table 5 lists the parameters following elementary effects screening that had the most significant impact on CVli and CVki tissue concentrations and are ranked according to eFAST.

Figure 4 shows the GSA for tissue concentrations of CVli and CVki as a Lowry plot (McNally et al., 2011). The Lowry plot shows the total effect of a parameter S_T, which is comprised of the main effect S_M (green bar) and any interactions with other parameters S_i (brown bar) given as a proportion of variance (range 0 – 1 on y axis) (McNally, et al., 2011). The ribbon (light blue), representing variance due to parameter interactions, is bounded by the cumulative sum of the S_M (lower bound) and the minimum of the cumulative sum of the S_T (upper bound). The most significant parameters that contributed to variance throughout the 24-h simulation period range from left to right. The 10 most significant parameters dominated from 0 to 5 h and remained in that order throughout that period whereas some swapping of order was observed for the other less important parameters. Furthermore, the S_T for the top three parameters for CVli (upper panel), and CVki (lower panel) accounted for 61 and 55% of variance, respectively. Parameter sensitivities from 0 to 5 h, the period of maximum variance and the period for which the AUC was calculated for use in QIVIVE, remained at similar proportions of variance throughout the simulation period of 24 h. To reduce simulation time, the top eight parameters for CVli and CVki listed in Table 5 were used in ABC MCMC.

Similar observations were made for plasma BPA, BPAG and BPAS where the top ten parameters shown in Table 5 accounted for 54, 82 and 66% variance, respectively.

The physicochemical parameters required to run the VCBA for BPA are listed in Table 6. The raw data from the Chemistry dashboard alongside the transformed nominal and calculated free in vitro concentrations and the ratios of free to nominal concentrations for all assays using the VCBA are presented in Table 7 and Table 8. The bioavailable concentrations of BPA were predicted to be around 48–49% of the nominal concentration (in the presence of 5% serum) for all the selected assays.

The percentage activity measure of response for assay OT_ER_ERaERb_0480 included 9 negative and two positive values at the lowest in vitro concentrations (Table 8). This occurs when activities are normalised using controls. This is a range of minimal activity associated with controls where the mean was assigned to zero. In most ToxCast assays a threshold point based on all controls is assigned and activity below the threshold is not considered to indicate anything about the test chemical. Therefore, values in this range were assigned zero on the assumption that the interpretation of this dataset does not change (Filer et al., 2014) (Dr. John Wambaugh, personal communication). A constant of one was then added to the entire range of percentage activities which was required for BMD analysis which assumes a lognormal distribution.

As described in the Materials and Methods and in previous reports (McNally et al., 2018; Loizou et al., 2021), a two-stage approach was used to sample PORALDOSE that was consistent with the free in vitro experimental data and a modest degree of model uncertainty. The latter was accounted for through accepting simulations within 5% of the target in vivo concentration. The first phase involves rejection sampling and is illustrated in Figure 5. Panel A of Figure 5 typically illustrates concentration-response profiles from 5,000 simulations, whereas panel B shows the concentration-response profiles from the retained simulations that were within 7.5% of the target concentration (0.011 mg/L, in this example). Plots like these were obtained for each in vitro concentration for each of the four assays listed in Tables 7, 8. In the second phase an ABC MCMC algorithm was applied for more efficient sampling of parameter space, in this case within a tighter threshold of 5%, consistent with a given in vitro target concentration. This two-stage process was repeated for each in vitro concentration with acceptance rates between 7.5 and 35%. Subsequent analysis was based upon results from the retained samples and pooled over the four chains run for QIVIVE for each in vitro concentration.

The posterior means and 95% credible interval for the in vivo dose-response relationships for PORALDOSE (µg/kg BW/Day) are provided in Table 9 for each target concentration in the assay.

The summary statistics for PORALDOSE were extrapolated from hepatic tissue concentrations (CVli) estimated from the in vitro concentration-response datasets; ATG_ERE_CIS_up (estrogen receptor activation), ATG_PXR_TRANS_up (pregnane X receptor binding), and from kidney tissue concentrations (CKli) estimated from OT_ER_ERaERb_0480 (estrogen receptor activation) and OT_ER_ERaERa_1440 (estrogen receptor activation).

The mean in vivo dose responses shown in Table 9 were used to calculate a BMDL₁₀ (lower limit of the 95% confidence interval on the BMR equivalent to a 10% effect size) for each in vitro assay (Figure 6). The mean, 2.5 and 97.5% percentile BMDL₁₀ values calculated for PORALDOSE were compared with the HED adjusted BMDL₁₀ and t_TDI values calculated by EFSA for relative kidney weight change in the mouse (Table 10). The mean BMDL₁₀ values for PORALDOSE showed 48-fold variability, ranging from 2.7 to 1300 µg/kg BW/day (Table 10). However, these BMDL₁₀ values may justifiably be adjusted by dividing with a chemical specific adjustment factor (CSAF) for inter-individual variability in pharmacokinetics (Bhat et al., 2017). Such a CSAF can be calculated by dividing the posterior 97.5% quantile by the median for a parameter that is responsible for the greatest variability in each output. In this case, the CSAF calculated for the fraction of BPA bound to plasma proteins, FB_BPA, is 2.52 (Table 3). The BMDL₁₀ of 516 µg/kg BW/day for the ATG_PXR_TRANS_up was similar to the HED adjusted BMDL₁₀ derived by EFSA. The BMDL₁₀ of 329 µg/kg BW/day is about 54% the HED adjusted BMDL₁₀ derived by EFSA. The application of an additional UF to an in vivo BMDL₁₀ calculated using an isolated human cell line to derive a t_TDI is a matter for discussion but could reduce these values further.

The BMDL₁₀ of 1.07 µg/kg BW/day with application of the CSAF or 2.7 µg/kg BW/day, without application of the CSAF, derived from the ATG_ERE_CIS_up assay was similar to the EFSA t_TDI of 4 µg/kg BW/day.

The mean in vivo BMDL₁₀ values for BPA calculated for the four ToxCast assays is considerable, ranging from 2.7 to 1300 µg/kg BW/day or 1.1–516 µg/kg BW/day with application of a CSAF of 2.52. The ratios of the higher to the lowest values of the latter (i.e., 516/1.1, 381/1.1 and 329/1.1) ranged from 299 to 469-fold which were far higher than the differences in the in vitro AC50 values which ranged from 3.2 to 43.1-fold (Table 10).

Also, the ratios of BMDL₁₀ calculated from free versus nominal concentrations were 0.48 for three of the four assays. Therefore, the ratio of free to nominal BMDL₁₀ is consistent with the ratio of in vitro free to nominal concentrations in three from four assays.

Typical predicted in vivo dose-response curves for PORALDOSE extrapolated from hepatic tissue concentration (CVli) for the in vitro datasets ATG_ERE_CIS_up (estrogen receptor activation), pregnane X receptor binding ATG_PXR_TRANS_up, and kidney tissue concentration (CVki) for the in vitro datasets OT_ER_ERaERb_0480 (estrogen receptor activation) and OT_ER_ERaERa_1440 (estrogen receptor activation) are shown in Figures 6 and 7.