Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

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

  • Nadia Quignot ,

    nadia.quignot@utc.fr

    Affiliations Chair of Mathematical Modeling for Systems Toxicology, Bioengineering Department, Université de Technologie de Compiègne, Compiègne, France, Models for Ecotoxicology and Toxicology, Institut National de l’Environnement industriel et des Risques, Verneuil-en-Halatte, France

  • Frédéric Y. Bois

    Affiliations Chair of Mathematical Modeling for Systems Toxicology, Bioengineering Department, Université de Technologie de Compiègne, Compiègne, France, Models for Ecotoxicology and Toxicology, Institut National de l’Environnement industriel et des Risques, Verneuil-en-Halatte, France

Abstract

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 estrogen-dependent 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 17β-hydroxysteroid-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 (E2) 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.

Materials and Methods

Test Chemicals

Atrazine (CAS number 1912-24-9, purity 97.1%) was provided by TCI Europe (Zwijndrecht, Belgium); methoxychlor (CAS number 72-43-5, purity >95%), HPTE (CAS number 2971-36-0, purity 97%), and bisphenol A (CAS number 80-05-7, purity 99%) were purchased from Sigma Aldrich Chemical Co. (Saint-Quentin-Fallavier, France); vinclozolin (CAS number 50471-44-8, purity 99.5%) was from Greyhound Chromatography (Birkenhead, UK); vinclozolin M2 (CAS 83792-61-4, purity >98%) was from Interchim (Montluçon, France).

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 µg 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 37°C in a humidified atmosphere with 5% CO2. 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 serum-free and phenol red-free culture medium) and four chemicals (atrazine, bisphenol A, HPTE, and vinclozolin M2) at 10 µM 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 (E1), and E2 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, E1, and E2, and modulation of steroidogenic enzyme expression by FSH or an EDC. In vitro, transport includes GC uptake and secretion of A, T, E1, and E2. 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).

thumbnail
Figure 1. Overview of the computational model for steroidogenesis last metabolic steps in a rat granulosa cell.

The transcription and translation events for the three last major enzymes involved in estradiol synthesis, and sex steroid synthesis itself, are modeled, with relevant FSH control, endocrine disrupting chemical (EDC) modulation, or methoxychlor (MXC) aromatase competitive inhibition. Steroids can be transported in and out of cell. In vitro, the exterior compartment corresponds to the culture medium; in vivo it corresponds to the ovary tissue (see Figure 2). Aliases (repeated species labels) are used for clarity but correspond in fact to a unique species.

https://doi.org/10.1371/journal.pone.0053891.g001

thumbnail
Figure 2. Overview of the compartments used to model in vitro (A) or in vivo (B) hormone transports.

In vitro (A), the exterior compartment corresponds to the culture medium. In vivo (B), the ovary tissue is subdivided into three compartments: granulosa cells, “other cells” for thecal and interstitial cells, and extracellular space.

https://doi.org/10.1371/journal.pone.0053891.g002

Metabolic reactions.

mRNA and protein synthesis. Cyp19 and Hsd17b1 mRNA quantities in GCs (εmRNA in pg/cell, with ε = Cyp19, Hsd17b1, or Hsd17b2) depend on their synthesis with baseline rate νmRNA,ε (pg/min). This rate is eventually altered by an EDC X (inducing fold-change fX (unitless)), upregulated by FSH (pg/cell) (with slope factors κ (/pg FSH), and affected by experimental variability (due to differences in cell pre-treatment, modeled by a variability factor σL (arbitrary unit)); mRNA levels depend also on their degradation, with rate constant δmRNA (/min):(1)(2)

Fold-change for species εmRNA was obtained from experimentally measured mRNA levels and computed as:(3)

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:(4)

For the three enzymes ε (in pg/cell), the following mass-balance equation, with synthesis rate constant νprot,ε (/min) and degradation rate constant δprot (/min), applies:(5)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 (δmRNA,ε and δprot,ε) 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 νmRNA.ε and νprot.ε. The value of νmRNA.ε was computed for the three enzymes ε as:

thumbnail
Table 1. Granulosa cell specific mRNA and protein initial values used.

thumbnail
Table 2. Model parameter values (for one cell) obtained from direct measurements on granulosa cells in vitro or from the published literature values.

https://doi.org/10.1371/journal.pone.0053891.t002

(6)Similarly, assuming that one mRNA gets translated into one protein, the value of νprot.ε was computed for the three enzymes ε as:(7)

Steroid biotransformation

The relevant enzymatic reactions in GCs, catalyzed by Cyp19, Hsd17b1, and Hsd17b2, were modeled by the following competitive Michaelis-Menten metabolic terms αi, where λε,Z (pmoles/min/pg enzyme) and ξε,Z (pmoles) denote respectively Vmax and Km parameters for enzyme ε and substrate Z (A, T, E1, or E2).

Methoxychlor metabolite HPTE inhibits aromatase activity directly and competitively [20]. To model that effect, the parameter fM in equations 9 and 12 below represents the fold-change of the aromatase Km for its substrate ZCyp19,Z), observed in vitro. Since aromatase activity is inversely proportional to its Km, this fold-change fM corresponds to the inverse of fold-change for aromatase enzymatic activity between treated and control cells. Fold-change for Km parameters ξCyp19,A and ξCyp19,T corresponds to:(8)

The conversion of A into E1 by aromatase takes into account T competition for the enzyme (the steroids are subscripted with GC, denoting the intra-cellular quantities):(9)

Conversion of A into T by Hsd17b1, with E1 competition:(10)

Conversion of T into A by Hsd17b2, with E2 competition:(11)

Conversion of T into E2 by aromatase, with A competition:(12)

Conversion of E1 into E2 by Hsd17b1, with A competition:(13)

Conversion of E2 into E1 by Hsd17b2, with T competition:(14)

In order to model the isotopic measurement of tritiated water (T2O) production during the conversion of tritiated A to E1 (see in vitro experimental section), we need the formation rate of T2O, which is simply:(15)

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), E1 (pmoles), E2 (pmoles), and FSH (pg), simple diffusion kinetics were assumed. The hormone quantity in a GC (XGC) has a rate of change equal to:(16)where Kin,X (ml/min) is the rate of medium (“med”) uptake by the GC, Kout,X (ml/min) the rate of excretion by the GC, Xmed (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), Vmed (ml) the volume of culture medium for one GC (total volume divided by the number of cells), and VGC (ml) the volume of one GC. Kin,X was computed by dividing Kout,X by the extra- over intra-cellular partition coefficient Roi,X (unitless), given in the literature [26]:

(17)Conversely, the hormone quantity for one cell in the medium (Xmed) has a rate of change equal to:(18)

The cellular kinetics of A, T, E1, and E2 quantities depend on the entry in and exit from the cell and on their metabolism by Cyp19, Hsd17b1, or Hsd17b2:(19)(20)(21)(22)where the αi are the Michaelis-Menten metabolic terms described in equations 9 to 14. The diffusion and transport of T2O in the in vitro system was not modeled, as the total quantity of T2O 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 (Kin) and exit rate constant (Kout) for a cell, on the hormone concentrations in each cell, and on the number of cells (NGCs or Nothers). The differential equations for the “GC” and the “other cell” compartments are:(23)(24)where Vext and Vothers are the volumes of the extracellular and “other cell” compartments, respectively.

The differential equation for quantity Xext in the extracellular compartment is:(25)where Qinput,X (pmoles or pg/min) is the rate of input of hormone X in the ovary (coming from blood), Fov (ml/min) the efflux of X from the ovary (clearance by blood flow), and Vov,diestrus the ovarian volume at diestrus (which was set at 0.05 ml [27]). For mimicking the female estrus cycle in vivo, Qinput,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). Qinput,X is determined as:

thumbnail
Figure 3. Experimental data vs predictions for FSH and sex steroid hormones in normal cycling rat.

The black line represents mean model predictions with 95% confidence interval (grey band); points represent our experimental observations (mean of 10 measurements ± standard deviation).

https://doi.org/10.1371/journal.pone.0053891.g003

(26)where Qbase,X (pmoles or pg/min) is the constant baseline concentration of hormone X, Qscale,X (unitless) the constant scale for hormone X magnitude, and Qshape,X (pmoles or pg/min) the variable magnitude of hormone X (adjusted to match the known hormone concentrations).

The time courses of NGCs, Vext, and Vothers 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.

thumbnail
Table 3. Prior distributions of the model parameters (for one granulosa cell) to be calibrated by MCMC sampling.

https://doi.org/10.1371/journal.pone.0053891.t003

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 νmRNA.ε. That random effect was modeled with a variability factor σL (see equations 1, 2, and 4), assumed to be log-normally distributed around a mean μσ, with variance Σ1. The hyperparameters μσ and Σ1 were in turn assigned vague prior distributions (Table 3). The individual random effects σL (one per data set used, see Information S1), μσ, and Σ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 (Σ2) and hormone measurements (Σ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 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 (α1 to α6, equations 9 to 14) as a function of time, to determine the predominant reactions for the conversion of A to E2.

thumbnail
Table 4. Summary statistics of the parameter posterior distributions after Bayesian calibration of the in vitro model.

https://doi.org/10.1371/journal.pone.0053891.t004

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”.

thumbnail
Table 5. Model parameter distributions used to describe in vivo variability (in addition to those of Table 4).

https://doi.org/10.1371/journal.pone.0053891.t005

Predictive Simulations of Endocrine Disruption

To evaluate the capacity of the above model to predict in vivo effects of EDCs on E2 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 Km fold-changes fX (equations 13, 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 E2 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.

thumbnail
Table 6. Modulation (fold-change) of steroidogenic enzymes mRNA levels and aromatase enzymatic activity following exposure of granulosa cells to selected chemicals.

https://doi.org/10.1371/journal.pone.0053891.t006

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].

Results

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 ( 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. Vmax and Km 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 Vmax prior for A and T, posterior distributions revealed a 3-fold higher Vmax for T. On the contrary, according to posterior distributions, the aromatase Km for A is 5-fold smaller than the one for T. The posterior distributions for Hsd17b1 Vmax and Km were modified by a factor 1 to 2 for ξHsd17b1,E1, ξHsd17b1,A, and λHsd17b1,A, and by a factor 3 for λ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 (Σ2), and 48% for hormone measurements (Σ3).

Flux Analyses of In Vitro and In Vivo Experiments

Figure 4A, B shows the results of A, E1, T, and E2 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 E1 increased from 7.29×10−9 pmoles/min/cell (without FSH) to 8.72×10−8 pmoles/min/cell (with FSH). The other reaction relative values show that the preferential pathway for E2 synthesis in GCs in vitro is conversion of A into E1, which is then converted into E2, both with or without FSH.

thumbnail
Figure 4. Flux analyses of in vitro and in vivo experiments.

Graphs A and B represent the in vitro flux analysis of steroid hormones conversion at 48 h after addition of 200 nM A into the medium, without or with FSH 20 ng/ml. Graphs C, D, and E illustrate the in vivo flux analysis of steroid hormones conversion at several times of the estrus cycle (corresponding to diestrus, proestrus, and estrus stages). The aromatization reaction of A into E1 is taken as the reference reaction for each condition. The flux values for that reference were 7.29×10−9 pmoles/min/cell in vitro without FSH, 8.72×10−8 pmoles/min/cell in vitro with FSH, 6.09×10−9 pmoles/min/cell in vivo in the diestrus stage, 6.17×10−9 pmoles/min/cell in the proestrus stage, and 5.10×10−9 pmoles/min/cell in the estrus stage of the estrous cycle. Values for the other reactions in each condition are relative to the corresponding reference. Arrow thicknesses are proportional to the flux absolute values.

https://doi.org/10.1371/journal.pone.0053891.g004

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 E1 increased from 5.10×10−9 pmoles/min/cell at the estrus stage to 6.09×10−9 pmoles/min/cell at the diestrus stage and to 6.17×10−9 at the proestrus stage. Those in vivo results, which show a preferential pathway for E2 synthesis trough E1 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 E2 in vivo than in vitro, or the relative importance of the retroconversions of T and E2 to A and E1, 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 E2, whereas the values for E1 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 E2 levels following EDC oral exposure. E2 levels were significantly increased in atrazine-treated females, whereas no statistically significant alteration of E2 was observed with bisphenol A, methoxychlor, and vinclozolin treatment. As far as vinclozolin-treated rats are concerned, one of those showed an elevated E2 ovarian concentration.

thumbnail
Figure 5. Experimental data vs predictions of estradiol levels in control and EDC-treated female rats at the diestrus stage.

Experimental data are represented by points (n = 8 for control data, n = 4 for EDC-treated animals data). Statistical distributions of the model predictions are represented by boxplots (showing the distribution quartiles). Control is for atrazine 200 mg/kg, bisphenol A 200 mg/kg, and vinclozolin 100 mg/kg; control 2 is for letrozole 5 mg/kg. ATZ: atrazine; BPA: Bisphenol A; MXC: methoxychlor; VCZ: vinclozolin; LET: letrozole.

https://doi.org/10.1371/journal.pone.0053891.g005

In vivo data extracted from the literature showed a significant decrease of E2 in letrozole-treated rat ovaries, compared to control [24].

Predictions for Ovarian Estradiol Concentrations in EDC-treated 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 E2 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 E2 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 E2 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-pituitary-gonadal (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 E2; Watanabe et al.’s model [39] simulates synthesis and feedback loops for T, 11-ketotestosterone, E2, and vitellogenin plasma concentrations in male fathead minnow; Li et al.’s model [40] simulate E2, 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 cross-validate 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 E1 and E2 production by the ovary in vivo. The results showed that the model was able to accurately simulate ovarian E2 concentrations during normal cycling in female rats. The results for E1 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 E1 by aromatase (in spite of its conversion into T by Hsd17b1) seems due to differences in Km values of androstedione for aromatase and Hsd17b1, rather than to differences in Vmax values.

The flux analyses indicate that the preferential pathway for E2 synthesis involves E1 both in vitro and in vivo. They also point out the need to perform toxicity testing experiments under FSH-controlled conditions.

Flux analyses show clear differences between in vitro and in vivo conditions. For example, steroid inactivation reaction fluxes (T to A and E2 to E1) are ranged from 10−3 to 10−6 pmoles/min/cell in vitro, and ranged from 10−4 to 10−12 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 Km or mRNA levels on the basis of transcriptomic and enzymatic activity data obtained in vitro for GCs. We limited our predictions to six hours post-exposure, 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 E2 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 atrazine-treated 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 E2 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 hypothalamic-pituitary-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 E2 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 anti-androgenic metabolite M2 [4], and subject to variable amplification 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 E2 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 EDC-treated 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.

https://doi.org/10.1371/journal.pone.0053891.s001

(DOC)

Acknowledgments

The authors would like to thank Karen Watanabe of the Oregon Health & Science University for her helpful comments and discussion.

Author Contributions

Conceived and designed the experiments: NQ FYB. Performed the experiments: NQ FYB. Analyzed the data: NQ FYB. Contributed reagents/materials/analysis tools: NQ FYB. Wrote the paper: NQ FYB.

References

  1. 1. Maranghi F, Mantovani A (2012) Targeted toxicological testing to investigate the role of endocrine disrupters in puberty disorders. Reprod Toxicol 33: 290–296.
  2. 2. Holloway AC, Anger DA, Crankshaw DJ, Wu M, Foster WG (2008) Atrazine-induced changes in aromatase activity in estrogen sensitive target tissues. J Appl Toxicol 28: 260–270.
  3. 3. Fan W, Yanase T, Morinaga H, Gondo S, Okabe T, et al. (2007) Herbicide atrazine activates SF-1 by direct affinity and concomitant co-activators recruitments to induce aromatase expression via promoter II. Biochem Biophys Res Commun 355: 1012–1018.
  4. 4. Molina-Molina JM, Hillenweck A, Jouanin I, Zalko D, Cravedi JP, et al. (2006) Steroid receptor profiling of vinclozolin and its primary metabolites. Toxicol Appl Pharmacol 216: 44–54.
  5. 5. Kelce WR, Monosson E, Gamcsik MP, Laws SC, Gray LE (1994) Environmental hormone disruptors: evidence that vinclozolin developmental toxicity is mediated by antiandrogenic metabolites. Toxicol Appl Pharmacol 126: 276–285.
  6. 6. Cummings AM (1997) Methoxychlor as a model for environmental estrogens. Crit Rev Toxicol 27: 367–379.
  7. 7. Gaido KW, Maness SC, McDonnell DP, Dehal SS, Kupfer D, et al. (2000) Interaction of methoxychlor and related compounds with estrogen receptor alpha and beta, and androgen receptor: structure-activity studies. Mol Pharmacol 58: 852–858.
  8. 8. Hiroi H, Tsutsumi O, Momoeda M, Takai Y, Osuga Y, et al. (1999) Differential interactions of bisphenol A and 17beta-estradiol with estrogen receptor alpha (ERalpha) and ERbeta. Endocr J 46: 773–778.
  9. 9. Petkov PI, Temelkov S, Villeneuve DL, Ankley GT, Mekenyan OG (2009) Mechanism-based categorization of aromatase inhibitors: a potential discovery and screening tool. SAR QSAR Environ Res 20: 657–678.
  10. 10. Odum J, Ashby J (2002) Detection of aromatase inhibitors in vitro using rat ovary microsomes. Toxicol Lett 129: 119–122.
  11. 11. Sandhoff TW, McLean MP (1996) Hormonal regulation of steroidogenic acute regulatory (StAR) protein messenger ribonucleic acid expression in the rat ovary. Endocrine 4: 259–267.
  12. 12. Lovekamp TN, Davis BJ (2001) Mono-(2-ethylhexyl) phthalate suppresses aromatase transcript levels and estradiol production in cultured rat granulosa cells. Toxicol Appl Pharmacol 172: 217–224.
  13. 13. EPA (2007) Integrated Summary Report on Aromatase. Environmental Protection Agency. Available: http://www.epa.gov/endo/pubs/aromatase_isr.pdf. Accessed 2012 August 27.
  14. 14. Stocco C (2008) Aromatase expression in the ovary: hormonal and molecular regulation. Steroids 73: 473–487.
  15. 15. Akinola LA, Poutanen M, Vihko R (1996) Cloning of rat 17 beta-hydroxysteroid dehydrogenase type 2 and characterization of tissue distribution and catalytic activity of rat type 1 and type 2 enzymes. Endocrinology 137: 1572–1579.
  16. 16. Stokes WS (2004) Selecting appropriate animal models and experimental designs for endocrine disruptor research and testing studies. Ilar J 45: 387–393.
  17. 17. Charles GD (2004) In vitro models in endocrine disruptor screening. Ilar J 45: 494–501.
  18. 18. Holme JA, Dybing E (2002) The use of in vitro methods for hazard characterisation of chemicals. Toxicol Lett 127: 135–141.
  19. 19. Wang C, Hsueh AJ, Erickson GF (1981) LH stimulation of estrogen secretion by cultured rat granulosa cells. Mol Cell Endocrinol 24: 17–28.
  20. 20. Quignot N, Desmots S, Barouki R, Lemazurier E (2012a) A comparison of two human cell lines and two rat gonadal cell primary cultures as in vitro screening tools for aromatase modulation. Toxicol In Vitro 26: 107–118.
  21. 21. Lephart ED, Simpson ER, McPhaul MJ (1992) Ovarian aromatase cytochrome P-450 mRNA levels correlate with enzyme activity and serum estradiol levels in anestrous, pregnant and lactating rats. Mol Cell Endocrinol 85: 205–214.
  22. 22. Goldman JM, Murr AS, Cooper RL (2007) The rodent estrous cycle: characterization of vaginal cytology and its utility in toxicological studies. Birth Defects Res B Dev Reprod Toxicol 80: 84–97.
  23. 23. Quignot N, Tournier M, Pouech C, Cren-Olive C, Barouki R, et al. (2012b) Quantification of steroids and endocrine disrupting chemicals in rat ovaries by LC-MS/MS for reproductive toxicology assessment. Anal Bioanal Chem 403: 1629–1640.
  24. 24. Sinha S, Kaseta J, Santner SJ, Demers LM, Bremmer WJ, et al. (1998) Effect of CGS 20267 on ovarian aromatase and gonadotropin levels in the rat. Breast Cancer Res Treat 48: 45–51.
  25. 25. Zheng X, Price CA, Tremblay Y, Lussier JG, Carriere PD (2008) Role of transforming growth factor-beta1 in gene expression and activity of estradiol and progesterone-generating enzymes in FSH-stimulated bovine granulosa cells. Reproduction 136: 447–457.
  26. 26. Breen MS, Breen M, Terasaki N, Yamazaki M, Conolly RB (2009) Computational model of steroidogenesis in human H295R cells to predict biochemical response to endocrine-active chemicals: model development for metyrapone. Environ Health Perspect 118: 265–272.
  27. 27. Quignot N, Arnaud M, Robidel F, Lecomte A, Tournier M, et al. (2012c) Characterization of endocrine-disrupting chemicals based on hormonal balance disruption in male and female adult rats. Reprod Toxicol 33: 339–352.
  28. 28. Bois FY (2009) GNU MCSim: Bayesian statistical inference for SBML-coded systems biology models. Bioinformatics 25: 1453–1454.
  29. 29. Gelman A, Rubin DB (1992) Inference from Iterative Simulation Using Multiple Sequences (with Discussion). Statistical Sciences 7: 457–511.
  30. 30. Liebermeister W, Klipp E (2005) Biochemical networks with uncertain parameters. Syst Biol (Stevenage) 152: 97–107.
  31. 31. Massey FJ (1951) The Kolmogorov-Smirnov test for goodness of fit. J Amer Statistical Assoc 46: 68–78.
  32. 32. Funahashi A, Matsuoka Y, Jouraku A, Morohashi M, Kikuchi N, et al. (2008) CellDesigner 3.5: A Versatile Modeling Tool for Biochemical Networks. Proceedings of the IEEE 96: 1254–1265.
  33. 33. R Development Core Team (2011) R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing, http://www.R-project.org/.
  34. 34. Gordon A, Garrido-Gracia JC, Aguilar R, Guil-Luna S, Millan Y, et al. (2009) Ovarian stimulation with FSH reduces phosphorylation of gonadotrope progesterone receptor and LH secretion in the rat. Reproduction 137: 151–159.
  35. 35. Andersen ME, Thomas RS, Gaido KW, Conolly RB (2005) Dose-response modeling in reproductive toxicology in the systems biology era. Reprod Toxicol 19: 327–337.
  36. 36. Villeneuve DL, Larkin P, Knoebl I, Miracle AL, Kahl MD, et al. (2007) A graphical systems model to facilitate hypothesis-driven ecotoxicogenomics research on the teleost brain-pituitary-gonadal axis. Environ Sci Technol 41: 321–330.
  37. 37. Villeneuve DL, Garcia-Reyero N, Martinovic-Weigelt D, Li Z, Watanabe KH, et al. (2012) A graphical systems model and tissue-specific functional gene sets to aid transcriptomic analysis of chemical impacts on the female teleost reproductive axis. Mutat Res 746: 151–162.
  38. 38. Breen MS, Villeneuve DL, Breen M, Ankley GT, Conolly RB (2007) Mechanistic computational model of ovarian steroidogenesis to predict biochemical responses to endocrine active compounds. Ann Biomed Eng 35: 970–981.
  39. 39. Watanabe KH, Li Z, Kroll KJ, Villeneuve DL, Garcia-Reyero N, et al. (2009) A computational model of the hypothalamic-pituitary-gonadal axis in male fathead minnows exposed to 17alpha-ethinylestradiol and 17beta-estradiol. Toxicol Sci 109: 180–192.
  40. 40. Li Z, Kroll KJ, Jensen KM, Villeneuve DL, Ankley GT, et al. (2011) A computational model of the hypothalamic: pituitary: gonadal axis in female fathead minnows (Pimephales promelas) exposed to 17alpha-ethynylestradiol and 17beta-trenbolone. BMC Syst Biol 5: 63.
  41. 41. Soffker M, Tyler CR (2012) Endocrine disrupting chemicals and sexual behaviors in fish - a critical review on effects and possible consequences. Crit Rev Toxicol 42: 653–668.
  42. 42. Fang H, Tong W, Shi LM, Blair R, Perkins R, et al. (2001) Structure-activity relationships for a large diverse set of natural, synthetic, and environmental estrogens. Chem Res Toxicol 14: 280–294.
  43. 43. Andersen ME, Clewell HJ, 3rd, Gearhart J, Allen BC, Barton HA (1997) Pharmacodynamic model of the rat estrus cycle in relation to endocrine disruptors. J Toxicol Environ Health 52: 189–209.
  44. 44. Bertram R, Li YX (2008) A mathematical model for the actions of activin, inhibin, and follistatin on pituitary gonadotrophs. Bull Math Biol 70: 2211–2228.
  45. 45. Harada N, Honda SI, Hatano O (1999) Aromatase inhibitors and enzyme stability. Endocr Relat Cancer 6: 211–218.
  46. 46. Auvray P, Nativelle C, Bureau R, Dallemagne P, Seralini GE, et al. (2002) Study of substrate specificity of human aromatase by site directed mutagenesis. Eur J Biochem 269: 1393–1405.
  47. 47. Hargrove JL, Hulsey MG, Summers AO (1993a) From genotype to phenotype: computer-based modeling of gene expression with STELLA II. Biotechniques 15: 1096–1101.
  48. 48. Hargrove JL (1993b) Microcomputer-assisted kinetic modeling of mammalian gene expression. Faseb J 7: 1163–1170.
  49. 49. Renwick AG, Soon CY, Chambers SM, Brown CR (1981) Estradiol-17 beta dehydrogenase from chicken liver. J Biol Chem 256: 1881–1887.
  50. 50. Plowchalk DR, Teeguarden J (2002) Development of a physiologically based pharmacokinetic model for estradiol in rats and humans: a biologically motivated quantitative framework for evaluating responses to estradiol and other endocrine-active compounds. Toxicol Sci 69: 60–78.
  51. 51. Odum J, Tinwell H, Van Miller J, Joiner R, Ashby J (2001) The uterotrophic activity of nonylphenol in the rat is not mediated by aromatase enzyme induction. Toxicol Lett 118: 165–169.
  52. 52. Krekels MD, Wouters W, De Coster R (1990) Aromatase inhibition by R 76 713: a kinetic analysis in rat ovarian homogenates. Steroids 55: 69–73.
  53. 53. Ishikura S, Matsumoto K, Sanai M, Horie K, Matsunaga T, et al. (2006) Molecular cloning of a novel type of rat cytoplasmic 17beta-hydroxysteroid dehydrogenase distinct from the type 5 isozyme. J Biochem 139: 1053–1063.
  54. 54. Steckelbroeck S, Watzka M, Reissinger A, Wegener-Toper P, Bidlingmaier F, et al. (2003) Characterisation of estrogenic 17beta-hydroxysteroid dehydrogenase (17beta-HSD) activity in the human brain. J Steroid Biochem Mol Biol 86: 79–92.