High-Throughput Analysis of Ovarian Cycle Disruption by Mixtures of Aromatase Inhibitors

Background: Combining computational toxicology with ExpoCast exposure estimates and ToxCast™ assay data gives us access to predictions of human health risks stemming from exposures to chemical mixtures. Objectives: We explored, through mathematical modeling and simulations, the size of potential effects of random mixtures of aromatase inhibitors on the dynamics of women's menstrual cycles. Methods: We simulated random exposures to millions of potential mixtures of 86 aromatase inhibitors. A pharmacokinetic model of intake and disposition of the chemicals predicted their internal concentration as a function of time (up to 2 y). A ToxCast™ aromatase assay provided concentration–inhibition relationships for each chemical. The resulting total aromatase inhibition was input to a mathematical model of the hormonal hypothalamus–pituitary–ovarian control of ovulation in women. Results: Above 10% inhibition of estradiol synthesis by aromatase inhibitors, noticeable (eventually reversible) effects on ovulation were predicted. Exposures to individual chemicals never led to such effects. In our best estimate, ∼10% of the combined exposures simulated had mild to catastrophic impacts on ovulation. A lower bound on that figure, obtained using an optimistic exposure scenario, was 0.3%. Conclusions: These results demonstrate the possibility to predict large-scale mixture effects for endocrine disrupters with a predictive toxicology approach that is suitable for high-throughput ranking and risk assessment. The size of the effects predicted is consistent with an increased risk of infertility in women from everyday exposures to our chemical environment. https://doi.org/10.1289/EHP742


Introduction
Concern is growing worldwide over the negative human health and environmental impacts of chemical pollutants that can interfere with the production, metabolism, and action of natural hormones, the so-called endocrine-disrupting chemicals (EDCs). In humans, EDCs have been linked to reproductive disorders (Sweeney et al. 2015), abnormal or delayed development in children (Schug et al. 2015), changes in immune function (Rogers et al. 2013), and cancer (Birnbaum and Fenton 2003). Exposure to mixtures of EDCs may result in effects that can depart from mere summation (Kortenkamp 2007), and human subgroups (e.g., women) may not be sufficiently protected against mixtures of EDCs by current regulatory limits (Kortenkamp 2014).
Each menstrual cycle in women involves hormonal regulation of follicular growth and maturation resulting in ovulation of a single oocyte (Falcone and Hurd 2013). The cycle is controlled by coordinated stimulations and inhibitions along the hypothalamuspituitary-ovarian axis. Gonadotropin-releasing hormone (GnRH), secreted by the hypothalamus, stimulates the secretion of gonadotropins [follicle-stimulating hormone (FSH) and luteinizing hormone (LH)] by the anterior pituitary gland. Those hormones, in turn, regulate the secretion of ovarian hormones, such as estradiol (E2) or progesterone (P4). Exposures to EDCs that interfere directly or indirectly with any of these hormones can eventually induce infertility or other pathological outcomes. Aromatase is critical because it irreversibly converts testosterone to E2 and androstenedione to estrone, maintaining the dynamic balance between androgens and estrogens.
The objective of the present work was to explore predictively the effects of exposure to large-scale (i.e., potentially real-life) mixtures of aromatase inhibitors on the dynamics of menstrual cycling in women. We input exposure estimates from ExpoCast (Wambaugh et al. 2013) and biological effect data from ToxCast™ (Dix et al. 2007) to coupled pharmacokinetic (PK) and ovarian cycle models; this provided a quantitative mechanistic link between exposure to mixtures of EDCs and their potential adverse effects on the menstrual cycle in women. We compared the expected effects of exposures to single EDCs, as is usually considered by risk assessment provisions in different regulations, with estimated effects of cumulative and concurrent exposures.

Workflow Overview
The overall computational workflow is pictured in Figure S1. Briefly, after selecting the chemicals of interest, we sampled millions of random mixtures of chemicals using the exposure estimates provided by ExpoCast (Wambaugh et al. 2013). Both constant and time-varying exposure scenarios of an adult woman were considered. A pharmacokinetic model of intake and disposition was then used to estimate the blood concentration (over 2 y) for each chemical present in each mixture. The resulting aromatase activity inhibition was estimated using the Hill's doseresponse model parameters provided by ToxCast™ (Dix et al. 2007). A mathematical model of the hypothalamus-pituitaryovarian hormonal events [based on a study by Chen and Ward (2014)] was used to predict the levels of E2, P4, and other quantities characterizing the ovarian cycle, for a reference cycle and following exposure to the mixtures generated. Monte Carlo sampling (Bois et al. 2010) was used to propagate uncertainties in exposure, kinetics, and dose-response relationships up to ovarian cycle perturbation.

Databases, Chemical Selection, and Mixture Sampling
ExpoCast (Wambaugh et al. 2013) provides exposure estimates, with measures of uncertainty, for 1,936 chemicals. Those exposure estimates were obtained using far-field, mass-balance human exposure models (USEtox ® and RAIDAR).
In ToxCast™ (11 December 2013 release), the Tox21aromatase-inhibition assay is a cell-based assay measuring CYP19A1 (aromatase) gene activity via a fluorescent protein reporter gene. Chemicals acting on aromatase mRNA synthesis, degradation, or translation, or on aromatase itself, should give positive results in this assay (Chen et al. 2015). For each chemical x assayed, ToxCast™ provides the geometric mean and standard error for the parameter values AC 50,x (half-maximal response), W x (exponent), B x (baseline value) and T x (maximum value) of a Hill function fitted to the concentration-inhibition data (scaled using the positive and negative controls' data): In ToxCast™, 1,102 chemicals are identified as aromatase inhibitors by the Tox21-aromatase-inhibition assay (on MCF-7 human breast cells) with a fitting that "mate basic requirements of Hill model with some minimal confidence in T and B." Among those, 256 chemicals (matching either by CAS number or by chemical name) also had exposure estimates in ExpoCast. However, cytotoxicity has been shown to induce many false positive results in ToxCast™ (Judson et al. 2016). Of the 256 chemicals mentioned above, we kept only the 86 that had an AC 90 (90% of maximal response) for aromatase inhibition lower than their cytotoxicity AC 10 (10% of maximal response) (as measured by the ToxCast™ proliferation decrease assay on T47D human breast cells). The virtual mixtures generated included all of those 86 chemicals.

Exposure Modeling
ExpoCast provided the molecular mass, geometric mean, and lower and upper 95% confidence limits of the exposure rate (milligrams per kilogram per day) for each chemical present in the randomly generated mixtures. For constant exposure modeling over 10 mo, we sampled a rate for each chemical from the corresponding log-normal distribution, but with a standard deviation (SD) scaled by the square root of the number of days of exposure simulated (because a constant exposure rate should be a time average level in that case) and converted it to micromoles per kilogram per minute.
More realistic time-varying exposures over a 2-y period were also modeled [similarly to the report of Bertail et al. (2010)] using exposure windows of random length and intensity. We first sampled the number n of exposure windows for each chemical in the mixture from a scaled exponential distribution (with a rate parameter equal to 5). That yielded on average 145 exposure events over 2 y (median: 100 events, first quartile: 40 events, third quartile: 200 events). The n exposures' start and end times were sampled uniformly over the 2-y period. The exposure rate during each of the n exposure windows was randomly sampled from the log-normal distribution given by ExpoCast, with an SD scaled by the number of days of the exposure window considered (or unscaled if the exposure lasted less than a day).
To obtain a lower bound on the effect of mixtures, we similarly simulated random nonoverlapping exposures to the 86 EDCs selected (i.e., each person was exposed to the 86 chemicals, in random order, at random times, but to only one chemical at a time).

Pharmacokinetics Modeling
For each chemical in each mixture, a one-compartment PK model was used to estimate its internal concentration (in micromoles) at steady state in the case of constant oral exposure or at any point in time in the case of varying oral exposures. Steady-state internal concentrations for chemical x were calculated as: where F x is the bioavailability of x (unitless), E x is its exposure rate (in micromoles per kilogram per minute, sampled as indicated above), and K e,x is its total body clearance rate constant (per minute). A body density of 1 was assumed. For time-varying exposures, internal concentrations were obtained as a function of time by numerical integration of the following differential over a 2-y period (with an initial value set to zero): We used quantitative structure-activity relationships (QSAR) to obtain central estimates of F x and K e,x for each of the 86 aromatase inhibitors considered. The robustness of the prediction was evaluated by examining compounds from the training set similar to the target substances, together with literature data and references. F x central estimates were obtained at several oral dose levels and were linearly interpolated between dose levels as needed. Beyond the dose rates of 0.001 to 10 mg=d, the F x value at the closest bound was used.
To take into account the uncertainty affecting the QSARestimated PK parameters, we randomly sampled F values from beta distributions (naturally bounded between 0 and 1), with parameters calculated such that the distribution modes corresponded to the interpolated value of F x with a coefficient of variation (CV) of 20% (for null F x modes, a and b were set to 1 and 50, yielding a median at 0.01, a first quartile at 0.006, and a third quartile at 0.03, approximately). K e,x values were lognormally sampled with a geometric mean equal to the central estimates obtained by QSAR and a geometric SD corresponding to a factor of 3.

Ovarian Cycle Model
We adapted the menstrual cycle model presented by Chen and Ward (2014). The model describes the inhibitory and stimulatory effects of hormones E2 and P4 on the hypothalamus-pituitary axis in women ( Figure 1). The equations and definitions of all parameters used in the model are given in the Supplemental Material (see "Menstrual cycle model equations"; see also Tables S3 and S4). In the original model, E2 and P4 were assumed to be instantaneously in equilibrium between blood and the ovaries. Instead, we described the kinetics of E2 and P4 using differential equations (Equations 15-19 in the Supplemental Material, "Menstrual cycle model equations"). That modification had practically no impact on the time course of the model variables during a normal cycle [equilibrium between blood and ovaries is fast, as assumed by Chen and Ward (2014)], but it allowed us to coherently integrate the dynamic aspect of estradiol synthesis inhibition by EDCs. The three additional parameters (blood and ovarian volumes, ovarian blood flow) were obtained from the literature (see Table S3).
An additional variable, the ratio of disrupted over basal E2 synthesis rate constants (ED CYP19 ), was introduced to link the internal doses of chemicals in mixtures to aromatase inhibition. ED CYP19 was calculated using Hill's model, parameterized with the chemical-specific values provided by ToxCast™, as a cumulative product of remaining activity for each of the m chemicals of the mixture considered: For constant exposures, C x was set to the steady-state internal concentrations C x,ss . For time-varying exposures, C x was computed by integration as explained above. The parameter B x was set to zero because a positive inhibition with no dosage would not make sense. AC 50,x , W x , and T x values were randomly sampled using the mean and standard error provided by ToxCast™. AC 50,x was sampled from a log-normal distribution (its logarithm is actually the ToxCast™ fitted value). W x and T x were sampled from truncated normal distributions. Truncation was from 0 to 10 for W x (values beyond 10 would be found for some chemicals for which W x is poorly identified, but have no biological meaning). Truncation was from 0 to 100 for letrozole's T x (the positive control). For the other chemicals, truncation was from 0 to 10,000 over T x for letrozole to properly rescale the ToxCast™ T x values between 0 and 100. ED CYP19 was entered as an input to the ovarian cycle model, which was then solved to obtain the time profile of its output variables over 2 y of simulated time. For constant exposures, we computed the square root of the sum of the squared Euclidean distances between a reference E2 concentration (ED CYP19 set at zero) and the perturbed concentrations (at a fixed set of times) as a summary measure of disruption.

Software Used
The ACD/Labs Percepta platform modules ACD/Oral Bioavailability and ACD/PK Explorer were used for the prediction of oral bioavailability (F) and the total body clearance rate constant (K e ), respectively (see Supplemental Material, "PK modules of ACD/ Labs"; see also Table S1 and Figures S2 and S3). GNU MCSim v5.6.5 (www.gnu.org/software/mcsim) (Bois 2009) was used to build the ovarian cycle model. R v3.1.1 (R Development Core Team) with the parallel, deSolve, and EnvStats packages was used for database processing, numerical integration of the models, and graphics.

Estimates of Internal Dose
The relationship between constant exposure rates and steadystate internal concentrations for the 86 EDCs considered indicates that exposures ranged from 10 −8 lmole=kg=d to 10 −3 lmole= Figure 1. Regulatory pathways of the human menstrual cycle as implemented in the model. During the follicular phase (1: germ cells; 2: developing follicle; 3: mature follicle; 4: ovulation; 5: corpus luteum formation; 6: corpus luteum degradation), negative feedback by estradiol (E2) reduces follicle-stimulating hormone (FSH) secretion, leading to the selection of one follicle for ovulation. Gonadotropin-releasing hormone (GnRH) secretion is promoted by E2 and inhibited by progesterone (P4), inducing a luteinizing hormone (LH) peak and consecutive ovulation. E2 is mainly produced by follicles and corpus luteum and P4 by corpus luteum. kg=d and that the resulting steady-state internal concentrations ranged from 10 −13 lM to 10 −3 lM (see Figure S4). The exposure rates and pharmacokinetic parameters were Monte Carlo sampled as described above. For any single EDC, uncertainty is approximately a factor of 10 for exposures and approximately a factor of 1,000 for the resulting internal concentrations. For time-varying exposures, Figure S5 shows an example of a simulated random 2-y time course of internal concentration for lindane. Such profiles were obtained for each chemical in each simulated mixture.

Cycle Model Behavior
Our implementation of the ovarian cycle model proposed by Chen and Ward (2014) correctly reproduces their results. Human data from McLachlan et al. (1990) and Welt et al. (1999) on LH, FSH, E2, and P4 normal cycles are correctly simulated except for McLachlan FSH data, for which the baseline levels are not well matched. There is a large intra-and inter-subject variability in hormonal levels across women in those data sets (see Figure S6). We took the model-simulated normal cycling of E2 as the "reference cycle" in the following. Constant exposure scenarios result, at steady state, in a constant level of aromatase inhibition. In that case, perturbation depends only on that parameter (according to the model assumptions), so the distance between the perturbed and reference cycles is a useful measure of effect (see Figure 2). As aromatase inhibition increases, cycles become increasingly perturbed and exhibit chaotic features (hence the misalignment of the points in Figure 2). At 5% inhibition (95% of normal aromatase activity), cycles are shortened, baseline levels change little, and peak levels either increase or decrease less than proportionally except for LH. Simply put, the regulations dampen the effect of perturbation. At ∼ 10% inhibition, LH peaks disappear after approximately five cycles, and a major bifurcation in cycle patterns occurs: cycles are further shortened, baseline levels are much increased (doubled for E2 and P4, for example, even though E2 synthesis by aromatase is decreased), and peak levels mostly decreased; E2 distance to normal increases up to a maximum (Figure 2). At higher inhibition levels, the cycles increasingly dampen and disappear completely between 30% and 40% inhibition (see Figures S7-S10). Overall, according to this model, having ∼ 10% constant inhibition of aromatase activity in vivo leads to perturbations of the cycle, which is still under control and should be compatible with normal reproductive function. Beyond 10% inhibition, an actual disruption of the system seems to occur.

Effects of Single Chemicals
We first simulated 1,000,000 constant exposures to each of the 86 chemicals considered, taken individually. In that case, despite accounting for uncertainty in exposure levels and dose-response parameters, none induced >1% aromatase inhibition. Hence, none of those chemicals alone was able to induce a significant disruption of the ovarian cycle. Figure 3 places those chemicals on a map with the slope W of the Hill dose-response curve at AC 50 and the logmargin of exposure as coordinates. The margin of exposure was defined as the ratio of the 97.5th percentile of internal concentrations over AC 50 . The log 10 -margins of the chemicals studied ranged from −10 to −1:8, indicating that for all chemicals, the high end of internal exposure concentrations was at most 1% of AC 50 . In that case, Equation 1 shows that the logarithm of aromatase inhibition is approximately equal to the product of W times the log-margin of exposure. The color background of the map codes for the resulting risk index (i.e., the log 10 -inhibitions) and ranges from <−2 (1%) to approximately −80 (10 −78 %), much too low to elicit changes in ovarian cycles such as those in Figure 2. Therefore, no effects can be expected from typical exposures to those chemicals when  . Map of the aromatase inhibitors studied over a plane defined by dose-response slope W and log-margin of exposure (see text). Colors vary linearly with powers of 10 of aromatase inhibition resulting from W and margin of exposure combinations, from red (10 −2 ) to blue (10 −80 ). For the chemicals-numbers correspondence, see Table S1. considered alone. Table 1 gives the list of the 10 chemicals for which individual risk is the highest. Note that letrozole is the reference chemical for the Tox21 aromatase inhibition assay, which is consistent with its high rank. The others are found in therapeutic drugs, agrochemicals, food contaminants, consumer products, and other materials.

Effects of Mixtures of Chemicals
We generated 1,000,000 hypothetical mixtures of the 86 aromatase inhibitors studied and evaluated their global effect on E2 synthesis and the resulting ovarian cycle disruption. Figure 4 shows a histogram of the resulting inhibition levels. Depending on the (random) composition of the mixtures, and given the uncertainties on exposures and effect parameters, responses ranged from 0% to 100% inhibition, but on average were very high. Such inhibition levels may lead to perturbation of the ovarian cycle, according to the model. Real-life exposures to chemicals do not usually occur at constant levels. We simulated time-varying exposures (see Figure  S5) to investigate the resulting effects on aromatase inhibition and ovarian cycle disruption. To define a lower bound on mixture effects, we first simulated fluctuating levels of exposure to EDCs, but without concomitant exposures to them. Interactions can still occur in that case because of storage in the body or because of persistent effects on the ovarian cycle. Figure 5 shows that only 0.3% of the simulated exposures caused >10% average aromatase inhibition. The maximum inhibition found was close to 50%. More realistic exposure scenarios do not prohibit concomitant exposures. In that case (Figure 6), the distribution of simulated time-averaged inhibitions is shifted toward greater effects, and average inhibitions >20% are not uncommon (yet they do not reach the extreme levels observed in Figure 4). Because inhibition changes with time, the distances between normal and perturbed cycles do not follow the pattern shown in Figure 2 (distances can be much larger), and the link between estradiol inhibition and cycle disruption is harder to establish. Examination of the time course of the dominant follicle mass (F) for 1,000 random simulations of mixtures of the 86 chemicals shows that perturbed cycles typically have a baseline shifted up (which may or may not return to normal) and an irregular succession of peaks (corresponding to ovulation) (Figure 7). The 1,000 simulations examined can be classified into four groups. In group 1 (17% of the samples), the cycles are practically normal with no baseline shifts and at most one or two missing ovulations. In group 2 (73% of cases), baseline shifts are always present but without major irregularities in ovulation. Group 3 (7% of cases) has systematic baseline shifts and frequent or prolonged anovulations. Such cycling would clearly impair fertility. In group 4 (3% of cases), disruption is catastrophic or total. Figure S11 shows the corresponding plots for E2 time courses. Judging by these plots, E2 profiles can have a shifted baseline even in normal ovulation profiles, but otherwise, the patterns are rather similar.

Discussion
We linked ToxCast™ data and ExpoCast estimates of exposures to (mixtures of) aromatase inhibitors and estimated their effects on the ovarian cycle in women. To our knowledge, this is the first application of computational toxicology and high-throughput testing to assessment of the combined effects of exposures to a large number of EDCs.
Our approach is predictive, and we had to make many assumptions and simplifications. ToxCast™ and ExpoCast are incomplete, and a full inventory of all the EDCs to which women are exposed is not yet available. Therefore, we were only able to look at a subset of the potential EDCs. We used human exposure estimates reported by Wambaugh et al. (2013), who gave summary statistics for the distribution of exposures to individual chemicals. This information allowed us to take the correspondingly large uncertainty into account via Monte Carlo simulations. However, we cannot differentiate between uncertainty and variability in those exposure estimates, and we cannot identify subgroups of sensitive individuals. We cannot even focus on women, our target population. More sophisticated exposure models (Isaacs et al. 2014) could help in that respect, but they still deal only with single chemicals and provide no data or estimates on coexposures. Depending on age, occupation, socioeconomic status, ethnicity, and health condition, we are exposed to different cocktails of chemicals in our diet, workplace, environment, and so forth (Tornero-Velez et al. 2012). We modeled coexposures by random sampling, either at a constant level or, more realistically, with time-varying exposure profiles and hence with timevarying mixture complexity. That approach is still imperfect, and we had to guess about the distribution of the number of exposure windows, for example. We respected the distribution of population exposure levels documented by ExpoCast, but lacking coexposure information, our estimates might be lower or higher than in reality. Efforts are ongoing to collect relevant data in, for example, the European Total Diet Study (Vin et al. 2014). An analysis of such data (Traoré et al. 2016) shows that among 153 synthetic chemicals studied in seven typical French diets (food associations), three are aromatase inhibitors according to ToxCast™: zearalenone, triadimenol, and lindane. In this regard, mixtures of ≥86 aromatase inhibitors may seem unrealistic, but only food contaminants were studied by Traoré et al. (2016).
We searched for the 86 selected aromatase inhibitors in a database of consumer products marketed in the United States (Gabb and Blake 2016). Briefly, this database was constructed by scraping product information from online retailers and currently contains 53,743 products. Twelve of the 86 aromatase inhibitors were detected in 5,701 products (representing 11% of the products in the database). [It is worth noting that none of the 86 aromatase inhibitors is among the volatile fragrance chemicals detected in consumer products (Steinemann 2015;Steinemann et al. 2011), so aromatase inhibitors are unlikely to be hidden in generic "fragrance" or "flavor" designations on consumer product labels.] Two-way combinations of these chemicals were found in 220 products, and the three-way combination of carminic acid, FD&C blue no. 1, and retinol was found in 3 products. These findings may not seem to indicate a large problem, but it is an incomplete view of combinatorial exposure. Consider that 3,660 of 4,501 makeup products (81%) in the database contain at least one of the aromatase inhibitors evaluated (carminic acid, retinol, and artificial colors are common ingredients in makeup) and that a typical consumer uses several products each day, possibly even several makeup products. This increases the likelihood of combined exposures. In addition, no readily available data address associations for all near-and far-field exposures for the 86 aromatase inhibitors, which are used in industrial or agricultural processes (18), consumer product formulations (12), biocidal applications (38), and pharmaceutical drugs (18). These usage categories are likely to be independent, so focusing on a few known associations would only give partial answers and would underestimate global risk. We are striving for a more extensive picture. To address the potential overestimation of mixture effects when generating purely random associations, we present the results of a very optimistic exposure scenario (with no coexposures at all). This scenario gives a lower bound estimate: 0.3% of exposures would lead to >10% average aromatase inhibition in women.
ToxCast™ aromatase inhibition data were obtained by exposing cells in vitro. We had no easy way to assess the in vitro kinetics of the substances assayed, and reconstruction methods (Armitage et al. 2014) require input data that we did not have. We assumed that the nominal assay concentrations were those actually experienced by the cells and that equivalent extracellular concentrations in vivo would lead to the same levels of aromatase inhibition. That is a typical assumption, but it is not necessarily correct (Coecke et al. 2013). To obtain extracellular concentrations in vivo, we estimated bioavailability and total clearance with QSAR methods and input them in a simple one-compartment PK model with oral exposure only (even though inhalation or dermal exposures might be more relevant). Again, more sophisticated [physiologically based pharmacokinetic (PBPK)] models and additional PK data would give more precise and more accurate predictions (El-Masri et al. 2016;Wambaugh et al. 2015).
On the effect side, we only considered chemicals for which the ToxCast™ dose-response parameters were estimated with reasonable confidence. This was a conservative choice, and additional chemicals would have been included if more relaxed criteria had been chosen (at the cost of lower confidence in the results). A concern with large databases such as ToxCast™ is the quality assurance for the data provided. Aromatase inhibition may not be the most sensitive toxicity end point, for example, or the action may be due to a burst effect of cytotoxicity (Judson et al. 2016). From the original set of 256 aromatase inhibitors for which we had exposure estimates and ToxCast™ data, 170 were excluded based on cytotoxicity. Those were mostly weak inhibitors (data not shown). For the screen, we kept only the 86 substances that had an aromatase inhibition AC 90 lower than their cytotoxicity AC 10 as measured by the ToxCast™ proliferation decrease assay on T47D human breast cells, a cell type similar to the MCF-7 used by the aromatase inhibition assays (Aka and Lin 2012). We preferred that screen to the omnibus z-score criterion, which aggregates cytotoxicity results from different cell types and species. For the remaining 86 substances, we may still have downplayed other types of toxicity. Cytotoxicity, it should be noted, is not a negligible end point, and an evaluation of the cytotoxicity of mixtures would be interesting in its own right.
The evidence provided by ToxCast™ is also not perfectly predictive of in vivo outcomes in humans. For example, the Tox21 aromatase inhibition assay uses the MCF-7 breast cancer cell line, which might not respond as normal ovary cells would. In addition, our model of the ovarian cycle (Chen and Ward 2014) describes only approximately the complex dynamic interactions between ovarian follicular growth and hormonal homeostasis. The hypothalamus and pituitary gland are treated as a single compartment, and the description of the central hormonal controls is simplified. More complex models have been proposed (Hendrix et al. 2014), but they still make many assumptions and do not seem to offer dramatically better performance. In addition, we note that we treated the parameters of the ovarian cycle as constant, when in fact they are affected by both uncertainty and variability in response to EDCs, adding to the tails of the distributions of our results. However, we did not have sufficient information to define statistical distributions for those parameters.
In terms of results, an obvious question is that of the "bad" actors, that is, the chemicals responsible for the predicted effects. The answer is provided by Figure 3 (and partially by the top-ten Table 1) because here, the impact of individual chemicals on aromatase is only conditioned by internal dose and inhibition potency. Figure 3 is a useful prioritization tool. It is rather simple to construct and does not require running the ovarian cycle (the PK model is needed). However, it does not give an answer in terms of magnitude of effect at the subject level. For that, we need the whole-body ovarian cycle model.
One of the consistent features of E2 cycle perturbation that we found is that the baseline (interovulation) levels of E2 tend to increase (to approximately twice the normal level) in response to aromatase inhibition (which implies a lower rate of E2 synthesis by aromatase). That counterintuitive feature of the complex cycle dynamic is induced by central nervous system (CNS) feedback. Beyond a certain level of inhibition, the control of E2 remains in effect (peaks are still observed) but moves the baseline to a higher value. We do not have confirmation that this is the case in women exposed to EDCs, but that would be interesting information and a potential biomarker of effects. We also note that we lack good measures of perturbation for such complex systems. We used visual inspection to classify cycles for 1,000 time-varying exposures ( Figure 7). Perturbation analysis of more cycles would require more sophisticated tools. Finally, many other perturbation pathways exist for ovarian cycle disruption that were not accounted for (e.g., actions mediated by the androgen or estrogen receptor). The simplicity of our model also precludes investigation of synergistic or antagonistic effects that could result from metabolic or toxicodynamic interactions (Cheng and Bois 2011).

Conclusion
High-throughput data collection requires high-throughput analysis, extrapolation, and decision tools if we want to avoid a bottleneck and accumulation of unused data. We developed such a tool, making use of our increasing understanding of toxicity mechanisms.
This exercise in prediction suggests large data gaps. Our knowledge of exposure to actual mixtures is minimal except in a few cases (e.g., tar, tobacco smoke). For the chemicals studied here, quantitative knowledge of their routes of exposure and PK parameters is also lacking. Our knowledge of endocrine disruption mechanisms is still in its infancy. We should check, for example, the inhibition potential of the 86 substances studied here with better tests and with better characterization of the in vitro fate of the chemicals. Nevertheless, we now have a prioritized list and some reason to deepen our investigations. We should also research relevant human biomarkers of exposure and effects for validation of the results and for better actions.
The basic assumption of regulatory practice and most risk assessments is that keeping individual chemicals under control with reasonable safety factors will keep their joint effects at bay. That may not be the case. Some women have various levels of ovarian dysfunction that can be caused by various internal and external factors. Environmental chemical exposures add on to that background (National Research Council 2009;Zeise et al. 2013). We found that even though individual chemicals are likely "safe" as used now, their joint effects, when they exceed a few dozen in number, can lead to severe disruption of the ovarian cycle in women in a sizable number of cases. Obviously, this study does not provide definite proof that some fertility problems can be caused by real-life exposures to EDCs. Nevertheless, risk assessment practice and regulations should start thinking of the problem of mixtures not as an unsolvable one, but as needing a clearly laid out research agenda. In the case in point, the simple graphical map of internal dose and potency we propose could already be used for prioritization. Given the magnitude of the range of the predicted effects on ovarian cycling, while waiting for confirmation of our results, a cautionary attitude should be adopted.