A Computational Model to Predict Rat Ovarian Steroid Secretion from In Vitro Experiments with Endocrine Disruptors

A finely tuned balance between estrogens and androgens controls reproductive functions, and the last step of steroidogenesis plays a key role in maintaining that balance. Environmental toxicants are a serious health concern, and numerous studies have been devoted to studying the effects of endocrine disrupting chemicals (EDCs). The effects of EDCs on steroidogenic enzymes may influence steroid secretion and thus lead to reproductive toxicity. To predict hormonal balance disruption on the basis of data on aromatase activity and mRNA level modulation obtained in vitro on granulosa cells, we developed a mathematical model for the last gonadal steps of the sex steroid synthesis pathway. The model can simulate the ovarian synthesis and secretion of estrone, estradiol, androstenedione, and testosterone, and their response to endocrine disruption. The model is able to predict ovarian sex steroid concentrations under normal estrous cycle in female rat, and ovarian estradiol concentrations in adult female rats exposed to atrazine, bisphenol A, metabolites of methoxychlor or vinclozolin, and letrozole.


Introduction
Humans may be exposed to numerous chemicals that impact endocrine activity, and notably alter androgen/estrogen balance [1]. Among environmental chemicals, atrazine, vinclozolin, methoxychlor, and bisphenol A were found to be of particular concern. Atrazine, a triazine herbicide which has been widely used in agriculture and is persistent in surface water, has been described in several in vitro studies to increase estrogen through elevation of aromatase levels and activity [2,3]. The fungicide vinclozolin has been documented for the anti-androgenic activity of its metabolite M2 in vitro [4] and in vivo [5]. Methoxychlor is an organochlorine pesticide of known estrogenic activities in vitro and in vivo [6]; its metabolite 2, 2-bis-(p-hydroxyphenyl)-1, 1, 1-trichloroethane (HPTE) displays estrogenic, anti-estrogenic, and anti-androgenic capacities in vitro [7]. Bisphenol A, a plasticizer, was clearly defined as an estrogenic agent due to its capacity to bind estrogen receptor with an EC50 in the sub-micromolar range [8]. As far as drugs are concerned, a good example of pharmacologically-designed endocrine modifier may be letrozole [9]. This potent and highly specific nonsteroidal competitive aromatase inhibitor, used for estrogendependent breast cancer, has been characterized by a half maximal inhibitory concentration (IC50) of 7 nM [10].
A potential target for endocrine disrupting chemicals (EDCs) is steroidogenesis. In females, sex steroids are synthesized primarily in the ovaries and derived from cholesterol through a series of biochemical reactions [11]. Among steroidogenic enzymes, cytochrome P450 aromatase (Cyp19), which catalyses the final irreversible conversion of androgens to estrogens in granulosa cells (GCs), appears to be a key target. Aromatase disruption is often associated with EDC toxicity [12], and several assay guidelines recommend testing chemicals for that endpoint [13]. Aromatase expression is regulated by follicle-stimulating hormone (FSH), through multiple signaling pathways including cyclic adenosine monophosphate (cAMP)-dependent regulatory events [14]. In GCs, the final steps of steroidogenesis are also mediated by 17bhydroxysteroid-dehydrogenases (Hsd17b1 and Hsd17b2), which catalyze the conversion of inactive sex steroids to active ones via Hsd17b1 or vice-versa by Hsd17b2 [15].
Assessing EDC toxicity is a challenge, given the complexity of the endocrine system and despite the increasing development of data on its workings. Most standardized ''regulatory'' tests developed to study EDC toxicity involve rats. Those in vivo tests naturally integrate hormone metabolism and feedback loops. They typically look at relevant integrated toxicity endpoints, such as impact on fertility [16]. In vitro models have also been extensively developed: they are faster, cheaper, and they spare animal lives [17]. They help the researcher to elucidate toxic mechanisms in a simple isolated system and, when performed on human cells, they avoid difficult interspecies transpositions.
Both characterization and quantification of toxicity mechanisms are necessary for a reliable quantitative in vitro to in vivo extrapolation (QIVIVE) [18]. In order to improve QIVIVE for endocrine toxicity, we developed and parameterized a dynamic systems biology model of the final steps of steroidogenesis in rat ovaries. We calibrated our mathematical model in a Bayesian framework on the basis of in vitro experimental data obtained from rat granulosa primary cell cultures. For cross-validation, the in vitro model was transposed to an in vivo context and predictions were compared with in vivo hormone dosage data obtained in control animals. We finally used our model to predict the effects of five selected EDCs on gonad estradiol (E 2 ) secretion, based on in vitro data following exposure to atrazine, bisphenol A, methoxychlor metabolite HPTE, vinclozolin metabolite M2, and letrozole. These chemicals were chosen based on their known endocrine activity in vitro and in vivo.

In Vitro Experiments
Rat GC isolation and in vitro culture. Immature (21 days old) Sprague-Dawley female rats (certified virus-free) were purchased from Janvier (Le Genest-Saint-Isle, France). They were housed with a 12 h light and 12 h dark cycle and received food and water ad libitum. All procedures were reviewed and approved by the Institutional Animal Care and Use Committee of INERIS. All animals were 26 days old at the start of treatment. Each animal was injected subcutaneously with diethylstilbestrol (DES; Sigma Aldrich Chemical Co., Saint-Quentin-Fallavier, France) dissolved in corn oil (100 mg/0.1 ml) every day for 3 days to increase the number of GCs. On the third day, the animals were sacrificed by a lethal intraperitoneal pentobarbital injection. Five animals were sacrificed for each experiment. The ovaries were harvested, and the associated fat, oviduct, and bursa ovary removed; the samples were placed in ice-cold medium 199 (M199; Sigma Aldrich), and punctured several times with a 26-gauge needle until the antral follicles ruptured and released the GCs. The GC-rich medium was centrifuged (200 g) for 5 min to obtain a GC pellet, which was resuspended in Dulbecco's modified Eagle medium/Ham's F-12 nutrient mix (DMEM/F-12; Sigma Aldrich) containing 5% fetal bovine serum, 100 mg streptomycin per ml, and 100 IU penicillin per ml. The cells (300,000/ml) were plated into 12-well culture plates (2 ml/well)), and grown at 37uC in a humidified atmosphere with 5% CO 2 . The cells were allowed to attach for 72 h prior to treatment to minimize any effects due to in vivo DES priming [19].
GC treatment. We performed two experimental studies: a baseline (control) study with measurements at 4 h, and an ''EDC study'' with control (0.1% dimethyl sulfoxide, DMSO, in serumfree and phenol red-free culture medium) and four chemicals (atrazine, bisphenol A, HPTE, and vinclozolin M2) at 10 mM in a final concentration of 0.1% DMSO culture medium, with measurements at 4 h. The chemical concentration was chosen on the basis of relevant literature [20]. Cellular viability was determined by trypan blue exclusion staining, visual inspection for morphology, and cellular attachment. mRNA level and direct aromatase activity measurements. mRNA levels and direct aromatase activity were quantified according to previously described methods [20]. Briefly, mRNA was extracted from the cells then reverse transcribed. Target fragments were amplified by real-time polymerase chain reaction. Aromatase enzymatic activity was measured on microsomal fractions of GCs with the tritiated water release assay [21]. These experimental data were expressed as ''fold difference'' between treated and control conditions. Differences of single doses from controls were statistically analyzed with a Mann-Whitney non-parametric test. Differences with a P value of less than 0.05 were considered to be statistically significant.

In Vivo Experiments
The female Sprague-Dawley rats used were approximately 8 weeks old at the start of chemical exposure. Estrous cycle staging was done with vaginal smears collected twice a day and classified microscopically as diestrus, proestrus, estrus, or metestrus [22]. We performed two experimental studies: a baseline (control) study, measuring ovarian steroid concentrations across the estrous cycle, and an ''EDC study'' where each animal in diestrus stage was administered a test chemical or vehicle by gavage (atrazine 200 mg/kg, dissolved in 0.5% methylcellulose; bisphenol A or methoxychlor at 200 mg/kg, dissolved in corn oil; vinclozolin 100 mg/kg, dissolved in corn oil). The animals were sacrificed six hours after treatment; ovaries were harvested, weighed, and homogenized in PBS-buffered water for tissue dosages. Atrazine, bisphenol A, methoxychlor metabolite HPTE, vinclozolin metabolite M2, testosterone (T), androstenedione (A), estrone (E 1 ), and E 2 were detected and quantified in whole ovaries by liquid chromatography with tandem mass spectrometry detection (LC-MS/MS) [23]. Differences between treated and control animals were statistically analyzed with a Mann-Whitney non-parametric test. Differences with a P value of less than 0.05 were considered to be statistically significant.

Model Chemical
We choose to include an additional compound to further test and cross-validate our mathematical model. Letrozole appeared to be a very good choice, in the sense that it is pharmacologically designed to specifically inhibit aromatase, which is one of the main target described in our computational model. This compound was not tested on our in vitro and in vivo systems, but experimental data were gathered from the literature [10,24].

Mathematical Model Development
Model overview. The model describes the final metabolic and transport steps of the steroidogenesis pathways in rat GCs (Figures 1 and 2). Metabolic steps include synthesis and degradation of Cyp19, Hsd17b1, and Hsd17b2 mRNAs and proteins, conversion of A into T, E 1 , and E 2 , and modulation of steroidogenic enzyme expression by FSH or an EDC. In vitro, transport includes GC uptake and secretion of A, T, E 1 , and E 2 . In vivo, transport also includes entry of A and T in ovaries, and exchange of hormones between extracellular space, GCs, and other kinds of cells ( Figure 2).
Metabolic reactions. mRNA and protein synthesis. Cyp19 and Hsd17b1 mRNA quantities in GCs (e mRNA in pg/cell, with e = Cyp19, Hsd17b1, or Hsd17b2) depend on their synthesis with baseline rate n mRNA, e (pg/min). This rate is eventually altered by an EDC X (inducing fold-change f X (unitless)), upregulated by FSH (pg/cell) (with slope factors k (/pg FSH), and affected by experimental variability (due to differences in cell pre-treatment, modeled by a variability factor s L (arbitrary unit)); mRNA levels depend also on their degradation, with rate constant d mRNA (/min):

LCyp19 mRNA
Lt~n mRNA:Cyp19 |s L |f X ,Cyp19,t : LHsd17b1 mRNA Lt~n mRNA:Hsd17b1 |s L |f X ,Hsd17b1,t : Fold-change for species e mRNA was obtained from experimentally measured mRNA levels and computed as: In contrast, the Hsd17b2 mRNA quantity in GCs is not assumed to be strongly controlled by FSH [25] nor affected by EDCs, and the corresponding equation is simply: For the three enzymes e (in pg/cell), the following mass-balance equation, with synthesis rate constant n prot,e (/min) and degradation rate constant d prot (/min), applies: Our experiments on GCs [20] gave us Hsd17b1 and Hsd17b2 initial mRNA and protein quantities, relative to Cyp19. We translated them to absolute values (pg/cell) on the basis of the initial quantities of Cyp19 mRNA and protein in GCs obtained from the literature ( Table 1). We assumed that these values were steady-state values, in the absence of FSH stimulation, EDC alteration, or experimental variability. Values for the mRNA and protein degradation rate constants (d mRNA,e and d prot,e ) were found in the literature ( Table 2). Using the above steady-state assumption, we set equations 1, 2, and 4 for mRNA quantities, and equation 5 for protein quantities, equal to zero and rearranged them for n mRNA.e Similarly, assuming that one mRNA gets translated into one protein, the value of n prot.e was computed for the three enzymes e as: Aromatase protein synthesis (/min) u prot.Cyp19 6000 b Hsd17b1 protein synthesis (/min) u prot.Hsd17b1 6300 b Hsd17b2 protein synthesis (/min) u prot.Hsd17b2 6000 b Maximal reaction rates V max (pmoles/min/pg enzyme) Hsd17b2, T R A reaction l Hsd17b2,T 6.65610 28c Hsd17b2, E 2 R E 1 reaction l Hsd17b2,E2 7.91610 28c

26c
A extra-over intra-cellular partition coefficient (unitless) Steroid biotransformation. The relevant enzymatic reactions in GCs, catalyzed by Cyp19, Hsd17b1, and Hsd17b2, were modeled by the following competitive Michaelis-Menten metabolic terms a i , where l e,Z (pmoles/min/pg enzyme) and j e,Z (pmoles) denote respectively V max and K m parameters for enzyme e and substrate Z (A, T, E 1 , or E 2 ).
Methoxychlor metabolite HPTE inhibits aromatase activity directly and competitively [20]. To model that effect, the parameter f M in equations 9 and 12 below represents the foldchange of the aromatase K m for its substrate Z (j Cyp19,Z ), observed in vitro. Since aromatase activity is inversely proportional to its K m , this fold-change f M corresponds to the inverse of fold-change for aromatase enzymatic activity between treated and control cells. Fold-change for K m parameters j Cyp19,A and j Cyp19,T corresponds to: f M,j Cyp19,Z ,t~M easured aromatase activity in control cells t Measured aromatase activity in treated cells t ð8Þ The conversion of A into E 1 by aromatase takes into account T competition for the enzyme (the steroids are subscripted with GC, denoting the intra-cellular quantities): Conversion of A into T by Hsd17b1, with E 1 competition: Conversion of T into A by Hsd17b2, with E 2 competition: Conversion of T into E 2 by aromatase, with A competition: Conversion of E 1 into E 2 by Hsd17b1, with A competition: Hsd17b1,E 1 : Conversion of E 2 into E 1 by Hsd17b2, with T competition: Hsd17b2,E 2 : In order to model the isotopic measurement of tritiated water (T 2 O) production during the conversion of tritiated A to E 1 (see in vitro experimental section), we need the formation rate of T 2 O, which is simply: The parameters of the above equations are listed in Table 2.
Transport kinetics. The model was first developed to simulate in vitro conditions, and then adapted to model in vivo conditions. While the GC internal workings remained the same, different exchanges with the environment had to be described ( Figure 2).
Transport kinetics in vitro. The in vitro model is divided in two compartments: GCs and culture medium ( Figure 2A). For A (pmoles), T (pmoles), E 1 (pmoles), E 2 (pmoles), and FSH (pg), simple diffusion kinetics were assumed. The hormone quantity in a GC (X GC ) has a rate of change equal to: where K in,X (ml/min) is the rate of medium (''med'') uptake by the GC, K out,X (ml/min) the rate of excretion by the GC, X med (pmoles or pg) the hormone quantity for one GC in the medium (total quantity divided by the number of cells used in a given assay), V med (ml) the volume of culture medium for one GC (total volume divided by the number of cells), and V GC (ml) the volume of one GC. K in,X was computed by dividing K out,X by the extra-over intracellular partition coefficient R oi,X (unitless), given in the literature [26]: Conversely, the hormone quantity for one cell in the medium (X med ) has a rate of change equal to: The cellular kinetics of A, T, E 1 , and E 2 quantities depend on the entry in and exit from the cell and on their metabolism by Cyp19, Hsd17b1, or Hsd17b2: where the a i are the Michaelis-Menten metabolic terms described in equations 9 to 14. The diffusion and transport of T 2 O in the in vitro system was not modeled, as the total quantity of T 2 O formed was directly measured. Transport kinetics in vivo. For in vivo simulations ( Figure 2B), the ovary was subdivided into three compartments: GCs, thecal and interstitial cells (''others''), and extracellular/vascular space (''ext''). The transport kinetics of hormone X in each cellular compartment depend on entry rate constant (K in ) and exit rate constant (K out ) for a cell, on the hormone concentrations in each cell, and on the number of cells (N GCs or N others ). The differential equations for the ''GC'' and the ''other cell'' compartments are: where V ext and V others are the volumes of the extracellular and ''other cell'' compartments, respectively. The differential equation for quantity X ext in the extracellular compartment is: where Q input,X (pmoles or pg/min) is the rate of input of hormone X in the ovary (coming from blood), F ov (ml/min) the efflux of X from the ovary (clearance by blood flow), and V ov,diestrus the ovarian volume at diestrus (which was set at 0.05 ml [27]). For mimicking the female estrus cycle in vivo, Q input,X for FSH and androgens were modeled as cyclic forcing functions, which were adjusted to give ovarian concentrations matching our in vivo physiological observations (see Figure 3). Q input,X is determined as: where Q base,X (pmoles or pg/min) is the constant baseline concentration of hormone X, Q scale,X (unitless) the constant scale for hormone X magnitude, and Q shape,X (pmoles or pg/min) the variable magnitude of hormone X (adjusted to match the known hormone concentrations). The time courses of N GCs , V ext , and V others during the estrous cycle were also modeled by forcing functions. The intracellular kinetic equations of the various hormones were the same as in the in vitro model (see metabolic reaction section).
Parameter value assignment and model calibration. Whenever possible, the model parameters were set to meaningful and physiologically based values that we directly measured in vitro or that we found in the published literature ( Table 2). The remaining model parameters (Table 3) were calibrated using in vitro experimental data that we developed ourselves (see above, in vitro data section), or that were published in the literature (Information S1). A Bayesian numerical approach, Markov Chain Monte Carlo (MCMC) simulations [28], was used.
The published in vitro data we used to calibrate the model included different cell pre-treatment protocols, which induced a large inter-study variability in baseline transcription rates n mRNA.e . That random effect was modeled with a variability factor s L (see equations 1, 2, and 4), assumed to be log-normally distributed around a mean m s , with variance S 1 . The hyperparameters m s and S 1 were in turn assigned vague prior distributions ( Table 3). The individual random effects s L (one per data set used, see Information S1), m s , and S 1 were calibrated together with the other parameters.  The other parameters to be calibrated were assigned a prior distribution (Table 3). We mostly used lognormal distributions with geometric means set at physiologically relevant values. The geometric standard deviations were set to 2 or 1.2 for the parameters for which we had better information ( Table 3). The data likelihoods were assumed to follow a lognormal distribution around the model predictions, a standard assumption with such measurements. The measurement error variances, which were assumed to be different between mRNA/protein quantities (S 2 ) and hormone measurements (S 3 ) (Table 3), were calibrated together with the other (physiological) parameters. A total of 24 parameters (11 physiological and 13 statistical) were MCMC sampled.
MCMC simulations (Metropolis-Hastings algorithm) were performed in triplicate chains of 20,000 iterations. For each model parameter sampled, convergence was evaluated using the last 10,000 iterations from each chain and the potential scale reduction criterion R R of Gelman and Rubin [29].

Flux Analyses of In Vitro and In Vivo Experiments
Maximum a posteriori probability estimates of the calibrated parameters (Table 4) were used to do metabolic flux analyses [30], computing the rate of each steroid biotransformation reaction (a 1 to a 6 , equations 9 to 14) as a function of time, to determine the predominant reactions for the conversion of A to E 2 .

Model Cross-validation Using In Vivo Data
In order to evaluate the predictive capacities of the model, we used random parameter vectors from their joint posterior distributions obtained by calibration with in vitro data (Table 4), and some other parameter distributions (Table 5), to simulate in vivo conditions. Table 5 includes parameters which were not calibrated from in vitro data because reasonable values were obtained for them in the literature, but which nonetheless have in vivo variability. We then simply compared the model predictions to the corresponding in vivo data. The cyclic entries of androgens and FSH in GCs and the time-varying number of ovarian cells were modeled as described in section ''Transport kinetics in vivo''.

Predictive Simulations of Endocrine Disruption
To evaluate the capacity of the above model to predict in vivo effects of EDCs on E 2 secretion on the basis of in vitro data, we ran a series of simulations of endocrine disruption by atrazine, bisphenol A, methoxychlor metabolite HPTE, vinclozolin metabolite M2, and letrozole over two estrous cycles. The mRNA and K m fold-changes f X (equations 1-3, 8, 9, and 12) were changed to their experimentally observed values (see Table 6), starting eight hours after the beginning of the second modeled diestrus. We then compared the in vivo E 2 quantities measured experimentally in EDC-treated females in diestrus with the model predictions. The hypothesis that the distributions of experimental data and model predictions were identical was statistically tested with a two-sample Kolmogorov-Smirnov test [31]. Differences with a P value of less than 0.05 were considered to be statistically significant.

Software Used
Cell Designer 4.2 [32] was used to produce Figure 1. Model simulations, MCMC simulations for model calibration, and flux analyses were performed with GNU MCSim v5.4.0 [28]. Statistical analyses and plots were performed with R, version 2.14.0 [33].

In Vitro Experimental Results
To evaluate and quantify how the selected EDCs affect aromatase and Hsd17b mRNA levels, as well as aromatase function, we exposed rat primary GCs (or microsomal fractions for direct aromatase activity) to atrazine, bisphenol A, methoxychlor metabolite HPTE, or vinclozolin metabolite M2. The chemical concentration used corresponded to the highest one found in rat ovaries following oral exposure to a high dose of each selected EDC [23]. None of the chemicals tested affected cell viability, as assessed with trypan blue exclusion staining and morphological evaluation. The purpose for measuring aromatase activity on microsomes (rather than in entire cells) was to discriminate a direct effect of chemicals at the functional protein level from an effect due to altered protein levels. Table 6 illustrates the fold-changes (relative to appropriate controls) for aromatase direct enzymatic activity, and aromatase and Hsd17b1 mRNA level modulation. In our experiments, Hsd17b2 mRNA levels were too low to be quantified. Atrazine, bisphenol A, and vinclozolin metabolite M2 did not affect aromatase direct activity, whereas HPTE decreased it by 11%.
Atrazine increased aromatase and Hsd17b1 mRNA levels with fold-inductions of 1.94 and 3.04, respectively. Bisphenol A increased 1.61-fold the amount of aromatase mRNA levels, but did not modify the Hsd17b1 mRNA levels. HPTE did not affect aromatase or Hsd17b1 mRNA levels. Vinclozolin up-regulated 3.13 and 1.61-fold aromatase and Hsd17b1 mRNA levels, respectively.
Experimental in vitro data for letrozole were found in the literature [10]. Hence, at the concentration of 50 nM, which corresponds to that found in rat ovaries after treatment with letrozole at 5 mg/kg, aromatase activity was decreased to 29% compared to control.

In Vivo Experimental Results: Baseline Study
Gonadal sex steroid and blood FSH concentrations in healthy cycling female rats were measured at several times, falling in different periods of the estrous cycle ( Figure 3). Those results are well in agreement with the published scientific literature [34].

Model Calibration
Twenty-four model parameters were jointly calibrated using MCMC simulations. The three chain simulations converged after 10,000 iterations ( R R was at most 1.01 for all sampled parameters). The posterior fit after parameter Bayesian calibration has an average absolute deviation of 18.82% between data and predictions. Table 4 presents summary statistics of the posterior distributions for the parameters calibrated. Those statistics are based on 30,000 iterations (the last 10,000 iterations from each of the three chains). For FSH effect on aromatase and Hsd17b1, while prior distributions were quite vague (see Table 3), posterior distributions indicated that the effect of FSH on aromatase is about four times higher than its effect on Hsd17b1. V max and K m for aromatase had informative priors and the posterior distributions were close to them; this probably may confirm the prior knowledge. However, although we used the same aromatase V max prior for A and T, posterior distributions revealed a 3-fold higher V max for T. On the contrary, according to posterior distributions, the aromatase K m for A is 5-fold smaller than the one for T. The posterior distributions for Hsd17b1 V max and K m were modified by a factor 1 to 2 for j Hsd17b1,E1 , j Hsd17b1,A , and l Hsd17b1,A , and by a factor 3 for l Hsd17b1,E1 . The average inter-study variability factor was about 40%. Study-specific variability factors ranged from 0.03 to about 5. The measurement error variances corresponded to a coefficient of variation of about 65% for mRNA/protein quantities (S 2 ), and 48% for hormone measurements (S 3 ). Figure 4A, B shows the results of A, E 1 , T, and E 2 in vitro interconversion flux analysis 48 h after addition of the substrate A (200 nM), with or without FSH (20 ng/ml). The flux value for the reference reaction A to E 1 increased from 7.29610 29 pmoles/ min/cell (without FSH) to 8.72610 28 pmoles/min/cell (with FSH). The other reaction relative values show that the preferential pathway for E 2 synthesis in GCs in vitro is conversion of A into E 1 , which is then converted into E 2 , both with or without FSH. Figure 4C, D, E shows the results of steroid hormone in vivo interconversion flux analysis at different times of the estrous cycle. The flux value for the reference reaction A to E 1 increased from 5.10610 29 pmoles/min/cell at the estrus stage to 6.09610 29 pmoles/min/cell at the diestrus stage and to 6.17610 29 at the proestrus stage. Those in vivo results, which show a preferential pathway for E 2 synthesis trough E 1 conversion, itself coming from A, are in accordance with the in vitro ones. Some differences between in vitro and in vivo flux analyses can be noted, like the greater conversion of T to E 2 in vivo than in vitro, or the relative importance of the retroconversions of T and E 2 to A and E 1 , respectively, in in vitro experiments compared to in vivo ones.

In Vivo Model Simulation
In order to evaluate the model accuracy, we set the GC model parameters in vivo to the values found by calibration with in vitro data. In vivo parameter uncertainty and variability were modeled by distributions of hormone inputs, clearances, mRNA/protein degradation and specific synthesis, and Hsd17b2 apparent kinetic constant parameters (Table 5). These distributions, used in inputs to Monte Carlo simulations, yielded predictive confidence intervals.
We compared the model-predicted ovarian steroid concentrations with the data from baseline experiments (Figure 3). The mean of our model predictions were within the 95% confidence interval of the model predictions. A quantitatively close profile for predicted data and experimental data was observed for E 2 , whereas the values for E 1 in the diestrus stage were somewhat under experimental data. Profiles for FSH and androgens are shown for informative purpose, since they were constructed (using forcing functions) to match the observed profiles.
In Vivo Experimental Results: EDC Study Figure 5 illustrates the distribution of experimentally measured ovarian E 2 levels following EDC oral exposure. E 2 levels were significantly increased in atrazine-treated females, whereas no statistically significant alteration of E 2 was observed with bisphenol A, methoxychlor, and vinclozolin treatment. As far as vinclozolintreated rats are concerned, one of those showed an elevated E 2 ovarian concentration.
In vivo data extracted from the literature showed a significant decrease of E 2 in letrozole-treated rat ovaries, compared to control [24].

Predictions for Ovarian Estradiol Concentrations in EDCtreated Female Rats
In vitro results with atrazine, bisphenol A, methoxychlor metabolite HPTE, and vinclozolin metabolite M2 showed a modulation in mRNA levels after four hours of chemical exposure; and only cells treated with HPTE and letrozole showed a significant decrease in aromatase enzymatic activity (Table 6). To further evaluate the predictive capacity of the model, we simulated E 2 concentrations in female rats exposed to atrazine, bisphenol A, methoxychlor, vinclozolin, and letrozole for six hours. After ''in vivo'' simulation with the mathematical model, we compared E 2 values predicted with those experimentally measured ( Figure 5). A two sample Kolmogorov-Smirnov test was performed for each pair of data (experimental versus predicted data for each treatment). It confirmed that the distributions of experimental data and model predictions were similar for control, atrazine, bisphenol A, methoxychlor, and letrozole treatments. Significantly different distributions were found only for vinclozolin treatment (p = 0.021).

Discussion
The model presented here offers a detailed description of some steroidogenic processes, focusing on what we felt to be the most important ones for in vitro to in vivo extrapolation. The Bayesian approach used for calibrating the model parameters permitted us to take into account both uncertainty and variability in experimental data, which is an asset for the relevance of the predictions. The in vitro and in vivo data we generated allowed us to finely calibrate and cross-validate the model, which was able to quantitatively predict E 2 ovarian concentration in physiological conditions or after exposure to selected EDCs. This model, in spite of its limitations, has many potential mechanistic or predictive applications, as we discuss in the following.

Model Development
In the context of EDC toxicity assessment, some authors developed systems biology models of the hypothalamic-pituitarygonadal (HPG) axis. Many of them are graphical systems models, which allow researchers to visualize and think more clearly about the impact of chemicals on the HPG axis (as reviewed by [35,36,37]. They can also provide a framework for integration of quantitative computational models, such as those of Breen et al. [38], Watanabe et al. [39], and Li et al. [40]. Breen et al. [38] proposed a steady-state model of fathead minnow ovarian synthesis and release of T and E 2 ; Watanabe et al.'s model [39] simulates synthesis and feedback loops for T, 11-ketotestosterone, E 2 , and vitellogenin plasma concentrations in male fathead minnow; Li et al.'s model [40] simulate E 2 , T, and vitellogenin plasma concentrations in female fathead minnow. Those models focused on fish as a target species since endocrine disruption is well documented in aquatic species [41]. However, the assessment of EDC toxicity for humans warrants the development of mammalian models. We chose to develop a computational model focusing on the last steps of steroidogenesis in the rat ovary. This choice seemed to be a good compromise between our purpose (make quantitative in vivo predictions for a mammal based on in vitro measurements), and the data available to calibrate and crossvalidate our model.

Model Calibration
The calibration of the model was done on the basis of several in vitro data sets, including our own. The diversity of protocols, in particular for cell pre-treatment, led us to model inter-study variability. Experiments reported in the literature were done to compare treatments with control conditions rather than to develop a computational model. For that reason, they lack endpoints such as time-response curves at several FSH levels, precursor hormone measurements, etc. In that sense, to develop a quantitative computational model forces one to identify the kind of data needed. Beyond answering the questions raised when developing the model, such a refinement of experimental design may yield new findings about cellular biology and toxicology in vitro. In any case, the model was able to account for the differences between studies and predicted the endpoints reasonably well. That can be considered as the first part of our model validation process.

Model Evaluation with Cross-validation
Results from a baseline in vivo study, without EDCs, were used to evaluate our model ability to predict some features of steroid synthesis in normal physiological conditions. We used model parameter values estimated by calibration of the in vitro data ''as is'', without adjustment, to simulate E 1 and E 2 production by the ovary in vivo. The results showed that the model was able to accurately simulate ovarian E 2 concentrations during normal cycling in female rats. The results for E 1 were less convincing, in particular during diestrus. We did not go as far as to model the ovarian steroid output, plasma concentrations, and the hypothalamic-pituitary (HP) feedback. That was beyond the scope of our work, but more importantly, modeling steroid output from the ovaries and sex steroid plasma concentrations would have required calibration of several more parameters and compromised the mere feasibility of model cross-validation.

Model Predictions and Biological Insight in Baseline Conditions
Updating the a priori parameter distributions into posteriors gives us some insight into features of the rat sex steroid synthesis network. For example, the preferred conversion of A into E 1 by aromatase (in spite of its conversion into T by Hsd17b1) seems due to differences in K m values of androstedione for aromatase and Hsd17b1, rather than to differences in V max values.
The flux analyses indicate that the preferential pathway for E 2 synthesis involves E 1 both in vitro and in vivo. They also point out Flux analyses show clear differences between in vitro and in vivo conditions. For example, steroid inactivation reaction fluxes (T to A and E 2 to E 1 ) are ranged from 10 23 to 10 26 pmoles/min/cell in vitro, and ranged from 10 24 to 10 212 pmoles/min/cell in in vivo conditions. Those differences can be explained by differences in hormone inputs to the system. Fluxes depend on reaction parameter values and hormone inputs applied. We showed that keeping parameter values equal in vitro and in vivo, and simply changing hormone inputs, is enough to explain flux differences between in vitro and in vivo conditions.

Model Predictive Capacity Evaluation with Selected EDCs
To further evaluate the model predictive capacity, we simulated in vivo steroid concentrations in the ovaries after chemical exposure and compared them to original experimental results. Simulations were performed by modifying aromatase K m or mRNA levels on the basis of transcriptomic and enzymatic activity data obtained in vitro for GCs. We limited our predictions to six hours postexposure, a period during which feedback regulation can be assumed to be negligible.
Results show that our model predictive capacity was different according to treatment. Model predictions were found to follow the same distributions as the experimental data, except for vinclozolin. However, Figure 5 shows nuances between treatments. The model predicted reasonably well the early ovarian response in E 2 concentration for adult female rats exposed to atrazine and letrozole. Atrazine and letrozole mechanisms of action can explain why their effects were the most clearly seen experimentally and the best predicted by the model after a few hours. Indeed, we have previously shown [27] that elevated aromatase mRNA expression (see also Table 6) and the subsequent increase in aromatase catalytic activity in atrazinetreated females explain a large part of the increase in estrogen levels. As far as letrozole is concerned, it was designed to be a specific aromatase inhibitor. The early ovarian responses in E 2 concentration for adult female rats exposed to bisphenol A, methoxychlor, or vinclozolin were less well predicted. The effects of bisphenol A, HPTE, or vinclozolin M2 on aromatase or Hsd17b1 did not explain the in vivo modulation of estrogen levels following treatment, although they can significantly affect enzyme mRNA levels in vitro. Instead we previously hypothesized that the main mechanisms of action are: a disruption of the hypothalamicpituitary-adrenal axis for methoxychlor and vinclozolin; a peripheral effect on conjugation/deconjugation metabolism processes for bisphenol A [27]. The model, which doesn't predict very well variations of E 2 concentrations following exposure to those three chemicals, may confirm that the effects on granulosa steroidogenesis are not predominant. Furthermore, vinclozolin predictions were less precise, and showed higher variability. That is actually an interesting feature: vinclozolin mechanism of action is known to be more complex, acting notably by its antiandrogenic metabolite M2 [4], and subject to variable amplifica- tion in the steroidogenesis pathway. The experimental data themselves showed higher variability for vinclozolin, although the small number of animals tested precludes strong conclusions.
Even if predictions for E 2 levels compared well with experimental values, the usefulness of the model could be improved. First, it does not account for EDC effects on androgen precursors, and can only predict effects for chemicals that act on the last steps of steroidogenesis. An improvement would be to add other pathways to the mathematical model, such as steroidogenic processes in thecal cells. The model may also integrate effects on steroid receptors, like the estrogen one, which is the target of numerous chemicals [42]. The model also lacks numerous feedbacks, in particular those mediated by the HP axis. Thereby, for now, the model predictions for steroid ovarian concentrations are of limited value for a complete analysis of endocrine disruption. Rat HP axis feedback models previously described [43,44] might be useful for coupling with ours.

Model Potential
Despite the limitations discussed above, the model perspectives are multiple. All the reaction parameters can be modulated to reflect changes observed in vitro, for example. That approach can be very useful for investigating mixture and chronic effects. It can also help formulate hypotheses and design experiments aimed at understanding the mechanisms of endocrine toxicity, notably for the effects which follow a non-monotonic dose-response, like those of EDCs. A model integrating feedback regulations would permit to describe further targets, such as the HPG axis, enzyme inhibition, or local gene expression effects.
Observations of alterations in ovarian functions at molecular and biochemical levels are useful for regulatory decision-making only if these changes can be translated into effects at higher biological levels of organization. The model is able to make quantitative predictions about steroid secretion based on data on the impact of chemicals on the last steps of ovarian steroidogenesis. Sex steroid concentration changes, even of low scale, account for a large part of effects in reproductive toxicology, but it is not sufficient. Integrated models, predicting multiple endpoints relevant for reproductive toxicology assessment, have been developed in the fathead minnow [39,40]. Since links between sex steroid concentration changes and reproductive toxicity are not clear in mammals, some work still has to be done.

Conclusions
The model developed was able to predict a very sensitive and integrative reproductive endpoint: ovarian sex steroid levels, from in vitro data. The results of flux analyses and predictions of EDCtreated females show that the model not only fits the data empirically, but also captures major features of the GC steroidogenesis network. We carefully limited the scope of our model to ovarian secretion in order to be able to cross-validate it with the data available. In some cases, investigating effects simply on gonads can be a powerful tool for understanding whole-body hormone disruption, in which case the model might be a valuable tool for toxicity assessment. While the predictive capacity of this mathematical model is still limited, it already has potential applications for improved evaluation of endocrine disruption following chemical exposure, in particular for low levels and mixtures of pollutants.

Supporting Information
Information S1 In vitro data used for parameter estimation. The table presents the endpoints experimentally measured by authors that allowed us to update our prior information. Data obtained are both from non-treated cells or FSH-induced conditions. All conditions and measurements were scaled down to one cell by dividing by the number of cells introduced in the assay system. (DOC)