# PBK-TD modelling of the gonadotropic axis: Case study with two azole fungicides in female zebrafish

Tu-Ky Ly, Edith Chadili, Olivier Palluel, Karyn Le Menach, Hélène Budzinski, Cleo Tebby, Nathalie Hinfray, Rémy Beaudouin

► To cite this version:

Tu-Ky Ly, Edith Chadili, Olivier Palluel, Karyn Le Menach, Hélène Budzinski, et al., PBK-TD modelling of the gonadotropic axis: Case study with two azole fungicides in female zebrafish. Aquatic Toxicology, 2025, 283, pp.107337. 10.1016/j.aquatox.2025.107337. ineris-05035539

HAL Id: ineris-05035539

https://ineris.hal.science/ineris-05035539v1

Submitted on 15 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# PBK-TD modelling of the gonadotropic axis: Case study with two azole fungicides in female zebrafish

Tu-Ky Ly<sup>a,b,c</sup>, Edith Chadili<sup>c</sup>, Olivier Palluel<sup>c</sup>, Karyn Le Menach<sup>d</sup>, Hélène Budzinski<sup>d</sup>, Cleo Tebby<sup>a</sup>, Nathalie Hinfry<sup>c</sup>, Rémy Beaudouin<sup>a,b,*</sup>

## ARTICLE INFO

Keywords:

PBK-TD

TKTD

Gonadotropic axis

HPG axis

Endocrine Disruptors

Zebrafish

qAOP

QST

## ABSTRACT

Endocrine disruptors (EDs) can disrupt the gonadotropic axis, which consists of the Hypothalamus-Pituitary-Gonads (HPG), notably by altering aromatase (cyp19a), a key enzyme regulating the endocrine system and reproductive function in fish. The effects of EDs can be predicted by integrating both toxicokinetic (TK) and toxicodynamic (TD) processes in order to relate adverse outcomes to external exposures. In this study, we developed a physiologically based kinetic-toxicodynamic model to simulate the disruption of the HPG axis (PBK-TD, hereafter named PBK-HPG) in female zebrafish exposed to either of two aromatase inhibitors, imazalil or prochloraz. The model was calibrated using Bayesian methods and supported by novel experimental data, including measurements of vitellogenin, 17β-estradiol, and 11-ketotestosterone levels, along with in vivo monitoring of the cyp19a1a gene in transgenic cyp19a1a-GFP ebrafish. Seamless integration of a PBK model within a TD model of the HPG-axis, provided the link between external exposure and internal levels of imazalil and prochloraz in key organs, allowing for mechanistic predictions of their inhibitory effects on gonadal aromatase. Our PBK-HPG model accurately predicted both baseline homeostasis and the effects of aromatase inhibition, with all endocrine endpoints including reproductive disruption, i.e., decreased egg production, falling within a twofold range of both experimental and literature data. Therefore, our PBK-HPG model could further support the development of a mechanistic qAOP with TK considerations. The model offers significant potential for improving environmental risk assessments of EDs and possibly other stressors across species.

## 1. Introduction

Quantitative Adverse Outcome Pathways (qAOPs) are of growing interest for environmental risk assessment (Jeong et al., 2025; Paini et al., 2022; Spînu et al., 2020). Among them, mechanistic qAOPs (Spînu et al., 2020) aim to model the progression of a mode of action by quantitatively link a series of key events (KEs) leading to adverse effects, based on the underlying physiological mechanisms (Jeong et al., 2025; Spînu et al., 2020). However, the scope of qAOPs is often limited to relating an adverse outcome to an internal exposure level in the organism under study at a given timepoint or assuming steady-state exposure: most existing qAOP models do not consider integrate either toxicokinetic (TK) or toxicodynamic (TD) aspects of the studied substance. Upstream characterization of the complete TK profile –

absorption, distribution, metabolism and excretion – using Physiologically-based kinetic (PBK) models and coupling these with a mechanistic TD model within a quantitative systems toxicology (QST) framework can enhance the predictive capability and regulatory relevance of qAOPs by explicitly connecting internal exposure dynamics with downstream biological effects (El-Masri et al., 2016; Spînu et al., 2020; Villeneuve et al., 2021). In this approach, the inclusion of a PBK-TD model in a qAOP framework, provides a robust approach to predict the dynamics of KEs, from the internal concentrations temporal dynamic into the different tissues to subsequent adverse effects.

One notable example of a mechanistic qAOP is the model developed to study cytochrome P450 aromatase inhibition, a key mode of action for many endocrine disruptors (EDs) (Conolly et al., 2017). This qAOP includes a mechanistic physiologically-based TD model of the

<footnote><sup>a</sup> Experimental Toxicology and Modeling Unit, INERIS, Verneuil en Halatte 65550, France</footnote>

<footnote><sup>b</sup> UMR-I 02 SEBIO, INERIS, Université de Reims Champagne-Ardenne, Université Le Havre Normandie, Normandie Univ, Verneuil en Halatte 65550, France</footnote>

<footnote><sup>c</sup> Ecotoxicology of Substances and Environments Unit, INERIS, Verneuil-en-Halatte 65550, France</footnote>

<footnote><sup>d</sup> UMR CNRS 5805 EPOC, Université de Bordeaux, Talence 33405, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: remy.beaudouin@ineris.fr (R. Beaudouin).</footnote>
gonadotropic axis (Cheng et al., 2016) as the first step of the qAOP. This model represents the hypothalamus, pituitary gland and gonads (HPG), which altogether is central to regulating the endocrine system and the reproductive function (Conolly et al., 2017). Aromatase, a key enzyme in this regulation, is responsible for converting androgens into estrogens in the gonads and brains of many teleost fish species, including zebrafish. It plays an essential role in maintaining the balance of sex steroid hormones (Hinfray et al., 2006a). Inhibition of aromatase disrupts this balance, particularly by reducing 17β-estradiol (E2) levels, which subsequently decreases vitellogenin (VTG) synthesis, a key factor in egg production (Yilmaz et al., 2018). This disruption leads to impaired reproduction, and the well-characterized mode of action serves as the foundation for this qAOP. Thus, the gonadotropic axis model is crucial, as it captures the entire sequence of events, from aromatase inhibition to the resulting decrease of VTG, thereby providing a thorough assessment of endocrine disruption effects (Conolly et al., 2017; Doering et al., 2024).

Several gonadotropic axis models have been developed for various fish species, including Sciaenid fish (Murphy et al., 2005; Schultz et al., 2001) and Salmonid fish (Gillies et al., 2016; Kim et al., 2006; Sundling et al., 2014). For the fathead minnow, the qAOP for aromatase inhibition includes a gonadotropic axis model which describes three primary compartments – brain, liver, and gonads – following a physiologically-based framework (Cheng et al., 2016; Conolly et al., 2017). These fathead minnow gonadotropic axis models have been tested with substances such as ethinylestradiol (Li et al., 2011; Watanabe et al., 2009), 17β-trenbolone (Li et al., 2011), and a selective non-steroidal aromatase inhibitor, fadrozole (Breen et al., 2016, 2013; Cheng et al., 2016). Latest versions of this gonadotropic axis model incorporates a TK model based solely on in vivo plasma concentration of fadrozole (Cheng et al., 2016). A better prediction of the bioaccumulation of the substance studied in key organs could broaden the domain of applicability of the gonadotropic axis model and enhance its predictive accuracy for a wider range of chemicals (Villeneuve et al., 2021).

In this context, the aim of the present study is to develop an integrative modeling approach within a QST framework to study EDs that act via aromatase inhibition towards reproductive impairment. Specifically, we aim to quantitatively link the internal concentrations of the studied substance in gonads to the aromatase inhibition and subsequent KEs, addressing a gap in current gonadotropic models. To this end, we developed a PBK-TD model of the gonadotropic axis, hereafter referred to as PBK-HPG model, for female zebrafish, for two azole fungicides which are aromatase inhibitors, prochloraz (PCZ) and imazalil (IMZ) (Hinfray et al., 2006a). The PBK part of the PBK-HPG model was based on a previously established zebrafish PBK model for the same substances (Ly et al., 2024).

The development of the PBK-HPG model was supported by both literature and original experimental data, obtained using female *cyp19a1a*-eGFP transgenic zebrafish exposed to IMZ and wild-type AB strain zebrafish exposed to either PCZ or IMZ at three nominal concentrations (3, 30, or 300 µg/L). Bioaccumulation of PCZ and IMZ in key organs, kinetics of circulating E2 and 11-ketotestosterone (11-KT), and levels of circulating VTG were measured in wild-type and *cyp19a1a*-eGFP transgenic zebrafish, respectively. Additionally, in *cy19a1a*-eGFP transgenic zebrafish, ovarian aromatase gene induction was monitored through in vivo GFP fluorescence measurements. The HPG axis homeostasis and its disruption was calibrated with a Bayesian approach, with both the data obtained in this study and literature data. The model's predictive accuracy was assessed using independent literature data.

## 2. Materials and methods

### 2.1. Experimental data: PCZ and IMZ exposure

To calibrate the PBK-HPG model, two experiments were conducted

to provide TK-TD data on unexposed, PCZ-exposed, and IMZ-exposed sexually mature female zebrafish (six-month-old). All experiments were conducted in accordance with Directive 2010/63/EU. The study design and experimental procedures were submitted to the French Animal Care and Experiment Committee for ethical review and received approval, garnering the official license n° 11,177. Fish handling, exposure structure, euthanasia, and tissue samples handling were similar to our previous study (Ly et al., 2024). In both experiments, concentrations of either PCZ or IMZ were measured by high-performance liquid chromatography-tandem mass spectrometry (LC/MSMS) in exposure water and lyophilized tissue samples. Detailed protocols for fish exposures and chemical analysis are provided in Supplementary Materials (Sections 1.1-1.3).

### 2.2. Experiment A: IMZ exposure (OECD TG229)

Experiment A aimed to assess the effects of IMZ on the HPG axis following OECD guideline 229 (TG 229) modified as described in De Oliveira et al. (2020). Transgenic *cyp19a1a*-eGFP zebrafish were acclimated for a week, followed by two weeks of reproductive monitoring to establish baseline egg production. Fish were then exposed to IMZ at three concentrations (3, 30, 300 µg/L) or in control conditions (DMSO, 0.02 %) for 21 days, in a flow-through system with six daily water renewal. Reproduction was monitored daily, and aromatase A gene expression was measured weekly via in vivo GFP fluorescence. After exposure, fish were euthanized for VTG measurement by ELISA, and IMZ bioaccumulation was quantified in liver and water samples.

### 2.3. Experiment B: effects of IMZ or PCZ on sex steroid concentrations

Experiment B investigated the kinetics of unbound E2 and 11-KT in plasma, i.e., deconjugated form and not bound to a transport protein, after short-term exposure to IMZ or PCZ, following the protocol described in Ly et al. (2024). Female wild-type AB zebrafish were exposed to 3, 30, or 300 µg/L of either compound or a control (DMSO, 0.02 %) for at most 14 days. Sampling occurred on the first day (T0), and at 4, 10, and 14 days of exposure, and after a 2-day depuration period. One aquarium was used for each timepoint and concentration level – each aquarium contained 10 fish (five males and five females). Unbound E2 and 11-KT were measured in pooled plasma samples by LC/MSMS, with quantification limits of 0.2 ng/mL for E2 and 0.4 ng/mL for 11-KT.

### 2.4. Statistical analysis

The detailed protocol for statistical analysis is provided in Supplementary Materials (Section 1.4). Analyses were performed using R Statistical Software (v4.3.3; R Core Team 2024). Experimental data were analyzed by fitting a linear mixed-effects model. The level of significance (Type I error; α) for all statistical analyses was set at *p* < 0.05; significant and indicated by asterisks.

### 2.5. Computational HPG axis model description

An overview of the model and a description of the major TD processes are presented in this section. Detailed equations and modeling assumptions are provided in Supplementary Materials (Section 4). All parameters of the PBK-HPG model are provided in Table 2.

### 2.6. Model overview

The HPG axis model is a PBK-TD model that simulates steroidogenesis of sex steroid hormones, driven by the secretion of gonadotropins, and that control the vitellogenesis (Fig. 1).

The structure of the PBK-HPG model was adapted from the PBK model described by Ly et al. (2024). Briefly, in the PBK model, IMZ or PCZ are absorbed through the gills or *per os* and distributed via the
<figure><img src="image_3.png" /><figcaption>Fig. 1. Schematic representation of the female zebrafish PBK-HPG axis model.</figcaption></figure>

plasma across 12 well-mixed compartments. Liver metabolization of IMZ and PCZ was identified as the primary route of elimination.

The TD part of the model uses the same structure as the PBK part: levels of luteinizing hormone (LH) and three steroid hormones, testosterone (T), E2, and 11-KT, are modelled in all 12 compartments with perfusion limited kinetics based on the blood flow (Eq. S2). VTG was only distributed in the venous and arterial blood, liver, gonads, and eggs. The partition coefficients (PC) for T, E2, and 11-KT were estimated using a QSAR model, based on their LogKow (n-Octanol/water partition coefficient) (Table 1) and refine for E2 during the calibration process (Table 2). PC values for LH were set to 1 by default, and PC values for VTG were calibrated (Table 2).

E2 and T kinetics in blood depend on binding to Sex Steroids-Binding Protein (SBP). Levels of the various forms of E2, T, ER, LH, LR, SBP were modelled, namely total E2 (unbound E2, E2 bound to estrogen receptor - ER - namely activated ER, E2 bound to SBP), total T (unbound T, T bound to SBP), total ER (unbound ER, activated ER), total SBP (unbound SBP, bound SBP), total LH (unbound LH, LH bound to LH receptor - namely activated LR), and total LR (unbound LR, activated LR). These variables follow the mass action law within a two-equation system as

<table><caption>**Table 1**
Organ:blood partition coefficients of T and 11-KT.</caption><thead><tr><th></th><th>T</th><th>11-KT</th></tr></thead><tbody><tr><td>Liver</td><td>1.52</td><td>1.83</td></tr><tr><td>Gonads</td><td>1.95</td><td>2.31</td></tr><tr><td>Brain</td><td>1.95</td><td>2.33</td></tr><tr><td>Richly perfused</td><td>1.74</td><td>2.05</td></tr><tr><td>Fat</td><td>1.53</td><td>1.53</td></tr><tr><td>Gastrointestinal tractus</td><td>1.13</td><td>1.37</td></tr><tr><td>Kidney</td><td>4.25</td><td>4.88</td></tr><tr><td>Skin</td><td>1.19</td><td>1.47</td></tr><tr><td>Poorly perfused</td><td>1.68</td><td>2.00</td></tr></tbody></table>

Organ:blood partition coefficients of T and 11-KT were calculated with a QSAR model based on their LogKow. LogKow, T = 3.32; LogKow, 11-KT = 2.

proposed by Murphy et al. (2005). Briefly, the quantity of the target, i.e., ER or SBP, and the quantity of hormone-target complex, i.e., activated ER and LR and bound SBP, depend on first-order dissociation reactions and second order association reactions (Eqs. (1 and 2)).
<table><caption>**Table 2**
List of the PBK-HPG model parameters.</caption><thead><tr><th rowspan="2">Parameters</th><th rowspan="2">Definition</th><th rowspan="2">Unit</th><th rowspan="2">Prior distribution</th><th rowspan="2">References</th><th colspan="2">Posterior distribution</th></tr><tr><th>MPV</th><th>95 % CI</th></tr></thead><tbody><tr><td colspan="7"><strong>E2 concentration</strong></td></tr><tr><td>ke_E2</td><td>Elimination constant of E2 in the gonads and liver</td><td>1/d</td><td>TN ~ (2.0,0.2) [1E-1;20]</td><td>-</td><td>21.5</td><td>[21.0;22.6]</td></tr><tr><td>Km_19a</td><td>Michaelis constant of aromatase A</td><td>pmol/mL</td><td>TN ~ (23.0,0.5) [1E-6;1E6]</td><td>Hinfray (2006)</td><td>0.52</td><td>[0.47;8.95]</td></tr><tr><td>Vmax_19a</td><td>Maximum rate of aromatase A</td><td>pmol/d/mL</td><td>TN ~ (61.0,0.5) [1E-6;1E6]</td><td>Hinfray (2006)</td><td>21.5</td><td>[4.8;23.0]</td></tr><tr><td>k1_E2SBP_ven</td><td>Association constant of E2 to SBP</td><td>mL/pmol/d</td><td>TN ~ (136.0,0.5) [1E-6;1E6]</td><td>Murphy et al. (2005)</td><td>118.6</td><td>[23.09;271.5]</td></tr><tr><td>kd_E2SBP_ven</td><td>Dissociation constant of E2 to SBP calculated with Eq. S7</td><td>1/d</td><td>Calculated</td><td>-</td><td>260.9</td><td>-</td></tr><tr><td>k1_E2ER_liver</td><td>Association constant of E2 to ER in the liver</td><td>ml/pmol/d</td><td>TN ~ (8E-3,0.5) [1E-6;1E6]</td><td>-</td><td>1.87</td><td>[0.35;3.97]</td></tr><tr><td>kd_E2ER_liver</td><td>Dissociation constant of E2 to ER in the liver</td><td>1/d</td><td>TN ~ (10.2,0.5) [1E-6;1E6]</td><td>Smith and Thomas (1990)</td><td>43.8</td><td>[39.0;50.3]</td></tr><tr><td>k1_E2ER_gonads</td><td>Association constant of E2 to ER in the gonads</td><td>ml/pmol/d</td><td>TN ~ (1.0E-3,0.5) [1E-6;1E6]</td><td>-</td><td>4.35</td><td>[0.9;10.0]</td></tr><tr><td>kd_E2ER_gonads</td><td>Dissociation constant of E2 to ER in the gonads</td><td>1/d</td><td>TN ~ (3.3E-2,0.5) [1E-6;1E6]</td><td>-</td><td>60.9</td><td>[46.5;72.0]</td></tr><tr><td>k1_E2ER_brain</td><td>Association constant of E2 to ER in the brain</td><td>ml/pmol/d</td><td>TN ~ (3.0E-4,0.5) [1E-6;1E6]</td><td>-</td><td>2.72</td><td>[0.39;3.99]</td></tr><tr><td>kd_E2ER_brain</td><td>Dissociation constant of E2 to ER in the brain</td><td>1/d</td><td>TN ~ (9.6E-4,0.5) [1E-6;1E6]</td><td>-</td><td>42.4</td><td>[30.4;46.8]</td></tr><tr><td>PC_fat_E2</td><td>Fat:blood partition coefficients of E2</td><td>No unit</td><td>TN ~ (1, 0.5) [1E-6;10]</td><td>-</td><td>1.53</td><td>[0.38;39.0]</td></tr><tr><td>PC_gonads_E2</td><td>Gonads:blood partition coefficients of E2</td><td>No unit</td><td>TN ~ (1, 0.5) [1E-6;10]</td><td>-</td><td>0.11</td><td>[0.007;1.02]</td></tr><tr><td>PC_kidney_E2</td><td>Kidney:blood partition coefficients of E2</td><td>No unit</td><td>TN ~ (1, 0.5) [1E-6;10]</td><td>-</td><td>3.30</td><td>[0.34;7.86]</td></tr><tr><td>PC_pp_E2</td><td>Poorly perfused:blood partition coefficients of E2</td><td>No unit</td><td>TN ~ (1, 0.5) [1E-6;10]</td><td>-</td><td>0.012</td><td>[0.005;0.92]</td></tr><tr><td>PC_liver_E2</td><td>Liver:blood partition coefficients of E2</td><td>No unit</td><td>TN ~ (1, 0.5) [1E-6;10]</td><td>-</td><td>1.49</td><td>[0.16;2.66]</td></tr><tr><td>PC_skin_E2</td><td>Skin:blood partition coefficients of E2</td><td>No unit</td><td>TN ~ (1, 0.5) [1E-6;10]</td><td>-</td><td>0.59</td><td>[0.07;1.86]</td></tr><tr><td>PC_GIT_E2</td><td>Gastrointestinal tractus:blood partition coefficients of E2</td><td>No unit</td><td>TN ~ (1, 0.5) [1E-6;10]</td><td>-</td><td>1.30</td><td>[0.07;1.79]</td></tr><tr><td>PC RP_E2</td><td>Richly perfused:blood partition coefficients of E2</td><td>No unit</td><td>TN ~ (1, 0.5) [1E-6;10]</td><td>-</td><td>1.60</td><td>[0.19;3.03]</td></tr><tr><td>PC_brain_E2</td><td>Brain:blood partition coefficients of E2</td><td>No unit</td><td>TN ~ (1, 0.5) [1E-6;10]</td><td>-</td><td>1.32</td><td>[0.26;3.54]</td></tr><tr><td colspan="7"><strong>11-KT concentration</strong></td></tr><tr><td>Km_11b</td><td>Michaelis constant of 11β-HSD</td><td>pmol/mL</td><td>TN ~ (53.7,0.5) [1E-6;1E6]</td><td>Watanabe et al. (2009)</td><td>50.0</td><td>[0.4;81.1]</td></tr><tr><td>Vmax_11b</td><td>Maximum rate of 11β-HSD</td><td>pmol/d/mL</td><td>TN ~ (240.0,0.5) [1E-6;1E6]</td><td>Watanabe et al. (2009)</td><td>3.4E-1</td><td>[2.9E-2;4.5]</td></tr><tr><td>ke_KT</td><td>Elimination constant of 11-KT in venous blood</td><td>1/d</td><td>TN ~ (2.4,0.5) [1E-6;1E6]</td><td>Watanabe et al. (2009)</td><td>1.32</td><td>[0.7;4.2]</td></tr><tr><td colspan="7"><strong>T concentration</strong></td></tr><tr><td>Km_chol_gonads</td><td>Michaelis constant for T de novo synthesis</td><td>pmol/mL</td><td>TN ~ (4.0,0.5) [1E-6;1E6]</td><td>-</td><td>5.37</td><td>[0.7;6.2]</td></tr><tr><td>Vmax_chol_gonads</td><td>Maximum rate for T de novo synthesis</td><td>pmol/d/mL</td><td>TN ~ (100.0,0.5) [1E-6;1E6]</td><td>-</td><td>37.2</td><td>[5.7;45.3]</td></tr><tr><td>Ratio_chol_LR</td><td>Ratio between cholesterol and LR for T synthesis</td><td>No unit</td><td>Fixed</td><td>-</td><td>1.83</td><td>-</td></tr><tr><td>k1_TSBP_ven</td><td>Association constant of T to SBP</td><td>mL/pmol/d</td><td>TN ~ (136.0,0.5) [1E-6;1E6]</td><td>Murphy et al. (2005)</td><td>78.6</td><td>[21.6;271.4]</td></tr><tr><td>kd_TSBP_ven</td><td>Dissociation constant of T to SBP calculated with Eq. S7</td><td>1/d</td><td>Calculated</td><td>-</td><td>141.5</td><td>-</td></tr><tr><td>ke_T</td><td>Elimination constant of T in the gonads</td><td>1/d</td><td>TN ~ (0.1,0.5) [1E-6;10]</td><td>-</td><td>0.004</td><td>[9.1E-4;0.03]</td></tr><tr><td>nT</td><td>Hill coefficient for T de novo synthesis</td><td>No unit</td><td>Fixed</td><td>-</td><td>1</td><td>-</td></tr><tr><td colspan="7"><strong>LH concentration</strong></td></tr><tr><td>P_LH_brain</td><td>Basal production of LH</td><td>pmol/d</td><td>Fixed</td><td>Cheng et al. (2016)</td><td>0.2</td><td>-</td></tr><tr><td>k1_LHLR_gonads</td><td>Association constant of LH to LR in the gonads</td><td>ml/pmol/d</td><td>Fixed</td><td>Miwa et al., 1994</td><td>1.1</td><td>-</td></tr><tr><td>kd_LHLR_gonads</td><td>Dissociation constant of LH to LR in the gonads</td><td>1/d</td><td>Fixed</td><td>Miwa et al., 1994</td><td>8.3</td><td>-</td></tr><tr><td>ke_LH</td><td>Elimination constant of LH in the blood</td><td>1/d</td><td>Fixed</td><td>Cheng et al. (2016)</td><td>2.4</td><td>-</td></tr><tr><td>Ki_E2ER</td><td>Inhibition constant of LH production by activated ER</td><td>pmol/mL</td><td>Fixed</td><td>Cheng et al. (2016)</td><td>2.6E-2</td><td>-</td></tr><tr><td colspan="7"><strong>VTG concentration</strong></td></tr><tr><td>Km_VTG_liver</td><td>Michaelis constant for VTG synthesis</td><td>pmol/mL</td><td>TN ~ (2.42E-4, 0.3) [1E-6;1E6]</td><td>-</td><td>3.4E-4</td><td>[3.4E-4;3.5E-4]</td></tr><tr><td>Vmax_VTG_liver</td><td>Maximum rate for VTG synthesis</td><td>pmol/d/mL</td><td>TN ~ (25.0, 0.3) [11.0;1E3]</td><td>-</td><td>860.2</td><td>[711.8;985.7]</td></tr><tr><td>nVTG</td><td>Hill index for VTG synthesis</td><td>No unit</td><td>TN ~ (1.0, 0.3) [1E-2;2.0]</td><td>-</td><td>29.7</td><td>[24.0;41.7]</td></tr><tr><td>PC_liver_VTG</td><td>Liver:blood partition coefficient for VTG</td><td>No unit</td><td>TN ~ (1.0E-2, 1.0) [1E-6;10.0]</td><td>-</td><td>7.1E-3</td><td>[5.8E-3;8.4E-3]</td></tr><tr><td>PC_gonads_VTG</td><td>Gonads:blood partition coefficient for VTG</td><td>No unit</td><td>TN ~ (1.0, 1.0) [1E-6;1E2]</td><td>-</td><td>9.8E-2</td><td>[6.8E-2;1.4E-1]</td></tr><tr><td>ke_egg</td><td>Gonads:eggs transfert constant of VTG</td><td>1/d</td><td>TN ~ (0.1, 1.0) [1E-6;10]</td><td>-</td><td>0.2</td><td>[0.1;0.3]</td></tr></tbody></table>

(continued on next page)
<table><caption>**Table 2 (continued)**</caption><thead><tr><th rowspan="2">Parameters</th><th rowspan="2">Definition</th><th rowspan="2">Unit</th><th rowspan="2">Prior distribution</th><th rowspan="2">References</th><th colspan="2">Posterior distribution</th></tr><tr><th>MPV</th><th>95 % CI</th></tr></thead><tbody><tr><td colspan="7"><strong>Prochloraz TD parameters</strong></td></tr><tr><td>Ki_AF</td><td>PCZ inhibition constant for gonadal aromatase</td><td>pmol/mL</td><td>TN ~ (100, 0.5) [10;1000]</td><td>-</td><td>255.6</td><td>[215.0;538.7]</td></tr><tr><td>Smax_aro</td><td>Maximum induction ratio of aromatase</td><td>No unit</td><td>TN ~ (2.0,0.3) [1;10]</td><td>-</td><td>1.68</td><td>[1.09;2.46]</td></tr><tr><td>EC_50_aro</td><td>Half maximal effective concentration for aromatase induction</td><td>No unit</td><td>TN ~ (1.0,0.3) [1;10]</td><td>-</td><td>0.62</td><td>[0.29;0.84]</td></tr><tr><td>K_aro</td><td>Constant for aromatase induction</td><td>No unit</td><td>TN ~ (0.5,0.3) [1E-1;10]</td><td>-</td><td>0.17</td><td>[0.07;0.28]</td></tr><tr><td colspan="7"><strong>Imazalil TD parameters</strong></td></tr><tr><td>Ki_AF</td><td>IMZ inhibition constant for gonadal aromatase</td><td>pmol/mL</td><td>TN ~ (100, 0.5) [10;1000]</td><td>-</td><td>30.1</td><td>[19.4;76.9]</td></tr><tr><td>Smax_aro</td><td>Maximum induction ratio of aromatase</td><td>No unit</td><td>TN ~ (1.8,0.3) [0,5;10]</td><td>-</td><td>1.5</td><td>[1.0;2.2]</td></tr><tr><td>EC_50_aro</td><td>Half maximal effective concentration for aromatase induction</td><td>No unit</td><td>TN ~ (0.5,0.5) [1E-2;5]</td><td>-</td><td>0.55</td><td>[0.28;0.83]</td></tr><tr><td>K_aro</td><td>Constant for aromatase induction</td><td>No unit</td><td>TN ~ (0.2,0.5) [1E-2;2]</td><td>-</td><td>0.05</td><td>[0.02;0.13]</td></tr></tbody></table>

TN stands for the truncated normal law (prior, coefficient variation, lowest bound, highest bound); U stands for the uniform law (lowest bound, highest bound). MPV stands for the Most Probable Value (mode). 95 % CI stands for 95 % Credibility Interval. The rationale behind each assumed prior value is detailed in Supplementary Materials (Section 4.2).

$$ \frac{d(Q_{R,i})}{dt} = v_i \times \left( k_{dRh,i} \times C_{Rh,i} - k_{1Rh,i} \times C_{R,i} \times \frac{C_{h,i}}{PC_{h,i}} \right) \quad (1) $$

$$ \frac{d(Q_{Rh,i})}{dt} = v_i \times \left( k_{1Rh,i} \times C_{R,i} \times \frac{C_{h,i}}{PC_{h,i}} - k_{dRh,i} \times C_{Rh,i} \right) \quad (2) $$

where $Q_{R,i}$ and $Q_{Rh,i}$ is the quantity of target (i.e., receptor or protein) and hormone-target complex (pmol), $v_i$ is the volume of the compartment $i$ (mL), $kd_{Rh,i}$ is the first-order dissociation constant of hormone-target complex in the compartment $i$ (1/d), $k_{1Rh,i}$ is the second-order association constant of the hormone-target complex in the compartment $i$ (mL/pmol/d), $C_{Rh,i}$ is the concentration of the hormone-target complex in the compartment $i$ (pmol/mL), $C_{R,i}$ is the concentration of the target in the compartment $i$ (pmol/mL), $C_{h,i}$ is the concentration of the free hormone in the compartment $i$ (pmol/mL), $PC_{h,i}$ is the partition coefficient of the hormone in the compartment $i$ (no unit).

A first-order general elimination process (Eq. S5) was used for 11-KT, and LH in the venous compartment, accounting for multiple elimination routes. T and E2 elimination were added in the liver.

## 2.7. Physiological processes

This section presents the major processes included within the model. Gonadotropins, which regulate steroid hormone synthesis in the ovary, were represented solely by LH. LH synthesis was modeled with zero-order kinetics in the brain, regulated by negative feedback from activated ER (Eq. (S6)). Steroidogenesis included the biosynthesis of T, E2, and 11-KT, using Hill functions as proposed by Li et al. (2011). We proposed that T biosynthesis was stimulated by activated LR and its ratio to cholesterol concentration in the gonads (Eq. (S9)). 11-KT biosynthesis was induced by 11β-hydroxysteroid dehydrogenase, depending on T concentration in the gonads (Eq. (S9–10)). E2 biosynthesis was under the control of aromatase A, depending on T concentration in the gonad. Vitellogenesis was modeled with a Hill function, dependent on activated ER levels in the liver as proposed by Li et al. (2011) (Eq. (S13)). In addition, we considered that synthesized VTG was passively distributed from the blood into the gonads, and the accumulated VTG was eliminated from the gonads at spawning every two days (Beaudouin et al., 2015) (Eq. (S14)).

## 2.8. Aromatase inhibition

The disruption of the gonadotropic axis was modeled as being a result of the inhibition of the aromatase A enzyme activity. Competitive

binding between T and azole fungicides was modelled with a Michaelis-Menten kinetics for competitive inhibition, to lower the rate of conversion in the presence of azole fungicides, as proposed by Breen et al. (2013). In addition, we added the induction of cyp19a1a gene expression, in response to the decrease in aromatase enzyme activity, by adjusting aromatase Vmax<sub>19A</sub> via the level of aromatase present relative to unexposed steady state. The modulating parameter of Vmax<sub>19A</sub>, representing the relative quantity of aromatase A compared to the unexposed condition (Q<sub>rel,aro</sub>), was driven by variations in unbound E2 concentration in the blood, as captured within a three-equation system (Eqs. (3–5)) and incorporated in Eq. (6) representing E2 synthesis in gonads. A decrease in unbound E2 levels due to competitive inhibition by azole fungicides would result in an increase in aromatase levels.

$$ \Delta_{E2} = C_{E2,unexposed} - C_{E2,exposed} \quad (3) $$

$$ I = 1 + \frac{Smax_{aro} \times \Delta_{E2}}{SC50_{aro} + \Delta_{E2}} \quad (4) $$

$$ \frac{d(Q_{rel,aro})}{dt} = K_{aro} \times (I - Q_{rel,aro}) \quad (5) $$

where $\Delta_{E2}$ is the difference between the venous E2 concentration during unexposed ($C_{E2,unexposed}$) to exposed ($C_{E2,exposed}$) period of time (pmol/mL), $I$ is the induction factor (no unit), $Smax_{aro}$ the maximum effect (no unit), $SC50_{aro}$ the concentration of E2 producing 50 % of the maximum effect (pmol/mL), $Q_{rel,aro}$ is the relative quantity of aromatase A compared to the unexposed condition (no unit), and $K_{aro}$ is the first-order rate of the response in terms of aromatase gene induction (1/d).

$$ \frac{d(Q_{synth,E2})}{dt} = \frac{Q_{rel,aro} \times Vmax_{19A} \times \frac{C_{T,i}}{PC_{T,i}}}{[Km_{19A} \times (1 + \frac{C_{AF,i}}{K_{iAF}})] + \frac{C_{T,i}}{PC_{T,i}}} \quad (6) $$

where $Q_{synth,E2}$ is the quantity synthesized of E2 in the gonads (pmol), $Q_{rel,aro}$ is the relative quantity of aromatase (no unit), $Vmax_{19A}$ is the maximum rate of aromatase A (pmol/d), $Km_{19A}$ is the substrate concentration at which the rate is equal to half Vmax of aromatase A (pmol/mL), $C_{T,i}$ is the concentration of T in the gonads or brain (pmol/mL), and $PC_{T,i}$ is the partition coefficient of T in the gonads or brain (no unit), $C_{AF,i}$ is the concentration of azole fungicides in the gonads or brain (pmol/mL), and $K_{iAF}$ is the inhibition constant of azole fungicide on the aromatase A (pmol/mL).
## 2.9. Model simulations, calibrations and evaluations

Fish growth was assumed to be null during the short 16- or 21-day observation period. As described in Ly et al. (2024), each experiment was meticulously described in the input file, including measured concentrations of PCZ or IMZ in exposure water, water temperature, oxygen saturation, aquarium volume, feeding and initial body weight of the fish. To ensure a physiological steady state in the model before exposure, all simulated exposures began after a burn-in period of 400 days. Calculations were performed using R version 3.6.1. and GNU MCSim v6.2.3 (Bois, 2009). MCSim model codes are provided under the DOI https://doi.org/10.5281/zenodo.13854340 with a Creative Commons Attribution 4.0 license.

A sensitivity analysis (SA) was conducted to guide model calibration as detailed in Supplementary Materials (Section 2), combining two successive analyses: (i) the Morris' method to identify a subset of the most influential parameters among the 51 parameters of the global model, (ii) the Sobol's method on this subset of parameters to quantify their impact on the concentrations of interest.

The model calibration was carried out using Bayesian method in two consecutive steps. To facilitate the handling and interpretation of experimental results in the model, notably during exposure period, E2 and VTG levels were expressed as their relative concentration to their own unexposed steady state. The first step involved calibrating the physiological parameters under baseline conditions, utilizing the data obtained in unexposed conditions throughout experiments. TK-TD parameters related to the HPG axis were calibrated based on measured concentrations of unbound E2 in plasma and gonads, SBP in blood, on levels of activated ER in the brain, liver, and gonads, and on measured VTG levels in blood, spawned eggs, and liver (Table S7, Section 4.3). The second step was the calibration of the TK-TD parameters of the HPG axis representing the effects of IMZ and PCZ exposure using data from experiments A and B with three concentration levels. All IMZ and PCZ TK parameters were fixed (Ly et al., 2024) (Table 2).

Finally, the calibrated PBK-HPG model was evaluated upon its capacity to predict data in zebrafish from the literature. A 2-fold range was considered, as recommended by the WHO for in vivo data for good predictions. The baseline condition was assessed on six publications (Table S9, Section 6.1). Effects of PCZ exposure were assessed on three publications (Table S10, Section 6.2). No empirical data were available for IMZ exposure in zebrafish.

# 3. Results

## 3.1. Experimental data

Due to a technical issue during experiment A (spillage of the entire stock solution into the exposure media) causing an abnormal concentration during the first three days of exposure in the low exposure level conditions (Table S2, Section 1.5), only results from the medium and high exposure concentrations are presented from experiment A. Results are presented in the sequence of KEs following azole exposure and subsequent aromatase inhibition: (i) exposure water concentration, (ii) internal concentrations, (iii) sex steroids hormones concentrations, (iv) VTG concentration, (v) *cyp19a1a* gene induction, and (vi) reproduction. Detailed experimental results are presented in Supplementary Materials (Section 1.5).

IMZ concentrations in water were measured during experiment A and no IMZ was detected in the control conditions. Fish were exposed to 28.8 ± 0.87 µg/L over the 21 days in the medium concentration level and 104.8 ± 1.9 µg/L for the high concentration level (mean ± SD, *n* = 3). The weekly measured data are reported in Table S2. These exposures resulted in bioaccumulation of IMZ in the liver of females of 1204.4 ± 132.7 ng/g for the medium concentration level and 4212.0 ± 799.0 ng/g for the high concentration level (mean ± SEM, *n* = 3), after 21 days of exposure (Table S4). PCZ and IMZ water and internal concentrations

from experiment B were presented in Ly et al. (2024).

The effects of aromatase inhibition on E2 and 11-KT levels were assessed with experiment B. The control conditions from both PCZ and IMZ experiments at all sampling times were used to determine physiological baseline levels of E2 in the plasma of female fish resulting in a mean concentration of 1.11 ± 0.11 ng/mL for unbound E2 (mean ± SEM, *n* = 10) (Table S5). Unbound 11-KT was not detected in any sample (LOQ = 0.4 ng/mL).

Over the 14 days of exposure to PCZ at the low concentration of exposure (2.9 µg/L), there was no trend in unbound E2 concentration, with a variation of 18.2% around a mean concentration of 1.13 ± 0.09 ng/mL (mean ± SEM, *n* = 5) (Table S5). At medium (26.0 µg/L) and high concentration levels (227.2 µg/L) of PCZ, a peak in unbound E2 plasma levels was observed on day 10, reaching 3 times and 2.1 times the initial concentration, respectively (Table S5). Excluding this peak, no trend was observed at the medium concentration level during and after exposure, with a variation of 31.9% around a mean concentration of 1.01 ± 0.16 ng/mL over all sampling times (mean ± SEM, *n* = 4) (Table S5). However, after 14 days of exposure to the high concentration of PCZ, unbound E2 plasma concentration decreased by 56.7%, and after the end of the exposure, the return to the physiological baseline was overshot: unbound E2 levels reached 2.26 times the initial concentration after two days of depuration (Table S5).

No trend was observed at the low concentration of IMZ exposure (2.2 µg/L), with a variation of 27.4% around a mean concentration of 1.05 ± 0.13 ng/mL (mean ± SEM, *n* = 5) (Table S5). At the medium (23.3 µg/L) and high concentration levels (238.9 µg/L) of IMZ, a downward trend in unbound E2 level was noted over the 14 days of exposure, with decreases of 25.0% and 77.2% at day 14, respectively (Table S5). In the high exposure condition, after two days of depuration, the return to the physiological baseline had been overshot: unbound E2 levels reached 2.03 times the initial concentration (Table S5). No overshoot was observed at the medium exposure level (Table S5).

The effects of aromatase inhibition on VTG levels in blood were assessed with a 21-day exposure period to IMZ in experiment A. VTG concentrations were 5.1 × 10⁷ ± 4.3 × 10⁶ ng/mL (mean ± SEM, *n* = 14) in the control group (Fig. S1). The medium exposure level (28.8 µg/L) resulted in a significant 47.0% decrease to 2.7 × 10⁷ ± 3.1 × 10⁶ ng/mL (mean ± SEM, *n* = 15) (*P* < 0.05), while the high exposure level (104.8 µg/L) led to an 84.3% decrease to 8.0 × 10⁶ ± 1.8 × 10⁶ ng/mL (mean ± SEM, *n* = 15) (*P* < 0.001) (Fig. S1).

The induction of *cyp19a1a* gene expression, in response to aromatase activity inhibition by IMZ, was monitored during experiment A. No statistically significant difference was observed between conditions at the beginning of the experiment (T0) and between any timepoints in the control condition (Fig. S2). No significant difference was observed between any timepoint and T0 within each exposure condition (Fig. S2). Still, an induction trend was observed at the medium concentrations (28.8 µg IMZ/L) with level 1.4 times higher after 21 days of exposure (T21) in comparison to T0 (Fig. S2).

Finally, the effects of IMZ on reproduction were assessed with experiment A. During the 9-day pre-exposure period, no statistical differences in egg production were observed between conditions (Figs. S3–4). At the medium concentration level (28.8 µg/L), there was no significant decrease in egg production, with 19.3 ± 5.1 eggs per female per day during exposure vs. 22.9 ± 4.4 eggs per female per day before exposure (mean ± SEM, *n* = 3) (Figs. S3–4). However, at the high concentration level (104.8 µg/L), a significant decrease in cumulative egg production per female was observed (represented by the slope of the linear regression compared before and after beginning of exposure), with 13.1 ± 4.3 eggs per female per day during exposure vs. 36.3 ± 7.8 eggs per female per day before exposure (mean ± SEM, *n* = 3) (Figs. S3–4, *P* < 0.001). This represents a decrease of 15.7% at the medium concentration level and 63.9% at the high concentration level when comparing pre-exposure period to the exposure period (Fig. S4).
## 3.2. Model sensitivity analysis

Detailed results of both Morris and Sobol analysis are shown in Supplementary Materials (Section 2). Briefly, using the Morris screening method for the global SA, the most influential parameters on 16 outputs were selected based on their $\mu^*$ and $\sigma$ values (≥ 0.5) (Fig. S5). This resulted in selecting 34 out of 51 parameters from the global SA for the Sobol method (Table S6).

The elimination constant of E2 in the liver, i.e., *ke*<sub>E2</sub>, is the most influential parameter on each of the following output: unbound E2 plasmatic concentration, the relative concentration of E2, unbound E2 or activated ER concentration in the liver, gonads and brain, unbound SBP and SBP bound to E2 or T in the blood and aromatase induction (Fig. S6). The $Vmax_T$ parameter involved in T biosynthesis has also a large influence (Fig. S6). The relative concentration of E2 was mostly influenced by the maximal induction ratio (*Smax*), the half maximal effective concentration for aromatase induction (*EC*50) and the inhibition constant of PCZ or IMZ (*Ki*<sub>AF</sub>) (Fig. S6). Concerning VTG concentration in blood, eggs' clutch, liver and relative expression of VTG concentration, the SA underlined the importance of the VTG transfer constant from gonads to eggs, i.e., *ke*<sub>egg</sub>, but also gonads and liver PC, i. e., PC<sub>gonads,VTG</sub>, PC<sub>liver,VTG</sub> (Fig. S6).

## 3.3. Modeling the baseline condition

Inclusion of data obtained both in exposed and unexposed fish in the calibration of the PBK-HPG model provided a sound model of the baseline condition. Table 2 shows the estimated biochemical parameters values obtained following model calibration based on endpoints listed in Table S7 (Section 4.3). The calibrated PBK-HPG model was able to capture the HPG axis homeostasis for a representative female zebrafish weighing 0.450 g, by predicting the mean concentration of unbound E2 concentration in the gonads (C<sub>E2, gonads</sub>), activated ER concentration in the liver (C<sub>E2ER, liver</sub>), gonads (C<sub>E2ER, gonads</sub>) and brain (C<sub>E2ER, brain</sub>), SBP concentration in venous blood (C<sub>SBP, venous blood</sub>) and VTG concentration in total blood (C<sub>VTG, total blood</sub>), spawned egg (C<sub>VTG, spawned eggs</sub>) and liver (C<sub>VTG, liver</sub>), the equilibrium ratio between unbound E2 or T and their bound fraction to SBP within a 2-fold range (Ratio<sub>T/TSBP</sub> and Ratio<sub>E2/E2SBP</sub>) (Fig. 2). The mean concentration of unbound E2 in the plasma (C<sub>E2, plasma</sub>) and VTG (C<sub>VTG, total blood</sub>) in the blood were perfectly predicted (Fig. 2). Each kinetic curve for every endpoint is presented in Supplementary Materials (Fig. S7, Section 4.3.2). Overall, all datapoints were predicted within a 2-fold range with 46 datapoints over 50 within a 1.15-fold range (92.0 %) (Fig. 2).

## 3.4. Modeling the disruption of the gonadotropic axis

In our model, the disruption of the gonadotropic axis homeostasis is related to the predicted IMZ or PCZ concentrations in the gonads. Substance-specific parameters are presented in Table S8 (Section 5). The validation of the model predictions for the bioaccumulation of either PCZ or IMZ in female zebrafish were presented in Ly et al. (2024). Liver IMZ concentration measurements from experiment A further validated the model's ability to predict IMZ concentrations in the gonads and brain, with all data points predicted within a 2-fold range and 4 out of 6 points being perfectly predicted (Fig S8). All parameters of the PBK-HPG model are listed in Table 2.

The model accurately captured the kinetic profiles of unbound E2 for all conditions to PCZ and IMZ exposure (Fig. 3). At low concentration (2.9 µg PCZ/L, 2.2 µg IMZ/L), there was no distinct profile leading to model predictions similar to the baseline condition (Fig. 3B). At medium concentrations (26.0 µg PCZ/L, 23.3 µg IMZ/L), the model successfully predicted a downward trend (Fig. 3C). At high concentration (227.2 µg PCZ/L, 238.9 µg IMZ/L), the model predicted a decrease to a bottom value of 0.2 corresponding to the relative concentration of E2, for both PCZ and IMZ over the 14-day exposure period (Fig. 3D). Most measured

points were slightly overestimated during the exposure phase, except at day 10 (Fig. 3). After a two-day depuration period, the model showed a compensation effect, as observed in the data, with an overshoot in E2 concentration before returning to baseline levels around 16 days post-exposure (Fig. 3D).

At medium concentration levels of both PCZ (17.2 µg/L) or IMZ (28.8 µg/L), the model predicted the decrease in VTG blood concentration at the end of the 21-day exposure with a fold 1.1 deviation (Fig. 4B). At high concentrations levels (172.6 µg PCZ/L, 104.8 µg IMZ/L), the model accurately predicted the median VTG blood concentrations for both PCZ and IMZ (Fig. 4C).

The kinetics of aromatase A gene induction were well predicted by the model: the median induction at the end of the 21-day exposure to IMZ was accurately predicted for the medium concentration level (Fig. 5B) and with a 1.1 fold deviation for the high concentration level (Fig. 5C). For PCZ exposure, the induction was predicted with a 1.15 fold deviation at the medium concentration level (Fig. 5E) and a 1.11 fold deviation at the high concentration level (Fig. 5F).

Following PCZ exposure, the model predicted a 55 % decrease in the eggs production at the medium concentration level and a 57 % decrease at the high concentration level. For IMZ exposure, the model over-predicted the decrease in eggs production at the medium concentration (28.8 µg/L), with a predicted decrease of 48 % in comparison to the 24 % observed, and at the high concentration (104.8 µg/L), with a predicted decrease of 82 % in comparison to the 48 % observed (Fig. 6).

## 3.5. Assessment of the predictive capacity of the model

The predictive accuracy of the calibrated model was assessed based on literature data for both the baseline condition and its disruption by PCZ. All model predictions were within a 2-fold range in both control (Fig. S9) and exposed conditions (Fig. S10). Detailed results are provided in Supplementary Materials (Section 6).

# 4. Discussion

To develop an integrative approach to study aromatase enzyme inhibition, a PBK model was combined with a TD model of the gonadotropic axis in the zebrafish. Unlike previous models, our PBK-HPG model precisely described the kinetics of the stressor to predict the dynamics of the KEs and adverse outcomes, allowing an enhanced level of prediction for disrupted endocrine pathways. This consideration resulted in excellent predictive performances of both experimental and literature data: both HPG axis homeostasis and effects of either exposure to IMZ or PCZ were predicted within a 2-fold range (Figs. 2–6).

## 4.1. Effects of IMZ and PCZ exposure

IMZ and PCZ are known to decrease circulating E2 and VTG levels, which are compensated by an induction of *cyp19a1a* gene expression, after short-term exposure in adult female zebrafish (Dang et al., 2015, 2018, 2016), and other fish species, e.g., fathead minnow and Japanese medaka (Ankley et al., 2009; Villeneuve et al., 2023; Zhang et al., 2008).

In our study, the plasmatic concentration of E2 decreased rapidly and reached its lowest value after a 14-day exposure period (a reduction of 57 % and 77 % at the highest dose levels for PCZ and IMZ, respectively, Table S5), consistent with previous zebrafish studies (Dang et al., 2018, 2016). Similar reductions were reported in fathead minnows, with over 80 % reduction after 24 h to 21 days at 165 µg IMZ/L and a 50 % decrease after 21 days at 24 µg IMZ/L (Villeneuve et al., 2023, 2021). Altogether, these results support the E2 kinetics observed in our experiments, accurately modelled following aromatase inhibition, with a rapid decline in E2 levels shortly after exposure, which remained stable over the 14-day period.

Along with the decrease in E2 levels, we observed a clear reduction in blood VTG after 21 days of exposure to the medium IMZ level
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Comparison between measured concentrations of physiological endpoint and model prediction of the PBK-HPG axis model after the calibration process. Experimental datasets used are indicated in the legend as for the different matrices where concentrations were measured. Boxplot were used to represent the distribution of E2 plasmatic concentration, ratio between unbound E2 or T and bound to SBP in the blood, and VTG blood concentration to capture the variability of these endpoints. The central line in boxplot represents the median. The box itself represents the interquartile range, which is the range between the first quartile at the bottom (Q1) and the third quartile at the top (Q3). Whiskers indicate the variability outside the upper and lower quartiles. Value of each endpoint and references are detailed in supplementary data (Section 4.3, Table S6). The solid line corresponds to the identity line, then the dashed line to 2-fold-range and the dotted line to the 5-fold-range.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Model predictions were compared to measured relative concentration of E2 (Exposed over Baseline) in the plasma at two concentrations of PCZ or IMZ from experiment B. A) correspond to the control condition. B) correspond to 2.2 µg PCZ/L, 2.9 µg IMZ/L. C) correspond to 26.0 µg PCZ/L, 23.3 µg IMZ/L. D) correspond to 227.2 µg PCZ/L, 238.9 µg IMZ/L. *n* = 1 per time sampling of 5 pooled plasma sample. Model predictions were represented with the 95 % credibility interval computed from the posterior distributions by the ribbon covering the mean prediction of the output. The solid (for PCZ) or dotted line (for IMZ) represents the Most Probable Value - MPV.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 4. Model predictions were compared to measured relative expression of VTG in the blood at two concentrations of PCZ (De Oliveira et al., 2020) or IMZ from experiment A. A) correspond to the control condition. B) correspond to 17.2 µg PCZ/L, 28.8 µg IMZ/L. C) correspond to 172.6 µg PCZ/L, 104.8 µg IMZ/L. n = 14–15 per condition. Boxplot represents the distribution of the data with the central line representing the median. The box itself represents the interquartile range, which is the range between the first quartile at the bottom (Q1) and the third quartile at the top (Q3). Whiskers indicate the variability outside the upper and lower quartiles. Model predictions were represented with the 95 % credibility interval computed from the posterior distributions by the ribbon covering the mean prediction of the output. The solid line (for PCZ) and dotted line (for IMZ) represents the Most Probable Value – MPV.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Model predictions were compared to measured fluorescence of gonadal GFP in cyp19a1a-eGFP transgenic female zebrafish exposed to two concentrations of PCZ (De Oliveira et al., 2020) or IMZ from experiment A. IMZ predictions are represented in A-C. PCZ predictions are represented in d-F. A, D) correspond to the control condition. B, E) correspond to 28.8 µg IMZ/L, 17.2 µg PCZ/L. C, F) correspond to 104.8 µg IMZ/L, 172.6 µg PCZ/L.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Observation and model predictions of the percentage of eggs per female per day after exposure relative to the pre-exposure level. Each observed point corresponds to one aquarium. The predicted point is represented by the Most Probable Value (MPV) with its 95 % credibility interval.</figcaption></figure>

(Fig. S1). In fathead minnows, no significant VTG change was observed after 24 h at 275 µg IMZ/L, but 60-hour exposure led to 37.5 % and 68 % decreases at 15 µg IMZ/L and 165 µg IMZ/L, respectively, persisting through 21 days (Villeneuve et al., 2023, 2021). This suggests that VTG reduction in zebrafish may also occur early, from 60 h to 21 days at high IMZ levels.

Despite E2 reductions, no significant GFP fluorescence induction indicating aromatase inhibition was observed at any IMZ exposure condition or time (Fig. S2). However, fathead minnow showed significant cyp19a1a mRNA inductions at 60 h (165 µg IMZ/L) and at 10 and 21 days (260 µg IMZ/L) (Villeneuve et al., 2023, 2021). These results might highlight the difference in species sensitivity to aromatase activity inhibition, as seen with fadrozole (Doering et al., 2019a, 2019b). However, even if effective, in vivo fluorescence detection of gonadal cyp19a1a expression is less sensitive in comparison to mRNA quantification by qPCR, as previously shown in De Oliveira et al. (2020). Thus, cyp19a1a mRNA quantification by qPCR could enhance statistical power in our study. Lastly, 104.8 µg IMZ/L exposure led to significant egg production decline, comparable to a 53 % reduction in fathead minnows (Villeneuve et al., 2023).

## 4.2. Modeling of the baseline condition

The development of the baseline condition of the HPG axis was based on physiological assumptions presented in Supplementary Materials (Section 4.2). Notably, special efforts have been made to account for censored data to maximize the information used to calibrate the model. The most commonly used method to handle values below detection limits is to substitute each non-quantified sample with a fraction of the detection limit or to set the non-quantified data to 0, or equal to this limit (Helsel, 2006). Depending on the methodology used, there is subsequent risk of misinterpreting results (Helsel, 2006).

In our study, a statistical model to handle censored data was defined to allow the non-quantified samples to vary between the limit of quantification and zero (Supplementary Materials, Section 3). This error model was applied to 11-KT plasma concentrations that were below detection limits of LC/MSMS in this study and to blood T concentrations retrieved from the literature measured by ELISA. The integration of this error model allowed us to inform the model parameter using of the information that 11-KT levels were below the detection limit. The 11-KT predictions given by our model are similar to levels measured in other studies in zebrafish using ELISA kit, with concentrations ranging from

0.08 ng/mL to 0.2 ng/mL (Christianson-Heiska et al., 2007; Liu et al., 2016, 2022). T concentrations were estimated with a theoretical quantification limit of 1 ng/mL, reducing uncertainty compared to ELISA measurements (Table S9). This approach maximized the available information and reduced uncertainty, contributing to the accurate prediction of hormones dynamics.

The model focused on LH dynamics due to its role in ovarian steroidogenesis, excluding FSH, which is critical for folliculogenesis (Truodeau, 2022). However, in view of the addition of an oogenesis dynamic model, it would have been interesting to describe FSH and LH distinctively, as FSH is essential to model folliculogenesis (Lin and Ge, 2009; Truodeau, 2022). This effort has been previously made in the Coho salmon and Rainbow trout (Gillies et al., 2016; Kim et al., 2006). Moreover, the addition of FSH kinetics would allow a better description of feedback mechanisms in fish as regulation can be either positive or negative by circulating steroids hormones on both FSH and LH in most teleost fish (Truodeau, 2022). In zebrafish, the physiological feedback appears to be positive but in response to aromatase inhibition, we considered a negative regulation on LH (Eq. S6), similar to other models developed, to capture adaptive responses to EDs exposure (Breen et al., 2016, 2013; Cheng et al., 2016; Kim et al., 2006; Li et al., 2011; Murphy et al., 2005). However, unlike the model developed in fathead minnow (Breen et al., 2013), this negative feedback does not appear to play a crucial role to the adaptive response in our model as the SA results rather highlighted T biosynthesis and E2 elimination (Fig. S6).

Furthermore, one major modeling assumption is that the model does not account for interindividual variability or circadian rhythm effects, as it predicts the mean concentrations for a representative median female zebrafish. However, in our study, the measured physiological concentrations of E2 and VTG highlighted a variation over time and at a given time, with a coefficient of variation (CV) of 31 % for E2, and 58 % for VTG, reflecting both biological and experimental variability (Table S5, Fig. S1). Consequently, a part of the gap between measured and predicted values could be explained by this variability.

Inter-individual variability on VTG levels in zebrafish can be extremely large. Whole blood VTG measurement variability has been estimated with an intra-laboratory CV of 35–1 19% accounting for interindividual variability, and an inter-laboratory CV of 54–97% (OCDE, 2006). This variability has also been studied in fathead minnow, with VTG concentrations varying by 50–110 % in females (Coady et al., 2014; Watanabe et al., 2007). Unlike synchronized fish, this interindividual variability cannot be attributed to the reproductive cycle (Gillies et al., 2016; Jensen et al., 2001; Watanabe et al., 2007). Nevertheless, the model was able to predict the mean whole blood VTG concentration of $5.1 \times 10^7$ ng/mL (Fig. S7K) in female zebrafish which falls within the literature range of $1.0 \times 10^6$ ng/mL to $8.1 \times 10^7$ ng/mL (De Oliveira et al., 2020; Hoffmann et al., 2006; OCDE, 2006).

Regarding time variation of the E2 concentrations in zebrafish, a variation of up to 33% was observed between morning and evening levels (Dang et al., 2016). Additionally, the reproductive cycle induces variations in sex steroids hormones with 45–100 % variability in E2 levels in the fathead minnow (Jensen et al., 2001; Watanabe et al., 2007) with similar results observed in the stickleback (Roufidou et al., 2018). This impact of the reproductive cycle could explain the variability found in our study. Nevertheless, the model was able to predict the mean E2 plasmatic concentrations measured in this study, $1.11 \pm 0.11$ ng/mL (Fig. S7A), which falls within the literature range in zebrafish of 0.25 ng/mL to 2.2 ng/mL (Table S9).

Nonetheless, modeling the levels in a median individual based on individual data helps to partially reduce uncertainties arising from quantification methods or variability between laboratories. Since the model accurately predicted mean E2 and VTG levels, which fall within the range of values reported in the literature, it can reliably predict mean concentrations from other laboratories, as highlighted in Fig. S9. This further reaffirms the model's representativeness of a median female zebrafish.
### 4.3. Modeling effects of exposure to IMZ or PCZ

To model the impact of IMZ or PCZ on the HPG axis, only their inhibitory effects on gonadal aromatase were considered. Their potential effects on brain aromatase were not included in our model, as it has been shown that only Aromatase A is responsible for the circulating E2 levels in plasma (Tang et al., 2017). Calibration of inhibition constants showed that IMZ is three times more potent than PCZ, consistent with in vitro data from rainbow trout, where IMZ exhibited threefold greater aromatase inhibition (IC50 = 3.75 pmol/mL) compared to PCZ (IC50 = 11.5 pmol/mL) (Hinfry et al., 2006b).

Unlike previous models using mRNA quantifications (Breen et al., 2016, 2013), our model incorporated aromatase A gene upregulation based on in vivo GFP fluorescence from cyp19a1a-GFP zebrafish, allowing real-time monitoring of gene induction in the same fish (Fig. 2). This innovative non-invasive method, though less sensitive, offers a non-invasive alternative that, when combined with VTG levels, enhances assessment of aromatase inhibition (De Oliveira et al., 2020).

While the model accurately captured the overall decline in E2 and VTG levels (Fig. 3–5), it slightly underpredicted the magnitude of these decreases during exposure. Furthermore, the rapid overshoot in E2 concentrations following depuration was not fully captured (Fig. 3). Future calibrations, incorporating additional data points for metabolism parameters such as Km, Vmax, and elimination constant, could improve the model’s ability to simulate the full kinetic profile of E2 and VTG during both exposure and recovery. Despite this, our E2 kinetic curves align well with previous models for fathead minnows, although our VTG decline is smoother at day 2 (Breen et al., 2016, 2013; Cheng et al., 2016).

The model includes a simplified VTG accumulation submodel that assumes linear accumulation in eggs between spawning events, differing from other HPG models with continuous accumulation (Cheng et al., 2016; Li et al., 2011; Watanabe et al., 2009). Neither approach fully captures the dynamics of vitellogenesis, during which the interval between clutches marks the end of vitellogenesis (Stage III to V, 350 to 650 µm) – a phase lasting longer than two days, with non-linear VTG accumulation in the eggs (Connolly et al., 2014). Although, the amount of VTG excreted can still be used to approximate the number of eggs laid, given that a mature egg contains an average of 30 µg of VTG (Selman et al., 1993). Consequently, the model can predict a decrease in the number of eggs laid due to the inhibition of E2 and VTG synthesis.

For PCZ exposure, the model predicts egg production decreases of 4.4 %, 54.6 %, and 56.6 % for low, medium, and high exposure levels, respectively, exceeding the mean 35 % reported in the literature (Dang et al., 2018; De Oliveira et al., 2020). In response to IMZ exposure, the model overpredicted the decrease in egg production rate by fold 2 at medium concentration and 1.7 at high concentration (Fig. 6). Although the simplified assumption on vitellogenesis enabled to accurately predict the effects on reproduction, incorporating a more mechanistic approach describing oogenesis dynamics, could further enhance model predictions (Ballif et al., 2024).

### 4.4. Applicability domain

In this study, the PBK-HPG model was developed and calibrated for two azole fungicides in zebrafish. However, although species-specific variations observed among teleost fish in regulation of the HPG axis, its general mechanism can be considered similar for most teleost fish (McArdle et al., 2020; Trudeau, 2022). The PBK model can be parameterized for IMZ and PCZ in other fish species as well (Ly et al., 2024). Thus, it would be interesting to evaluate the predictive accuracy of this model for other asynchronous teleost fish – Japanese medaka and fathead minnow. However, an adjustment factor would be necessary to account for differences in sensitivity to aromatase inhibition and for the differences in physiological basal levels between other teleost species and zebrafish (Doering et al., 2019b). In addition, this model could also

be applied to other aromatase-inhibiting azole fungicides. To broaden the model’s applicability domain, future refinement could integrate additional endocrine pathways, such as thyroid hormone regulation under the control of thyroid peroxidase, as its inhibition may also reduce E2 and VTG levels, as described in AOP 271 (currently under construction – “Inhibition of thyroid peroxidase leading to impaired fertility in fish”). Thus, it could enhance the model’s utility in assessing the effects of a broader range of EDs and support the development of qAOP networks for environmental risk assessment.

## 5. Conclusion

The PBK-HPG model successfully predicted HPG axis homeostasis and responses to stressors in zebrafish. This model widens the scope of existing HGP axis models since it takes into account the kinetics of the stressors in response to the exposure scenario and levels into a dynamic model of adverse effects from aromatase inhibition to reproductive impacts. This model provides a mechanistic quantitative relationship between internal concentrations of stressors and an ED mode of action at the target organ where the molecular initiating event occurs. Therefore, this PBK-HPG model could be further developed into a mechanistic qAOP with TK considerations by incorporating a detailed model of oocyte growth dynamics (Ballif et al., 2024) and an individual-based model (IBM) based on Dynamic Energy Budget theory (DEB-IBM) (Beaudouin et al., 2015). As a result, by offering a mechanistic quantitative link between internal concentrations of stressors and their adverse outcomes, out PBK-HPG model provides a valuable tool for environmental risk assessment and regulatory decision-making, particularly within an “Integrated Approaches to Testing and Assessment” (IATA) framework (Paini et al., 2022; Worth and Patlewicz, 2016).

### CRediT authorship contribution statement

**Tu-Ky Ly:** Writing – review & editing, Writing – original draft, Visualization, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Edith Chadili:** Investigation. **Olivier Palluel:** Investigation. **Karyn Le Menach:** Writing – review & editing, Methodology, Investigation. **Hélène Budzinski:** Writing – review & editing, Methodology. **Cleo Tebby:** Writing – review & editing, Validation, Software, Formal analysis. **Nathalie Hinfry:** Writing – review & editing, Validation, Supervision, Methodology, Investigation, Funding acquisition, Conceptualization. **Rémy Beaudouin:** Writing – review & editing, Supervision, Software, Methodology, Funding acquisition, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Supplementary materials

Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.aquatox.2025.107337.

### Data availability

Data will be made available on request.
References<fnref n="1" />

Ankley, G.T., Bencic, D.C., Cavallin, J.E., Jensen, K.M., Kahl, M.D., Makynen, E.A., Martinovic, D., Mueller, N.D., Wehmas, L.C., Villeneuve, D.L., 2009. Dynamic nature of alterations in the endocrine system of fathead minnows exposed to the fungicide prochloraz. *Toxicol. Sci.* 112, 344–353.

Ballif, G., Clément, F., Yvinec, R., 2024. Nonlinear compartmental modeling to monitor ovarian follicle population dynamics on the whole lifespan. *J. Math. Biol.* 9, 9.

Beaudouin, R., Goussen, B., Piccini, B., Augustine, S., Devillers, J., Brion, F., Pery, A.R.R., 2015. An individual-based model of Zebrafish population dynamics accounting for energy dynamics. *PLoS One* 10.

Bois, F.Y., 2009. GNU MCSim: bayesian statistical inference for SBML-coded systems biology models. *Bioinformatics* 25, 1453–1454.

Breen, M., Villeneuve, D.L., Ankley, G.T., Bencic, D., Breen, M.S., Watanabe, K.H., Lloyd, A.L., Conolly, R.B., 2016. Computational model of the fathead minnow hypothalamic-pituitary-gonadal axis: incorporating protein synthesis in improving predictability of responses to endocrine active chemicals. *Comparative biochemistry and physiology. Toxicol. Pharmacol.*: CBP 183-184, 36–45.

Breen, M., Villeneuve, D.L., Ankley, G.T., Bencic, D.C., Breen, M.S., Watanabe, K.H., Lloyd, A.L., Conolly, R.B., 2013. Developing predictive approaches to characterize adaptive responses of the reproductive endocrine axis to aromatase inhibition: II. *Comput. Model.*: *Toxicol. Sci.* 133, 234–247.

Cheng, W.Y., Zhang, Q., Schroeder, A., Villeneuve, D.L., Ankley, G.T., Conolly, R., 2016. Editor's highlight: computational modeling of plasma vitellogenin alterations in response to aromatase inhibition in fathead minnows. *Toxicol. Sci.* 154, 78–89.

Christianson-Heiska, I., Smeds, P., Granholm, N., Bergelin, E., Isomaa, B., 2007. Endocrine modulating actions of a phytosterol mixture and its oxidation products in zebrafish (Danio rerio) *Comp. Biochem. Physiol. C Toxicol. Pharmacol.* 145, 518–527.

Coady, K.K., Lehman, C.M., Currie, R.J., Marino, T.A., 2014. Challenges and approaches to conducting and interpreting the amphibian metamorphosis assay and the fish short-term reproduction assay. *Birth Defects Res. Part B: Dev. Reproduct. Toxicol.* 101, 80–89.

Connolly, M.H., Dutkosky, R.M., Heah, T.P., Sayler, G.S., Henry, T.B., 2014. Temporal dynamics of oocyte growth and vitellogenin gene expression in zebrafish (Danio rerio). *Zebrafish* 11, 107–114.

Conolly, R.B., Ankley, G.T., Cheng, W.Y., Mayo, M.L., Miller, D.H., Perkins, E.J., Villeneuve, D.L., Watanabe, K.H., 2017. Quantitative adverse outcome pathways and their application to predictive toxicology. *Environ. Sci. Technol.* 51, 4661–4672.

Dang, Y., Giesy, J.P., Wang, J.H., Liu, C.S., 2015. Dose-dependent compensation responses of the hypothalamic-pituitary-gonadal-liver axis of zebrafish exposed to the fungicide prochloraz. *Aquat. Toxicol.* 160, 69–75.

Dang, Y., Sun, Q., Letcher, R.J., Liu, C.S., 2018. Unexpected observations: exposure to aromatase inhibitor prochloraz did not alter the vitellogenin content of zebrafish ova but did inhibit the growth of larval offspring. *Environ. Sci. Technol.* Lett. 5, 629–634.

Dang, Y., Wang, J., Giesy, J.P., Liu, C., 2016. Responses of the zebrafish hypothalamic-pituitary-gonadal-axis PCR array to prochloraz are dependent on timing of sampling. *Aquat. Toxicol.* 175, 154–159.

De Oliveira, J., Chadili, E., Piccini, B., Turies, C., Maillot-Marechal, E., Palluel, O., Pardon, P., Budzinski, H., Cousin, X., Brion, F., Hinfrey, N., 2020. Refinement of an OECD test guideline for evaluating the effects of endocrine disrupting chemicals on aromatase gene expression and reproduction using novel transgenic cyp19a1a-eGFP zebrafish. *Aquat. Toxicol.* 220.

Doering, J.A., Dubiel, J., Stock, E., Collins, C.H., Frick, I., Johnson, H.M., Lowrey-Dufour, C.M., Miller, J.G.P., Xia, Z., Tomy, G.T., Wiseman, S., 2014. A quantitative adverse outcome pathway for embryonic activation of the aryl hydrocarbon receptor of fishes by polycyclic aromatic aromatic hydrocarbons leading to decreased fecundity at adulthood. *Environ. Toxicol. Chem.* n/a.

Doering, J.A., Villeneuve, D.L., Fay, K.A., Randolph, E.C., Jensen, K.M., Kahl, M.D., LaLone, C.A., Ankley, G.T., 2019a. Quantitative sensitivity to In vitro inhibition of cytochrome P450 aromatase (CYP19) activity among 18 freshwater fishes. *Toxicol. Sci.* 170, 394–403.

Doering, J.A., Villeneuve, D.L., Poole, S.T., Blackwell, B.R., Jensen, K.M., Kahl, M.D., Kittelson, A.R., Feifarek, D.J., Tilton, C.B., LaLone, C.A., Ankley, G.T., 2019b. Quantitative response-response pathways linking aromatase inhibition to decreased fecundity are conserved across three fishes with asynchronous oocyte development. *Environ. Sci. Technol.* 53, 10470–10478.

El-Masri, H., Kleinstreuer, N., Hines, R.N., Adams, L., Tal, T., Isaacs, K., Wetmore, B.A., Tan, Y.-M., 2016. Integration of life-stage physiologically based pharmacokinetic models with adverse outcome pathways and environmental exposure models to screen for environmental hazards. *Toxicol. Sci.* 152, 230–243.

Gillies, K., Krone, S.M., Nagler, J.J., Schultz, I.R., 2016. A computational model of the Rainbow Trout hypothalamus-pituitary-ovary-liver axis. *PLoS Comput. Biol.* 12, e1004874.

Helsel, D.R., 2006. Fabricating data: how substituting values for nondetects can ruin results, and what can be done about it. *Chemosphere* 65, 2434–2439.

Hinfrey, N., 2006. Etude de l'expression des cytochromes P 450 aromatases comme marqueurs biologiques d'une perturbation endocrinienne chez le poisson. *Thèse de doctorat. Toxicologie de l'environnement*. http://www.theses.fr/2006MET2021S.

Hinfrey, N., Palluel, O., Turies, C., Cousin, C., Porcher, J.M., Brion, F., 2006a. Brain and gonadal aromatase as potential targets of endocrine disrupting chemicals in a model species, the zebrafish (Danio rerio). *Environ. Toxicol.* 21, 332–337.

Hinfrey, N., Porcher, J.M., Brion, F., 2006b. Inhibition of rainbow trout (Oncorhynchus mykiss) P450 aromatase activities in brain and ovarian microsomes by various environmental substances. *Compar. Biochem. Physiol. C-Toxicol. Pharmacol.* 144, 252–262.

Hoffmann, J.L., Torontali, S.P., Thomason, R.G., Lee, D.M., Brill, J.L., Price, B.B., Carr, G. J., Versteeg, D.J., 2006. Hepatic gene expression profiling using genechips in zebrafish exposed to 17α-ethynylestradiol. *Aquat. Toxicol.* 79, 233–246.

Jensen, K.M., Korte, J.J., Kahl, M.D., Pasha, M.S., Ankley, G.T., 2001. Aspects of basic reproductive biology and endocrinology in the fathead minnow (Pimephales promelas). *Comparative biochemistry and physiology. Toxicol. Pharmacol.*: CBP 128, 127–141.

Jeong, J., Gasparyan, M., Choi, J., 2025. Advancing the quantitative understanding of adverse outcome pathways: current status, methodologies, and future directions. *Environ. Toxicol. Chem.*

Kim, J., Hayton, W.L., Schultz, I.R., 2006. Modeling the brain-pituitary-gonad axis in salmon. *Mar. Environ. Res.* 62, S426-S432.

Li, Z.H., Kroll, K.J., Jensen, K.M., Villeneuve, D.L., Ankley, G.T., Brian, J.V., Sepulveda, M.S., Orlando, E.F., Lazorchak, J.M., Kostich, M., Armstrong, B., Denslow, N.D., Watanabe, K.H., 2011. A computational model of the hypothalamic-pituitary-gonadal axis in female fathead minnows (Pimephales promelas) exposed to 17 alpha-ethynylestradiol and 17 beta-trenbolone. *BMC Syst. Biol.* 5.

Lin, S.-W., Ge, W., 2009. Differential regulation of gonadotropins (FSH and LH) and growth hormone (GH) by neuroendocrine, endocrine, and paracrine factors in the zebrafish—An in vitro approach. *Gen. Comp. Endocrinol.* 160, 183–193.

Liu, W., Chen, C., Chen, L., Wang, L., Li, J., Chen, Y., Jin, J., Kawan, A., Zhang, X., 2016. Sex-dependent effects of microcystin-LR on hypothalamic-pituitary-gonad axis and gametogenesis of adult zebrafish. *Sci. Rep.* 6, 22819.

Liu, X.S., Lu, X.X., Hong, J.B., Zhang, J., Lin, J.T., Jiang, M.Z., Liu, Q., Choi, Y.O., Zhang, J.J., 2022. Effects of long-term exposure to TDCPP in zebrafish (Danio rerio) - alternations of hormone balance and gene transcriptions along hypothalamus-pituitary axes. *Animal Models Experim. Med.* 5, 239–247.

Ly, T.-K., De Oliveira, J., Chadili, E., Le Menach, K., Budzinski, H., James, A., Hinfrey, N., Beaudouin, R., 2024. Imazalil and prochloraz toxicokinetics in fish probed by a physiologically based kinetic (PBK) model. *Environ. Sci. Pollution Res.*

McArdle, M.E., Freeman, E.L., Staveley, J.P., Ortego, L.S., Coady, K.K., Weltje, L., Weyers, A., Wheeler, J.R., Bone, A.J., 2020. Critical review of read-across potential in testing for endocrine-related effects in vertebrate ecological receptors. *Environ. Toxicol. Chem.* 39, 739–753.

Miwa, S., Yan, L., Swanson, P., 1994. Localization of two gonadotropin receptors in the Salmon Gonad by in vitro ligand autoradiography. *Biol. Reprod.* 50 (3), 629–642.

Murphy, C.A., Rose, K.A., Thomas, P., 2005. Modeling vitellogenesis in female fish exposed to environmental stressors: predicting the effects of endocrine disturbance due to exposure to a PCB mixture and cadmium. *Reprod. Toxicol.* 19, 395–409.

OCDE, 2006. Report of the validation of the 21-days fish screening assay for the detection of endocrine substances (phase 1B). OECD Series on Testing and Assessment. No. 61.

Paini, A., Campia, I., Cronin, M.T.D., Asturiol, D., Ceriani, L., Exner, T.E., Gao, W., Gomes, C., Kruisselbrink, J., Martens, M., Meek, M.E.B., Pamies, D., Pletz, J., Scholz, S., Schüttler, A., Spünu, N., Villeneuve, D.L., Wittwehr, C., Worth, A., Luijten, M., 2022. Towards a qAOP framework for predictive toxicology - linking data to decisions. *Computat. Toxicol.* 21, 100195.

Roufoud, C., Schmitz, M., Mayer, I., Sebire, M., Katsiadaki, I., Shao, Y.T., Borg, B., 2018. Hormonal changes over the spawning cycle in the female three-spined stickleback, Gasterosteus aculeatus. *Gen. Comp. Endocrinol.* 257, 97–105.

Schultz, I.R., Orner, G., Merdink, J.L., Skillman, A., 2001. Dose-response relationships and pharmacokinetics of vitellogenin in rainbow trout after intravascular administration of 17α-ethynylestradiol. *Aquat. Toxicol.* 51, 305–318.

Selman, K., Wallace, R.A., Sarka, A., Qi, X., 1993. Stages of oocyte development in the zebrafish, Brachydanio rerio. *J. Morphol.* 218, 203–224.

Smith, J.S., Thomas, P., 1990. Binding characteristics of the hepatic estrogen receptor of the spotted seatfish. *Cynoscion nebulosus*. *Gen. Comp. Endocrinol.* 77 (1), 29–42.

Spünu, N., Cronin, M.T.D., Enoch, S.J., Madden, J.C., Worth, A.P., 2020. Quantitative adverse outcome pathway (qAOP) models for toxicity prediction. *Arch. Toxicol.* 94, 1497–1510.

Sundling, K., Craciun, G., Schultz, I., Hook, S., Nagler, J., Cavileer, T., Verducci, J., Liu, Y., Kim, J., Hayton, W., 2014. Modeling the endocrine control of vitellogenin production in female rainbow trout. *Math. Biosci.* Eng.: MBE 11, 621–639.

Tang, H., Chen, Y., Liu, Y., Yin, Y., Li, G., Guo, Y., Liu, X., Lin, H., 2017. New insights into the role of estrogens in male fertility based on findings in aromatase-deficient zebrafish. *Endocrinology* 158, 3042–3054.

Trudeau, V.L., 2022. Neuroendocrine control of reproduction in teleost fish: concepts and controversies. *Annu. Rev. Anim. Biosci.* 10, 107–130.

Villeneuve, D.L., Blackwell, B.R., Blanksma, C.A., Cavallin, J.E., Cheng, W.Y., Conolly, R. B., Conrow, K., Feifarek, D.J., Heinis, L.J., Jensen, K.M., Kahl, M.D., Milsk, R.Y., Poole, S.T., Randolph, E.C., Saari, T.W., Watanabe, K.H., Ankley, G.T., 2023. Case study in 21st-century ecotoxicology: using In vitro aromatase inhibition data to predict reproductive outcomes in fish In vivo. *Environ. Toxicol. Chem.*

Villeneuve, D.L., Blackwell, B.R., Cavallin, J.E., Cheng, W.Y., Feifarek, D.J., Jensen, K.M., Kahl, M.W., Milsk, R.Y., Poole, S.T., Randolph, E.C., Saari, T.W., Ankley, G.T., 2021. Case study in 21st century ecotoxicology: using In vitro aromatase inhibition data to predict reproductive outcomes in fish In vivo. *Environ. Toxicol. Chem.*

<footnote n="1">¹ During the preparation of this work the authors used Scribans in order to grammar correct some sentences. After using this tool, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication</footnote>
predict short-term In vivo responses in adult female fish. Environ. Toxicol. Chem. 40, 1155–1170.

Watanabe, K.H., Jensen, K.M., Orlando, E.F., Ankley, G.T., 2007. What is normal? A characterization of the values and variability in reproductive endpoints of the fathead minnow, Pimephales promelas. Comparative biochemistry and physiology. Toxicol. Pharmacol.: CBP 146, 348–356.

Watanabe, K.H., Li, Z., Kroll, K.J., Villeneuve, D.L., Garcia-Reyero, N., Orlando, E.F., Sepulveda, M.S., Collette, T.W., Ekman, D.R., Ankley, G.T., Denslow, N.D., 2009. A computational model of the hypothalamic-pituitary-gonadal axis in male fathead

minnows exposed to 17 alpha-ethinylestradiol and 17 beta-estradiol. Toxicol. Sci. 109, 180–192.

Worth, A.P., Patlewicz, G., 2016. Integrated approaches to testing and assessment. Adv. Exp. Med. Biol. 856, 317–342.

Yilmaz, O., Patinote, A., Nguyen, T., Bobe, J., 2018. Multiple vitellogenins in zebrafish (Danio rerio): quantitative inventory of genes, transcripts and proteins, and relation to egg quality. Fish. Physiol. Biochem. 44, 1509–1525.

Zhang, X.W., Hecker, M., Jones, P.D., Newsted, J., Au, D., Kong, R., Wu, R.S.S., Giesy, J. P., 2008. Responses of the medaka HPG axis PCR array and reproduction to prochloraz and ketoconazole. Environ. Sci. Technol. 42, 6762–6769.