High-throughput screening tools facilitate calculation of a combined exposure-bioactivity index for chemicals with endocrine activity.

High-throughput and computational tools provide a new opportunity to calculate combined bioactivity of exposure to diverse chemicals acting through a common mechanism. We used high throughput in vitro bioactivity data and exposure predictions from the U.S. EPA's Toxicity and Exposure Forecaster (ToxCast and ExpoCast) to estimate combined estrogen receptor (ER) agonist activity of non-pharmaceutical chemical exposures for the general U.S. population. High-throughput toxicokinetic (HTTK) data provide conversion factors that relate bioactive concentrations measured in vitro (µM), to predicted population geometric mean exposure rates (mg/kg/day). These data were available for 22 chemicals with ER agonist activity and were estimated for other ER bioactive chemicals based on the geometric mean of HTTK values across chemicals. For each chemical, ER bioactivity across ToxCast assays was compared to predicted population geometric mean exposure at different levels of in vitro potency and model certainty. Dose additivity was assumed in calculating a Combined Exposure-Bioactivity Index (CEBI), the sum of exposure/bioactivity ratios. Combined estrogen bioactivity was also calculated in terms of the percent maximum bioactivity of chemical mixtures in human plasma using a concentration-addition model. Estimated CEBIs vary greatly depending on assumptions used for exposure and bioactivity. In general, CEBI values were <1 when using median of the estimated general population chemical intake rates, while CEBI were ≥1 when using the upper 95th confidence bound for those same intake rates for all chemicals. Concentration-addition model predictions of mixture bioactivity yield comparable results. Based on current in vitro bioactivity data, HTTK methods, and exposure models, combined exposure scenarios sufficient to influence estrogen bioactivity in the general population cannot be ruled out. Future improvements in screening methods and computational models could reduce uncertainty and better inform the potential combined effects of estrogenic chemicals. Published by Elsevier Ltd.

Introduction

Environmental exposures to multiple chemicals with endocrine activity are a concern for public health (Gore et al., 2015). Although humans are simultaneously exposed to many environmental chemicals (Rappaport, 2012; Wishart et al., 2015; Woodruff et al., 2011) that can act collectively to alter endocrine function (Kortenkamp, 2014), hazard is typically evaluated for just one chemical at a time (NRC, 2009a). This is due in part to the legal requirements for chemical risk assessment (for example, the Federal Insecticide Fungicide, Rodenticide Act, the Toxic Substances Control Act, the Safe Drinking Water Act) (Monosson, 2004). The uncertainty factors applied in chemical assessments for individual chemicals are not intended to account for the combined effects of multiple chemical exposures (Martin et al., 2013). Understanding the combined impacts of multiple exposures and incorporating these predictions into chemical evaluations will therefore be an important step towards translating emerging science to support regulations (NASEM, 2017; Gwinn et al., 2017).

Recent progress in the development of high-throughput assays and computational models for human endocrine activity provide tools to estimate combined endocrine bioactivity for hundreds of diverse chemicals. For example, a computational model of estrogen receptor (ER) bioactivity (Judson et al., 2015) integrates results from 18 different in vitro assays from the EPA’s Toxicity Forecaster (ToxCast) and the U.S. Federal Tox21 collaboration (Tice et al., 2013). The ER model predicts ER activity measured in validated guideline assays (Kleinstreuer et al., 2015) with high accuracy (Browne et al., 2015). High throughput toxicokinetics (HTTK) based on human in vitro protein binding and metabolism data are currently available for a third of the chemicals tested in ToxCast and can be used to convert ToxCast bioactivity data from an in vitro concentration to an equivalent in vivo exposure (Wetmore et al., 2015). Integration of bioactivity and kinetic data has been proposed for prediction of biological-pathway-altering doses in vivo to inform high-throughput chemical assessments (Judson et al., 2011). In parallel, EPA’s ExpoCast Systematic Empirical Evaluation of Models (SEEM) predicts human exposure to ToxCast chemicals based on chemical application (for example, use in consumer products) and production volume data (Wambaugh et al., 2014). Together, these high-throughput models of bioactivity, toxicokinetics and exposure provide an integrated platform to predict combined bioactivity of general population exposures to multiple chemicals acting through a common pathway of toxicity.

As regulatory agencies develop strategies for addressing effects of mixtures (Bell and Edwards, 2015; Bopp et al., 2015; EFSA, 2013; NRC, 2009b; CPSC, 2014; Kienzler et al., 2016), opportunity exists to incorporate high-throughput bioactivity and exposure predictions into these strategies. To that end, we explored the potential to apply existing high-throughput data to understand potential combined effects of multiple chemical exposures. We propose a combined exposure-bioactivity index (CEBI) that calculates the sum of exposure/bioactivity ratios for multiple chemicals with estrogen receptor bioactivity using high-throughput bioactivity, kinetic and exposure predictions from ToxCast, HTTK modeling, and ExpoCast at varying levels of confidence. The CEBI approach is intended to estimate overall bioactivity of combined chemical exposures predicted in humans or environmental media. Here we describe the potential combined estrogen receptor agonist bioactivity of exposures predicted in the general population, estimate the uncertainty surrounding current high-throughput predictions, and highlight future research needed to address the limitations of current predictive models.

Methods

Description of high-throughput data sources

The ToxCast ER model (Judson et al., 2015) integrates results of 18 independent in vitro high throughput screening assays to predict overall ER activity, agonist activity or antagonist activity on a scale of 0–1. A score of 1 corresponds to a chemical with equivalent efficacy and potency to 17-alpha ethinyl estradiol across all 18 assays, while a score 0 indicates no activity across all 18 assays. The current analysis focused on non-pharmaceutical, synthetic chemicals with ToxCast ER agonist scores ≥0.1. This threshold was found to identify ER active reference chemicals with a balanced accuracy ranging from 84 to 93% (Browne et al., 2015). Although the overall ER model provides a validated prediction of ER bioactivity, it does not provide quantitative predictions of chemical potency. Therefore, the minimum (most sensitive assay), and median (across all active assays) active concentrations from the 18 ToxCast/Tox21 assays were used to inform the current analysis. Concentrations associated with different levels of effect, the AC10 (10% of maximum effect), AC50 (50% of maximum effect) or ACC (lowest significant effect) for each chemical were identified from the concentration-response curves modeled for each of the 18 assays using the Hill equation: where the activity (A) for chemical i is a function of chemical concentration([C]), the estimated concentration at which 50% of maximum activity occurs (AC50,), the maximum efficacy (E), and Hill slope (h) (Judson et al., 2015). Maximum effects in most assays were defined relative to 17-beta estradiol positive controls as defined in (Judson et al., 2015). The median and minimum AC10, AC50, and ACC values (referred to in this manuscript as AC10min, AC10med, AC50min, AC50med, ACCmin, and ACCmed) across all assays with significant effects (as defined in (Judson et al., 2015) were used to provide a range of bioactivity values for each chemical, reflecting variability across assays and different levels of biological effect (Supplemental Table 1). The ToxCast data were analyzed as presented in the supplementary material of Judson et al. (2015). All data (including additional assay information) are available from the EPA website (https://www.epa.gov/chemical-research/exploring-toxcast-data-downloadable-data) under the heading “High-throughput Screening Data for Estrogen Receptor Model”.

The ExpoCast SEEM exposure model (Wambaugh et al., 2014) uses data on chemical production volume and chemical use categories to predict population-level, geometric mean intake rate (mg/kg/day) with defined levels of confidence (e.g., median or upper 95% confidence interval estimates) (Supplemental Table 1). Briefly, chemical use was determined algorithmically by counting occurrence of chemicals on >2000 lists that comprise the Aggregated Computational Toxicology Repository (ACToR: http://www.epa.gov/actor/). The model was calibrated by linear regression to describe mg/kg/day intake rates inferred from the National Health and Nutrition Examination Survey (NHANES) biomonitoring data on urinary concentrations of analyte chemicals (either the chemical itself or metabolites that have been mapped to the parent chemical) (CDC, 2012). A reverse toxicokinetics model assumed that the analyte molecules in urine measured by NHANES were the result of steady-state exposures to parent chemicals (Ring et al., 2017; Wambaugh et al., 2014). Confidence intervals on the geometric mean were determined based on the chemical-to-chemical variability in intake rate that was not explained by the SEEM model. For those chemicals included in the NHANES the actual inferred exposure rates were used (as reported in Ring et al. (2017)). For all other chemicals the SEEM predictions were used. The SEEM model is not designed to predict exposure to naturally occurring chemicals (such as genistein and zearalenone). Such chemicals were excluded from the current analysis.

R package “httk” v1.9.2 (https://CRAN.R-project.org/package=httk) was used with default settings to perform in vitro-in vivo extrapolation (Pearce et al., 2017). The script used is provided as supplemental material. The httk tool has previously been used to predict high throughput toxicokinetics (HTTK) (Ring et al., 2017; Sipes et al., 2017; EPA, 2014). The extrapolation relies on a simple toxicokinetic model to predict plasma concentration using in vitro data. For each chemical two parameters are measured: protein binding to human plasma samples and metabolism due to primary human hepatocytes in suspension (Rotroff et al., 2010; Wetmore et al., 2012; Wetmore et al., 2015). Total chemical clearance rate from the body can be estimated from these parameters (Pearce et al., 2017) and used in an algebraic equation to predict steady-state plasma concentration (Css) (Wilkinson and Shand, 1975). In the model, bioavailability was assumed to be 100% and hepatic metabolism was assumed to be restricted to the free fraction of chemical in plasma. A Monte Carlo population simulation was used to estimate the median and 95th percentile Css resulting from a constant 1 mg/kg bodyweight/day exposure rate. Css was used for in vitro-in vivo extrapolation via “reverse dosimetry” (Tan et al., 2007) conversion of ToxCast in vitro nominal concentrations (C, μM) to an equivalent in vivo dose rate (ED, mg/kg/day) (Wetmore et al., 2015):

In this analysis, only the median HTTK value is used. The 95th percentile data is reported in the Supplemental Material and all figures may be recreated using that value. The 95th percentile is by necessity a higher value, which results in lower equivalent in vivo dose rate predictions.

Prediction of ER agonist bioactivity of a mixture using a combined exposure-bioactivity index

We use bioactivity measured in ToxCast assays to calculate a ‘Bioactivity Quotient (BQ)’ analogous to a Hazard Quotient: where ER is the exposure rate (mg/kg/day) predicted for chemical i by the ExpoCast SEEM consensus model. To calculate ratios of exposure to bioactivity, Eq. (2) must be used to convert ToxCast bioactivity data to units comparable to SEEM exposure estimates (mg/kg/day). The particular bioactive in vitro concentration (C) from ToxCast that is used in Eq. (2) can vary (for example, AC50, AC10, or ACC). Likewise, the quantile used for HTTK (for example, median or 95th) Css can reasonably be varied. The chemical-specific HTTK data required for comparison of in vitro bioactivity to in vivo exposure predicted by SEEM were only available for 22 of the 59 non-pharmaceutical, synthetic compounds with ER agonist bioactivity in the ToxCast ER model. For chemicals lacking experimentally informed HTTK values, the geometric mean of HTTK values available for the subset of 22 ER bioactive chemicals was used to estimate Exposure/Bioactivity ratios. Using the geometric mean assumes that the HTTK estimates for chemicals with and without experimentally derived HTTK values have a similar overall distribution. This approach will over- or underestimate HTTK values for individual chemicals but offers an approximation of overall predictions for combined ER activity. Bioactivity quotients calculated from in vitro bioactivity data were then summed to calculate a ‘Combined Exposure-Bioactivity Index’ (CEBI) analogous to a Hazard Index (Fig. 1): where n = 59 includes all the chemicals analyzed here.

Fig. 1.

Method for calculating the Combined Exposure-Bioactivity Index (CEBI). For each chemical with ER bioactivity in ToxCast, bioactive concentrations (AC50, AC10, ACC) were determined for median and most sensitive (minimum) assay responses. Hypothetical concentration-response curves are shown above for illustration and do not reflect actual ToxCast data. ExpoCast produces median and upper bound (95% confidence interval) model estimates of population geometric mean exposure rates for each chemical. Hypothetical exposure distributions are included for illustration, but the true shape of the population exposure distribution and exposure rates for the 95th percentile of the population are unknown. Bioactivity values (μM) were multiplied by a high-throughput toxicokinetic (HTTK) conversion factor to put predicted exposures and bioactivity values into equivalent units. Ratios of exposure rate and HTTK-adjusted bioactivity are calculated for each chemical using all combinations of values. The CEBI is the sum of these individual chemical exposure/bioactivity ratios.

Assuming an additive model for estrogen receptor agonist activity (Kortenkamp, 2007), we calculated CEBIs for different combinations of bioactivity and predicted exposure values, reflecting varying levels of biological effect and model confidence. CEBIs were calculated for the subset of 22 chemicals with experimentally-informed HTTK values as well as for the full set of 59 ER active chemicals. Because the estrogen receptor bioactivity data from ToxCast ER assays are provided in terms of the percent of estradiol activity (i.e. percent maximum activity), this approach is similar to applying a toxic equivalency factor approach (EPA, 1994). To characterize the degree to which individual chemicals influence CEBIs, we calculated the percent of each CEBI contributed by each chemical.

Prediction of ER agonist bioactivity of a mixture using a concentration-addition model

Combined estrogen bioactivity was also calculated in terms of the percent maximum bioactivity of chemical mixtures in human plasma using a concentration-addition model. By inverting Eq. (2), Css defines the relationship between the exposure rate ER (mg/kg/day) and the predicted plasma concentration (C, μM) for each chemical at steady state (Pearce et al., 2017; Ring et al., 2017; Wambaugh et al., 2015).

For the concentration-addition model, the conversion factor was applied to convert mg/kg/day exposures predicted by SEEM from ExpoCast to an estimated corresponding plasma concentration (Eq. (5)) to estimate total receptor bioactivity of ‘mixtures’ in plasma based on ToxCast AC50 value. Concentration-addition models predict the percent maximum activity of a mixture of chemicals acting additively through a common mechanism of action and have been shown to predict the combined effects of estrogenic chemicals (Kortenkamp, 2007). Here, the combined ER bioactivity of 59 chemicals was estimated using a concentration-addition model (Eq. (6)) (Bermudez et al., 2010; Rider and LeBlanc, 2005). equals the % response of the mixture relative to the activity elicited by 17-beta estradiol (E2) (maximum receptor activity), C is the estimated steady-state plasma concentration based on SEEM exposure predictions of each individual chemical in the mixture (calculated in Eq. (3)), and 〈AC50,〉 is either the median or minimum concentration of the chemical causing a 50% maximum response observed across active ToxCast assays for each chemical. The average of the hillslope values from the dose response curve of chemicals in the mixture is assumed to be equal to one for all chemicals and hillslope is therefore not included in the equation. The concentration addition estimations were performed using both median and minimum ToxCast AC50s for each chemical and the estimated steady state plasma concentrations inferred from the predicted exposures of the geometric mean of the population and the upper bound 95% confidence interval of the ExpoCast SEEM model. This calculation yields a prediction of the percent ER binding activity of the mixture of 22 chemicals based on bioactivity predicted by ToxCast and steady state plasma concentrations based on HTTK-adjusted, ExpoCast SEEM predictions. We used the concentration-addition model to calculate ER bioactivity of predicted exposures to individual chemicals as well as for the mixture of 59 chemicals.

Results

Among the 1812 chemicals screened in the ToxCast ER model, 59 non-pharmaceutical synthetic chemicals had ER agonist bioactivity (defined here as a ToxCast ER agonist model score ≥ 0.1) (Chiu and White, 2006; Judson et al., 2015). Of these 59 chemicals, 22 had high-throughput toxicokinetic (HTTK) data available to facilitate comparisons between ToxCast data on in vitro bioactivity to ExpoCast SEEM predictions of in vivo human exposure. We calculated CEBIs (the sum of exposure/bioactivity ratios) for 22 chemicals with HTTK data (Eq. (4); Fig. 1) based on different levels of bioactivity and exposure model confidence (Table 1). The population 95th percentile Css is calculated and reported in the Supplemental Material but was not used here. A CEBI greater than one indicates that the predicted combined exposures are sufficient to exceed the defined level of bioactivity (i.e., concentrations exceeding the AC10, AC50, or ACC). For each combination of bioactivity and exposure values, the resulting CEBIs range from over four orders of magnitude below one to three orders of magnitude above one (Fig. 2), depending on the bioactivity values and exposure predictions used. CEBIs for the broader set of 59 chemicals with ER agonist bioactivity in ToxCast were approximated using the geometric mean of existing HTTK values to approximate HTTK for those chemicals lacking experimentally derived values (also Fig. 2). In the case of the salt 3,3′-Dimethylbenzidine dihydrochloride (which does not have HTTK data), the measured HTTK data for 3,3′-Dimethylbenzidine were used.

Table 1

Summary of parameters used to generate CEBIs for different levels of bioactivity and exposure.

High-throughput Toxicokinetic Estimates (HTTK)	ExpoCast SEEM Exposure Model Outputs	ToxCast Bioactivity Values
For chemicals with HTTK data: median HTTK model estimatesFor chemicals lacking HTTK data: geometric mean of median HTTK values for other chemicals	Median model estimate of the population geometric mean exposure rate	ACC (min)
		ACC (med)
		AC10 (min)
		AC10 (med)
		AC50 (min)
		AC50 (med)
	Upper bound (95% confidence interval) estimate of the population geometric mean exposure rate	ACC (min)
		ACC (med)
		AC10 (min)
		AC10 (med)
		AC50 (min) AC50 (med)

AC50 = Concentration at which 50% of maximum activity is observed.

AC10 = Concentration at which 10% of maximum activity is observed.

ACC = Concentration at which activity in significantly different from control.

min = minimum response across all ToxCast assays.

med = median response across all active ToxCast assays.

Fig. 2.

Range of Combined Exposure-Bioactivity Indexes (CEBI) estimated from bioactivity and exposure model outputs. For 22 non-pharmaceutical synthetic chemicals with ToxCast ER agonist activity and high-throughput toxicokinetic data available, CEBIs were calculated from all possible combinations of ExpoCast exposure estimates (median and upper bound population geometric mean exposure rates) and ToxCast bioactivity values (median and minimum AC50, AC10 and ACC across active ToxCast assays) (dark diamonds). To expand CEBI calculations to include all 59 non-pharmaceutical synthetic chemicals with ER agonist bioactivity in ToxCast (light diamonds), the geometric mean of HTTK values available for the subset of 22 chemicals were used to approximate exposure/bioactivity ratios for the remaining chemicals. CEBIs above one (indicated by the red line) indicate conditions where combined exposures exceed those expected to produce estrogen bioactivity at the designated level of biological effect (i.e. significant, 10%, or 50% receptor activity).

High throughput methods for addressing population variability are currently only available for toxicokinetics (Ring et al., 2017). Present methods do not allow analysis of population variability for bioactivity or exposure in a high throughput manner. Across the 22 chemicals with HTTK data, the ratio of the median to the 95th percentile most “susceptible” individuals (that is, those with higher plasma concentrations for the same intake rate) indicated that the susceptible individuals would require a dose 5 times lower than the median individual. This would uniformly shift all CEBIs slightly higher in Fig. 2. However, for consistency with exposure and bioactivity estimates the median HTTK was used for the analyses reported here.

The relative contribution of each chemical to total combined bioactivity varies with the levels of exposure model confidence and biological effect. CEBIs calculated from median (Fig. 3A) and upper bound (Fig. 3B) estimates of geometric mean exposure are dominated by a relatively small subset of the chemicals included in the analysis. The median geometric mean exposure estimates (Fig. 3a) identified phenyl paraben, kepone, 3,3-dimethylbenzidine dihydrochloride, bisphenol A, 4-methylpent-3-en-2-one, benzyl butyl phthalate and flumetralin as significant contributors to combined activity. When the upper bound geometric mean exposure estimates are used (Fig. 4b), the particularly large CEBIs generated from minimum AC50, AC10, or ACC bioactivity values are largely driven by hydramethylnon, 2,2-bis(4-hydroxyphenyl)-1,1,1-trichloroethane, 3,3-dimethylbenzidine dihydrochloride, and branched 4-nonyphenol, which have very low minimum bioactivity values in specific assays.

Fig. 3.

Relative contributions of 59 chemicals to Combined Exposure-Bioactivity Indexes (CEBI) based on different bioactivity values. Percent contribution of individual chemicals to each CEBI is shown for CEBIs calculated from a range of bioactivity values and (A) median geometric mean exposure estimates or (B) upper bound geometric mean exposure estimates. *Indicates chemicals for which no chemical-specific HTTK data were available.

Fig. 4.

Estimated bioactivity of ER bioactive chemicals at predicted plasma concentrations. A concentration-addition model was used to calculate the percent of maximum ER bioactivity expected for individual chemicals (only those with > 0.5% maximum ER activity are shown above) and for combined bioactivity of all 59 chemicals combined (“COMBINED” reported in the last column) based on exposure and bioactivity predicted by ExpoCast and ToxCast. HTTK conversion factors were used to estimate plasma concentrations for each chemical based on ExpoCast median and upper bound (95% confidence interval) estimates of the population geometric mean exposure rate. Percent maximum ER bioactivity at these predicted plasma concentrations was calculated based on the AC50 minimum (AC50min) and median (AC50med) across ToxCast assays, assuming a concentration-response curve Hill equation slope of 1 for all chemicals. *Indicates chemicals for which no chemical-specific HTTK data was available.

Alternatively, the CEBIs can be presented in terms of the percent of maximum ER bioactivity of chemical concentrations predicted in human plasma. Using a concentration-addition model, we calculated the percent maximum ER bioactivity expected for a mixture of all 59 chemicals at plasma concentrations extrapolated from median and upper bound exposure estimates using the HTTK conversion factor (Fig. 4). At steady-state plasma concentrations extrapolated from median estimates of geometric mean exposures, the concentration-addition model predicted additive bioactivity to be 1.5% of maximum based on the AC50min and 0.01% of maximum based on the AC50med. At steady-state plasma concentrations extrapolated from upper bound model estimates of geometric mean exposures, additive bioactivity was predicted to be 689% and 38% of maximum based on the AC50min and AC50med, respectively.

Discussion

Environmental chemical exposures occur in combinations, posing a challenge to toxicologists, exposure assessors, and regulators worldwide (Thomas et al., 2019; Kienzler et al., 2016). This analysis provides a model for applying high-throughput data to quantify combined bioactivity from diverse environmental chemical exposures using the estrogen pathway as an example. The range of CEBI predictions vary from several orders of magnitude above 1 (indicating potential for estrogen receptor bioactivity at predicted exposure levels) to several orders of magnitude below 1, depending on the assumptions made (Fig. 2). The wide range of results reflects limitations and uncertainties of current methods for hazard, exposure, and toxicokinetics (summarized in Table 2). Likewise, the contribution of any one chemical to the overall activity similarly varies depending upon assumptions (Fig. 3). However, certain chemicals contribute significantly more to combined activity than most others regardless of assumptions; these chemicals and mixtures thereof are obvious testing priorities.

Table 2

Selected Sources of Uncertainty and Variability.

Factor	Description	Uncertainty and variability captured in this analysis	Uncertainty and variability remaining
Bioactivity  ACC, AC10, AC50 (μM)	In vitro potency	• Uncertainty: Use of AC50, AC10 and ACC captures multiple levels of bioactivity• Uncertainty and Variability: Use of median and minimum bioactivity values across the 16 assays associated with ER agonism• Co-Variability: Mixture additivity	• Uncertainty and variability: The level of receptor occupancy required to trigger an effect• Uncertainty: Uncertainty in ToxCast curve-fits• Uncertainty: The extent to which activity of each chemical is altered in the presence of endogenous estrogen• Uncertainty and variability: Metabolic activation/deactivation of chemicals• Uncertainty and variability: In vitro distribution of chemicals• Uncertainty and Variability: ToxCast assays and ER model may not capture allaspects of ER agonist bioactivity• Variability: Differences in population in response to ER agonist activity• Co-Variability: Other mixture effects (synergy/antagonism)
HTTK conversion  Css (μM/1 mg/kg/ day)	Estimate of the steady state plasma concentration associated with a given rate of exposure (the conversion factor from dose to concentration)	• Variability: The HTTK model uses a Monte Carlo approach to accounts for some aspects of toxicokinetic differences between adults including distribution	• Uncertainty: Experimental noise in assays used to estimate rates of chemical metabolism and protein binding• Uncertainty: in the HTTK model• Variability: Population differences absorption, metabolism, and excretion.• Co-Variability: Correlations in physiology that impact toxicokinetics• Variability: Juvenile/geriatric populations
Exposure Estimate  GM (mg/kg/day)	Estimate of population geometric mean of constant daily exposure	• Uncertainty: Exposure rates include 95% confidence interval for predictions the population geometric mean	• Variability: Population variability in exposure, particularly high-end users and occupational exposures• Co-Variability: Exposure to some combinations of chemicals are much more prevalent than others

Depending on how the CEBI is calculated, it is possible that combined environmental chemical exposures are sufficient to have significant estrogen receptor bioactivity. The wide range of CEBI values from current model predictions reflect different degrees of biological effect and model uncertainty. If there were a balance between over- and under-prediction of exposure, toxicokinetics, and bioactivity across chemicals, median estimates might provide a reasonable approximation of the CEBI. Using the most potent bioactivity values for all 59 chemicals considered in this analysis provides larger CEBIs by likely overestimating for potential estrogen receptor agonist bioactivity. Using the upper bound of population geometric mean exposure for all chemicals assumes that the exposure estimates of all chemicals are biased in the same manner. Using the upper bound would also assume that individuals are exposed to all chemicals at or above population geometric mean exposures. The chance of any one person being at or above the median for all 59 chemicals might be extremely low, roughly 1 in 1017, given that for n uncorrelated chemicals, only 1/2n people are expected to have the median concentration. However, Fig. 3 shows that while the identities of the most significant contributors to a CEBI depend on the assumptions used, the most significant contributors were always a small subset of chemicals. Therefore, examination of co-occurrence should focus on these 2–6 chemicals for which 1.5% to 25% of the population might be expected to be above the median for uncorrelated exposures. The analysis here is for the geometric mean general population; occupationally exposed individuals and heavy consumer product users might have CEBIs that are substantially higher than population geometric means and require further studies.

In addition to the hazard index-like CEBI, we also used a mixtures toxicology approach, (i.e. the concentration-addition model) to predict overall ER bioactivity of a hypothetical mixture of 59 chemicals. Predictions of the concentration-addition model were generally consistent with calculated CEBIs. For example, AC50med bioactivity values and upper bound geometric mean exposures for the 59 chemicals yield a CEBI of 0.4, indicating that predicted combined exposures are about half the level predicted to have 50% of maximum activity. This is consistent with concentration addition model estimates that total bioactivity of the equivalent mixture would be 38% of maximum bioactivity.

Future model optimization

Several assumptions and limitations in the current models could lead to an underestimate of CEBIs. Ongoing refinement of ToxCast, HTTK, and ExpoCast data and models will likely reduce uncertainty by incorporating more endocrine mechanisms, eliminate current assumptions in toxicokinetic modeling, integrate additional exposure pathways, and address susceptible populations and sensitive life-stages. For example, the current ToxCast ER model is primarily based on nuclear receptor activity and the predictive ability for other mechanisms of estrogen disruption are uncertain. A comprehensive model of endocrine bioactivity would incorporate all relevant pathways, including steroidogenesis, metabolism, and nuclear and non-nuclear receptor signal transduction. Prediction of bioactive metabolites from inactive parent compounds can also be used to expand application to many thousands of chemicals (Pinto et al., 2016).

A TK model is necessary to compute a linear conversion factor relating in vitro bioactivity to predicted exposures. The linear conversion factor is only valid for conditions of first-order metabolism and steady-state plasma concentration (Wetmore et al., 2012) and assumes that bioactive in vitro concentrations are equal to bioactive plasma concentrations. Furthermore, the model used to predict steady-state plasma concentrations assumes 100% oral absorption, passive renal clearance, and first-order hepatic-only metabolism, and ignores tissue partitioning (Wilkinson and Shand, 1975). These simplifying assumptions allow a generic TK model (with the same form for all chemicals) to be rapidly parameterized for large numbers of chemicals using only measurements of intrinsic hepatic clearance and plasma protein binding, both of which can be measured using high-throughput in vitro assays (Rotroff et al., 2010; Wetmore et al., 2012; Wetmore et al., 2015). The errors introduced by these simplifying assumptions are generally expected to underestimate the bioactive dose for most chemicals (therefore overestimating the CEBI), but the degree of confidence in the specific model used here (HTTK, Pearce et al. (2017)) can in some cases be predicted from chemical properties (Wambaugh et al., 2015). Meanwhile, the SEEM model for chemical intake rates is calibrated to human biomarker data and therefore reflects exposure from all routes. It does not take into account the fraction absorbed, which may result in an underestimate of actual exposure (Wambaugh et al., 2014)

Here we used toxicokinetics for the median of the population. Population variability has been incorporated into the HTTK model predictions using a Monte Carlo approach that samples from individuals profiled by the NHANES, such that factors like age, gender, ethnicity, and body size inform co-varying properties such as cardiac output and liver volume (Ring et al., 2017). However, because the metabolism and plasma protein binding data used to parameterize the model were measured in pooled adult tissue samples, the model is limited in its ability to predict the effects of variability in metabolism and plasma protein binding across potentially sensitive subpopulations. Metabolic variability can be simulated by assuming an overall population coefficient of variation of 0.3 for metabolic rate and further assuming that a small fraction of the population were ultra-low metabolizers (Ring et al., 2017). However, metabolism for any one chemical may be a function of multiple enzymes, and the expression levels of those enzymes can co-vary significantly within and between populations (Bois et al., 2010; Dorne, 2007) so that no one individual is likely to be the 95th percentile most-sensitive for all chemicals.

SEEM predictions describe population geometric mean exposure rates based on chemical use and production volume data. The model is based upon inference from NHANES urine biomonitoring data, and the current ExpoCast model is built upon the assumption that the concentrations of chemical analytes in urine reflect steady-state exposure to chemicals with 100% absorption. These assumptions may underestimate overall exposure. In addition, though geometric means can be estimated for various demographic groups, the model does not yet capture population variability within those groups. The median and upper bound geometric mean exposure estimates reflect model uncertainty, but there are no predictions of the high end of the population exposure distribution (e.g. 95th percentile exposures). Although 95th percentile estimates cannot be predicted without additional knowledge of the shape of the exposure distribution, ongoing work to generate life-stage or demographic-specific exposure scenarios will help to identify more highly exposed individuals (Ring et al., 2017). A major limitation is that the current model assumes all chemical exposures are independent. In reality, many exposures are likely to be correlated (Kapraun et al., 2017; Tornero-Velez et al., 2012). As ExpoCast is iteratively refined to address these shortcomings, it may ultimately be able to support probabilistic exposure modeling to predict combined estrogenic exposures for specific populations. Additional research would be useful on combinations of estrogenic chemicals that are likely to have co-exposure.

There are several assumptions made in integration of these high throughput datasets to calculate the CEBI. Estrogen receptor activity of multiple chemicals is assumed to be additive (Kortenkamp, 2007). The additive nature of multiple chemicals could be demonstrated directly in ToxCast assays by generating concentration-response curves for defined mixtures. Furthermore, the current approach does not explicitly account for partial agonists. While future studies could account for partial agonists by incorporating differences in the shape and slope of concentration-response curves for each chemical (Howard et al., 2010) into the concentration-addition model, this would be complicated by differences in concentration-response relationships across ToxCast assays. Finally, the current analysis is limited to non-pharmaceutical, synthetic chemicals evaluated in ToxCast which may not fully represent the chemical diversity of the estimated 10,000 chemicals that might require screening in the EDSP (EPA, 2012). More comprehensive consideration of combined estrogenic exposures should account for estrogenic drugs and phytoestrogens as well as chemicals that have not yet been tested in ToxCast.

The CEBI approach could be applied to predict overall bioactivity of combined exposures measured directly in ecological or human biomonitoring samples. Rapidly improving analytical chemistry technology will soon facilitate non-targeted screening for an increasing proportion of the exposome (Wild, 2005), providing a powerful platform for characterization of combined exposures in specific populations (Smith et al., 2015). As methods to quantify the exposome are incorporated into biomonitoring, ToxCast bioactivity could be related directly to combined exposures measured in specific populations of interest. This would be a particularly valuable comparison because human plasma concentrations measured in biomonitoring studies could be directly related to in vitro concentrations (assuming in vitro concentration = plasma concentration), eliminating the need for HTTK models that introduce additional uncertainty. Such comparisons would facilitate evaluation of ‘pathway-based’ hypotheses for epidemiology, in which effects are associated with combined exposures targeting a pathway of interest as opposed to exposure to individual chemicals. In addition, by anchoring pathway-based in vitro data to epidemiological outcomes, this type of analysis could serve to further validate high-throughput models.

Given the well-characterized relationship between estrogen bioactivity and downstream phenotypic outcomes, the combined ER bioactivity predicted in this analysis provides a basis for inferring potential additive effects at an estimated median chemical intake for the U.S. population. Pathway based approaches currently under development ((NASEM, 2017; Gwinn et al., 2017; Fay et al., 2018; Browne et al., 2017; Kleinstreuer et al., 2018) leverage existing knowledge of biological pathways to integrate mechanistic data from high-throughput assays to in vivo phenotypic outcomes. Individual pathways and networks of pathways (Knapen et al., 2015) integrate data on endocrine bioactivity that converge on specific in vivo outcomes. As increasingly sophisticated predictive pathway models are developed and validated, they can be used in combination to inform chemical evaluations.

Consideration of combined effects of chemicals has often remained within the context of specific regulatory applications (Evans et al., 2016). For example, some proposed frameworks for pesticide combined risk are limited to pesticide active ingredients that act through the same mechanism of action (Bell and Edwards, 2015), but may not account for combined effects of non-pesticide chemicals that also target the same mechanisms or endpoints. Although the magnitude of uncertainty in current high-throughput predictions may make them unsuitable for direct application to regulatory decisions, the high-throughput models (ToxCast, HTTK, ExpoCast) used in the current analysis may provide a consistent platform for pathway-based evaluation of chemicals and includes a diverse set of relatively data-poor chemicals.

Conclusions

The current level of uncertainty in the respective predictive models for endocrine bioactivity, toxicokinetics, and exposure are relatively large. However, the results indicate that combined exposure scenarios sufficient to influence estrogen bioactivity in the general population cannot be ruled out. Future improvements in screening data and refinement of predictive models will be needed to better inform high-throughput combined assessments, beyond the current uses of high throughput and computational methods for screening individual chemicals potential for endocrine disruption. In particular, techniques are needed to address variability in both toxicodynamics (for example, sensitive populations) and exposure (for example, occupational or heavy consumer product use scenarios).