Derivation of a Human In Vivo Benchmark Dose for Bisphenol A from ToxCast In Vitro Concentration Response Data Using a Computational Workflow for Probabilistic Quantitative In Vitro to In Vivo Extrapolation

A computational workflow which integrates physiologically based kinetic (PBK) modelling; global sensitivity analysis (GSA), Approximate Bayesian Computation (ABC), Markov Chain Monte Carlo (MCMC) simulation and the Virtual Cell Based Assay (VCBA) for the estimation of the active, free in vitro concentration of chemical in the reaction medium was developed to facilitate quantitative in vitro to in vivo extrapolation (QIVIVE). The workflow was designed to estimate parameter and model uncertainty within a computationally efficient framework. The workflow was tested using a human PBK model for bisphenol A (BPA) and high throughput screening (HTS) in vitro concentration-response data, for estrogen and pregnane X receptor activation determined in human liver and kidney cell lines, from the ToxCast/Tox21 database. In vivo benchmark dose 10% lower confidence limits (BMDL10) for oral uptake of BPA (ng/kg BW/day) were calculated from the in vivo dose-responses and compared to the human equivalent dose (HED) BMDL10 for relative kidney weight change in the mouse derived by European Food Safety Authority (EFSA). Three from four in vivo BMDL10 values calculated in this study were similar to the EFSA values whereas the fourth was much smaller. The derivation of an uncertainty factor (UF) to accommodate the uncertainties associated with measurements using human cell lines in vitro, extrapolated to in vivo, could be useful for the derivation of Health Based Guidance Values (HBGV).


PBPK MODEL
A human PBPK model was developed to study the fate of BPA following single oral doses. The initial model was based upon the PBPK model for the plasticizer DPHP described in McNally et al. (2012) with adaptions made as necessary. The final model described entry of BPA through ingestion with absorption of BPA from the gut and gastro-intestinal (GI) tract and a simple model of the lymphatic system describing uptake of BPA via the lacteals in the 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 absorbed dose minus the lymphatic and hepatic dose components. The model described the metabolism of BPA to BPAG and 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.
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 GI tract 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 four hours (to represent transport in bile) before BPA (and metabolites) appeared in the gut 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. Because of coding enterohepatic recirculation, the PBPK model was solved as a system of delay differential equations (DDEs) The final structure of the model descripted above followed the iterative model development process (incorporating uncertainty and sensitivity analysis) documented in McNally et al. (2021).
Default values for masses and flows, metabolism parameters, partition coefficients were based on literature sources and algorithms (Brown et al. 1997;ICRP 2002;Mazur et al. 2010;Schmitt 2008). Default values for uptake and elimination rates parameters were based on corresponding terms in McNally et al. (2021) and subsequently refined through tuning to available data. The default model assumed 80% of the administered BPA was absorbed through the hepatic route, with a further 5% absorbed through the lymphatic route with the complementary 15% passing through unabsorbed. .

STATISTICAL ANALYSIS Parameter Distributions
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 Table A1. The tabulated values are therefore based upon expert judgement and represent conservative yet credible bounding estimates.

Uncertainty Analysis
Uncertainty analysis was conducted throughout the model development process in order to identify coding errors, assess for deficiencies in the current version of the model, and to assess the bounding behaviour of the model as a function of its uncertain parameters; a formalised process for studying model behaviour and thus reducing the likelihood of coding errors is desirable for any PBPK model developed using a bespoke (open source) code (Loizou et al. 2021). An efficient approach using Latin Hypercube sampling was used throughout the model development process. For each new version of the PBPK model a 500-point maxi-min Latin Hypercube design was generated in order to efficiently cover the parameter space defined by the distributions ascribed to model parameter, the PBPK model was run for each of these design points, and various outputs from PBPK models were studied to assess the quantitative behaviour of the current version of the model, and the parameter limits in use at the time. The final uncertainty analysis reported on in this work was based on a 500-point maxi-min Latin Hypercube design and the parameter limits given in Table 2. This analysis used uniform distributions for all parameters, defined by upper and lower limits given in Table A1, with a total of 74 parameters studied. The outputs reported on in this work at the plasma concentrations of BPA, BPAG and BPAS (nM).

Sensitivity Analysis
Sensitivity analysis was conducted throughout the model development process in order to study the key model output sensitivities for each version of the model under development. A two-phased global sensitivity analysis comprising of elementary effects screening followed by a variance-based approach is generally desirable (Loizou et al. 2015;McNally et al. 2011), however during the model development and calibration phase of research, elementary effects screening alone was used to allow a more rapid analysis. In the final analysis reported on in this document the sensitivity analysis focussed on plasma concentrations of BPA, BPAG and BPAS (nM) with results at three time points: 0.5, 2-and 5-hours following ingestion of BPA.
The Morris Test computes two measures of sensitivity: µ * is the average effect of a parameter on model output; σ is a measure of a non-linear effect of a parameter on the response under study and/or interactions between parameters. For each output under study and at each time point under study the maxima of each of the µ * and σ over all parameters was determined, and all parameters where µ * and σ were within 5% of the maxima were recorded. Any parameter that BPA, BPAG or BPAS were sensitive to, at any time point and for either measure of sensitivity was considered as a variable for taking for forward into calibrationparameters passing this initial filter were subsequently reviewed with decisions for inclusion/exclusion of parameters in the calibration model based on expert judgement. Parameters that did not show sensitivity based upon the within 5% criterion described above were fixed at central values during calibration. This analysis used uniform distributions for all parameters, defined by upper and lower limits given in Table A1, with 74 parameters studied. Five elementary effects were computed for each parameter, leading to a design of 375 model runs.

Calibration
Calibration of a subset of sensitive model parameters using the HBM data of Thayer et al. (2015) was attempted. A Bayesian approach (McNally et al., 2012) was followed. This requires the specification of a joint prior distribution for the parameters under study, which is refined through a comparison of PBPK model predictions and measurements within a statistical model. The resulting (refined) parameter space that is consistent with the prior specification and measurements is the posterior distribution.
In the first instance individual model calibrations were performed using the blood data from each of the 14 volunteers in turn. The final calibration model used in this work used data from all 14 individuals from the HBM study of Thayer et al. (2015) with data on three specific outputs, concentrations of BPA, BPAG and BPAS (nM) in plasma, formally compared within the calibration model.
The final form of the statistical calibration model is given in equations A1-A3, where CPlasma BPA ij , CPlasma BPAG ij and CPlasma BPAS ij denote measurement i (at time t i ) for individual j (for j in 1:14) for the three respective model outputs, and μ BPA (θ, ω j ) ij , μ BPAG (θ, ω j ) ij and μ BPAS (θ, ω j ) ij denote the predictions from the PBPK model for these outputs corresponding to parameters (θ, ω j ). The vectors θ and ω j denote the global parameters common to all individuals (suitable for partition coefficients etc.), and participant specific parameters (suitable for accounting for variability in the physiology and modelling the participant specific uptake of BPA etc.). Prediction error for BPA, BPAG and BPAS was modelled through standard deviations , and respectively. Normal distributions, truncated at zero were assumed for all three relationships.
The prior distributions for model parameters are given in

Software
The PBPK model was written using GNU MCSim and compiled. In uncertainty and sensitivity analysis and for runs to produce fits following calibrations the executable was run using scripts called from RStudio. MCMC runs were conducted using GNU MCSim called from a command line environment. The DiceDesign package of R was used for generating Latin Hypercube designs. GSA of model outputs (through elementary effects screening and eFAST) were conducted using the Sensitivity package of R. The reshape2 package of R was used for reshaping of data for plotting and other processing of results.

Assessment of model form
As described in the statistical analysis section an iterative process of model development was conducted with uncertainty and sensitivity analysis informing improvements to model structure. The final model was tested using a 500-point maxi-min Latin Hypercube design and the parameter limits given in Table S1. The simulations corresponding to these design points are shown for concentrations of BPA, BPAS and BPAG in blood (nM) over the period 0-48 hours. The broad range of the runs covered the HBM data reported in Thayer et al. (2015) with the model form considered suitable for the purposes of this work.

Sensitivity analysis
Results from elementary effects screening are reported for BPA, BPAS and BPAG in Tables S2 to S4 for three time points, 0.5, 2 and 5 hours following ingestion of BPA. The mu and sigma measures at each time point are normalised with all results at the time point divided by the maximum value of mu or sigma at that time point; values are therefore between 0 and 1, with 1 denoting the most important for a given measure at a particular time point. The parameters where any of the measures were within 5% of the maximum (i.e. any value above 0.05) are highlighted in bold and considered for inclusion in the calibration model. This analysis was conducted for BPA, BPAS and BPAG.
In total 29 global and a further 11 local (individual specific parameters) were taken forward to calibration.