Viviane David

Sandrine Joachim

Jean-Marc Porcher

Rémy Beaudouin
email: remy.beaudouin@ineris.fr

Modelling BPA effects on three-spined stickleback population dynamics in mesocosms to improve the understanding of population effects

Keywords: bisphenol A, individual-based model, risk assessment, mesocosms, three-spined sticklebacks

Bisphenol A (BPA), a well-known endocrine-disrupting chemical, is ubiquitously present in the aquatic environment. Its impacts at the population level on three-spined sticklebacks have been studied in artificial streams with low-dose BPA concentrations. The causes explaining the observed effects remained unclear. Here, we used an individual-based model coupled to a Dynamic Energy Budget model to (i) assess the potential of modelling to predict impacts at the population level using individual level laboratory ecotoxicological endpoints and (ii) provide insight on the mechanisms of BPA toxicity in these mesocosms. To do that, both direct and indirect effects of BPA on three-spined sticklebacks were incorporated in the model. Indeed, direct BPA effects on fish have been identified based on literature data whereas indirect effects on sticklebacks have been taken into account using sampling data of their prey from the exposed artificial streams. Results of the modelling showed that direct BPA effects on fish (impacts on gonad formation, growth, male reproductive behavior, eggs and larvae survival) mainly explained the three-spined stickleback population structure in the mesocosms, but indirect effects were not negligible. Hence, this study showed the potential of modelling in risk assessment to predict the impacts on fish population viability from behavioural and physiological effects measured on organisms.

Introduction

Bisphenol A (BPA) is of a scientific and public interest due to its potential adverse effects on exposed organisms and its ubiquity in the environment (Flint et al. 2012, Im and[START_REF] Im | Fate of Bisphenol A in Terrestrial and Aquatic Environments[END_REF]. This chemical is mainly used for the synthesis of polycarbonate and epoxy resins [START_REF] Beausoleil | Regulatory identification of BPA as an endocrine disruptor: Context and methodology[END_REF]. Furthermore, BPA is one of the most investigated substances for its endocrine disruptor properties and the regulatory consideration of this substance has been a challenge during the last decades [START_REF] Beausoleil | Regulatory identification of BPA as an endocrine disruptor: Context and methodology[END_REF].

Many studies have shown that BPA can impact different aquatic organisms as effects on growth, reproduction and survival were observed for macrophytes, mollusks, crustaceans, insects, zooplankton and amphibians (Flint et al. 2012, Canesi and[START_REF] Canesi | Environmental Effects of BPA: Focus on Aquatic Species[END_REF]. Furthermore, BPA can also impact different fish species as shown by several studies at the individual level (Flint et al. 2012, Canesi and[START_REF] Canesi | Environmental Effects of BPA: Focus on Aquatic Species[END_REF]. More specifically, BPA induces synthesis of vitellogenin at concentrations ranging from 10 to 1360 µg/L depending on the fish species [START_REF] Ishibashi | In Vivo Testing System for Determining the Estrogenic Activity of Endocrine-Disrupting Chemicals (EDCs) in Goldfish (Carassius auratus)[END_REF][START_REF] Hatef | Modulations in androgen and estrogen mediating genes and testicular response in male goldfish exposed to bisphenol A[END_REF][START_REF] Mihaich | Adult fathead minnow, Pimephales promelas, partial life-cycle reproductive and gonadal histopathology study with bisphenol A[END_REF][START_REF] Villeneuve | Ecotoxicogenomics to Support Ecological Risk Assessment: A Case Study with Bisphenol A in Fish[END_REF][START_REF] Allner | Kinetic determination of vitellogenin induction in the epidermis of and perciform fishes: Evaluation of sensitive enzyme-linked immunosorbent assays[END_REF][START_REF] Negintaji | Short-term induction of vitellogenesis in the immature male yellowfin seabream (Acanthopagrus latus) exposed to bisphenol A and 17β-estradiol[END_REF] and fish gonads could also present abnormal structure starting from 1 µg/L [START_REF] Mandich | In vivo exposure of carp to graded concentrations of bisphenol A[END_REF][START_REF] Molina | Endocrine-active compound evaluation: Qualitative and quantitative histomorphological assessment of zebrafish gonads after bisphenol-A exposure[END_REF], Chen et al. 2017). Effects on eggs and embryos were observed such as hatching delay, embryo mortality, embryo abnormality and disruption of the hatching success starting from 10 µg/L [START_REF] Shioda | Effect of certain chemicals on the reproduction of medaka (Oryziaslatipes)[END_REF][START_REF] Aluru | Bisphenol A in oocytes leads to growth suppression and altered stress performance in juvenile rainbow trout[END_REF][START_REF] Staples | Early life-stage and multigeneration toxicity study with bisphenol A and fathead minnows (Pimephales promelas)[END_REF][START_REF] Bhandari | Transgenerational effects from early developmental exposures to bisphenol A or 17α-ethinylestradiol in medaka, Oryzias latipes[END_REF][START_REF] Fan | Maternal Bisphenol A exposure impaired endochondral ossification in craniofacial cartilage of rare minnow (Gobiocypris rarus) offspring[END_REF]. Other impacts on physiological processes, such as growth, development and general metabolism, were also shown in fish (Flint et al. 2012, Canesi and[START_REF] Canesi | Environmental Effects of BPA: Focus on Aquatic Species[END_REF].

In a long-term experiment, effects of BPA on three-spined stickleback (Gasterosteus aculeatus) population dynamics have been tested in artificial streams (i.e. mesocosms) using environmentally realistic exposure conditions (1,10 and 100 µg/L) (de [START_REF] De Kermoysan | Effects of bisphenol A on different trophic levels in a lotic experimental ecosystem[END_REF].

Mesocosm experiments have the advantage to test a chemical in relevant environmental conditions and investigate the effects at different levels of biological organization [START_REF] Caquet | Mesocosms in ecotoxicology (1): Outdoor aquatic systems[END_REF], EFSA 2013). In de [START_REF] De Kermoysan | Effects of bisphenol A on different trophic levels in a lotic experimental ecosystem[END_REF], both individual and population level effects on three-spined sticklebacks were observed. Indeed, gonad morphology of females and males was impacted respectively at 1 µg/L and 10 µg/L respectively. Concerning population level endpoints, the fish length frequencies distribution shifted to higher lengths compared to the control at the 100 µg/L BPA treatment. However, the causes leading to the observed population effects remained unclear. Hence, to improve environmental risk assessment, it is essential to provide a robust understanding of ecosystem functioning and population dynamics under toxicant stress, which is one of the challenges of ecological modelling [START_REF] Goussen | Integrated presentation of ecological risk from multiple stressors[END_REF][START_REF] Forbes | A Framework for Predicting Impacts on Ecosystem Services From (Sub)Organismal Responses to Chemicals[END_REF].

Individual-based models (IBMs) allow to extrapolate ecotoxicological effects from individuals to populations [START_REF] Forbes | The extrapolation problem and how population modeling can help[END_REF]. Furthermore, coupling an IBM to a Dynamic Energy Budget (DEB) model offers great flexibility for integrating individual impacts on bioenergetic processes such as growth and reproduction [START_REF] Beaudouin | Individual-based model of Chironomus riparius population dynamics over several generations to explore adaptation following exposure to uranium-spiked sediments[END_REF][START_REF] Martin | Extrapolating ecotoxicological effects from individuals to populations: a generic approach based on Dynamic Energy Budget theory and individual-based modeling[END_REF]. A DEB-IBM has been developed to describe three-spined stickleback population dynamics in mesocosms under control conditions [START_REF] David | Modelling population dynamics in mesocosms using an individual-based model coupled to a bioenergetics model[END_REF]. The mesocosm experimental design considered in the model was the same as the one used in de [START_REF] De Kermoysan | Effects of bisphenol A on different trophic levels in a lotic experimental ecosystem[END_REF]. Thus, this model can be adapted to predict a toxicant exposure and the potential impacts on the biotic and abiotic compartments in order to improve the understanding of the population responses found in the artificial streams. Hence, this would help to validate the modelling approach for predicting toxicant effects which is increasingly needed in risk assessment [START_REF] Galic | Potential application of population models in the European ecological risk assessment of chemicals II: Review of models and their potential to address environmental protection aims[END_REF], Hommen et al. 2010[START_REF] Forbes | A Framework for Predicting Impacts on Ecosystem Services From (Sub)Organismal Responses to Chemicals[END_REF] Accordingly, we aimed at integrating BPA effects in the DEB-IBM developed by [START_REF] David | Modelling population dynamics in mesocosms using an individual-based model coupled to a bioenergetics model[END_REF] to (i) assess the potential of modelling to predict impacts on population viability using individual level ecotoxicological endpoints and (ii) provide insights on the mechanisms of BPA toxicity during the artificial stream experiment. Thus, both direct and indirect effects of BPA on three-spined sticklebacks were incorporated in the model. Direct effects were integrated via dose-response curves on physiological processes of fish whereas indirect effects were integrated by using prey abundances measured in the exposed artificial streams. Finally, the model predictions were compared to the experimental fish population data observed in the mesocosm experiment in order to assess the accuracy of the predictions and conclude on BPA impacts at the population level.

Material and methods

Mesocosm experiment

A complete detailed description of the experimental design and the data is given in de [START_REF] De Kermoysan | Effects of bisphenol A on different trophic levels in a lotic experimental ecosystem[END_REF]. Briefly, the experiment was performed using twelve artificial streams of 20-m long and 1-m wide, locating in the North of France (INERIS, Verneuil-en-Halatte, France). Three nominal concentrations of BPA were tested (1, 10 and 100 µg/L) with 3 replicates as well as three controls. Before exposure, each mesocosm was set up with sediments, macrophytes and macroinvertebrates between October and December 2011. On March 5 th , 2012, 15 female and 10 male sticklebacks (after called the founders) with similar length were introduced in each artificial stream (mean lengths of 42.95 ± 1.40 mm and 44.36 ± 2.50 mm, respectively for females and males). Exposure began on April 15 th , 2012 and ended on September 28 th , 2012. Treatment for each mesocosm was randomly chosen. BPA water concentrations were monitored every month at different locations (0, 5 and 19 m from the inlet of water). Water temperature was measured every 10 minutes using temperature sensors (HOBO0257, Prosensor, Amanvilliers, France). Abundance and diversity of macroinvertebrates and zooplankton communities were sampled every four weeks ( details on the methods used for temperature measurements and prey samplings are provided in [START_REF] David | Modelling historical mesocosm data: Application of a fish bioenergetics model in semi-natural conditions[END_REF]).

Mesocosms were emptied on October 1 st , 2012 and, all the fish of the 12 populations were recovered, measured and weighed. Fish with a length inferior to 26 mm were assumed to be immature and thus classified as juveniles. Sexes of fish with a length superior to 26 mm were determined by visual observations of the gonads. Descriptive variables of the populations were the population abundance (N.tot), the number of female and male founders at the end of the experiment (N.F.C00, N.M.C00), female and juvenile frequencies (F.F and F.J) and mean lengths and the coefficients of variation of the lengths of five categories of individuals: male (L.M.C00 and CV.M.C00) and female founders (L.F.C00 and CV.F.C00), males (L.M.CXX and CV.M.CXX) and females born in the mesocosms (L.F.CXX and CV.F.CXX) and juveniles (L.J, CV.J). Finally, observations of the maturity of male gonads(immature vs mature) were also made and the frequency of mature males (F.M.m) was calculated as well as the mean length of mature males and the coefficient of variation of their lengths (L.M.m and CV.M.m).

Model description

The DEB-IBM predicting the population dynamics of three-spined sticklebacks in mesocosms in control conditions was already published in [START_REF] David | Modelling population dynamics in mesocosms using an individual-based model coupled to a bioenergetics model[END_REF]. Thus, a full detailed description of the model according to the ODD protocol [START_REF] Grimm | The ODD protocol A review and first update[END_REF]) is provided there. However, all parameter values of the DEB-IBM in control conditions can be found in the supplementary Information (SI) (Table S1). The model in control conditions was calibrated on two mesocosm experiment datasets and was successfully evaluated on three others: all the prediction evaluations and datasets are detailed in [START_REF] David | Modelling population dynamics in mesocosms using an individual-based model coupled to a bioenergetics model[END_REF]. In particular, the control populations of the BPA mesocosm experiment were amongst the validation datasets.

Here, we adapted the model in control conditions to predict the stickleback population dynamics in the BPA mesocosm experiments by integrating the exposure scenarios and BPA direct and indirect impacts on sticklebacks. Briefly, the model includes relevant information of the environment such as temperature, photoperiod and prey density over time. The stickleback preys considered in the model are zooplankton, higher crustacea (Gammarus pulex and Asellus aquaticus) and a diptera family (Chirominidae). The inputs of the model for the simulations are the water temperature, the photoperiod and the food scenarios (given per day) as well as the lengths of the 15 female and 10 male founders that were introduced in each artificial stream.

The time step of the model was set to 1 day and the simulation runs lasted the same amount of time as the real BPA experiment (211 days). The model outputs were the endpoints which were determined for the experimental populations. The model was implemented under Netlogo 6.0.2 (language Logo). The script of the DEB-IBM is provided in SI.

Integration of the BPA exposure and ecological scenarios

Mean BPA concentrations for each treatment and sampling date were calculated using the measured concentrations made at the three sampling locations. Then, the daily BPA concentrations were linearly interpolated between the sampling dates. Thus, three BPA concentration scenarios were estimated, one for each BPA treatment (1, 10 and 100 µg/L) which gave the daily BPA concentration in water. These BPA scenarios were used as an additional input in the DEB-IBM.

The food scenarios were calculated using the measured abundance of macroinvertebrates and zooplankton monitored in each BPA treatment (averaged on the replicates), which was then converted into energy data and linearly interpolated between the sampling dates as described in [START_REF] David | Modelling historical mesocosm data: Application of a fish bioenergetics model in semi-natural conditions[END_REF]. Thus, changing the food scenario input of the model allowed to assess the indirect BPA effects on the prey density.

As no difference of the daily mean water temperature was found between the BPA treatments (de [START_REF] De Kermoysan | Effects of bisphenol A on different trophic levels in a lotic experimental ecosystem[END_REF], the temperature scenarios were not changed in the DEB-IBM between the BPA treatments and the control.

Integration of the direct BPA impacts on sticklebacks

Direct effects of BPA on the physiological processes of three-spined sticklebacks were integrated in the model using dose-response equations. To do that, BPA effects on fish were identified at the individual level by reviewing the literature (See table S2 in SI). Only the studies using an exposure via water and reporting measured concentrations in water were kept. Only one article was kept even if the measured concentrations were not provided as this study provided toxicological information directly measured on three-spined sticklebacks [START_REF] Jolly | Detection of the anti-androgenic effect of endocrine disrupting environmental contaminants using in vivo and in vitro assays in the three-spined stickleback[END_REF]).

Based on literature and the available data, four physiological processes of fish were impaired by BPA exposure: growth, gonad maturity, eggs and embryo survival and male reproductive behavior.

As shown in [START_REF] Lindholst | Metabolism of bisphenol A in zebrafish (Danio rerio) and rainbow trout (Oncorhynchus mykiss) in relation to estrogenic response[END_REF], the kinetics in the organism of BPA proceeds very fast and the steady state is reached within 24 h for zebrafish and rainbow trout. Consequently, we directly linked the BPA concentration in water to the effects on sticklebacks in the model as the internal concentration was assumed to be proportional to the external concentration.

Finally, as the end of the BPA exposure was very close to the emptying date of the mesocosms (three days difference), we did not introduce delayed effects of BPA between the end of the BPA exposure and the emptying date.

Growth inhibition

Several studies showed that BPA could affect growth of the F0 and F1 generations of several fish species such as medaka, fathead minnow, sheepshead minnow and rainbow trout [START_REF] Sohoni | Reproductive Effects of Long-Term Exposure to Bisphenol A in the Fathead Minnow (Pimephales promelas)[END_REF][START_REF] Warner Kara | Effects of 17α-ethinylestradiol and bisphenol a on vertebral development in the fathead minnow (Pimephales Promelas)[END_REF][START_REF] Yokota | Effect of bisphenol a on the early life stage in Japanese medaka (Oryzias latipes)[END_REF][START_REF] Sadoul | Bisphenol A in eggs causes development-specific liver molecular reprogramming in two generations of rainbow trout[END_REF]. Two dose-response curves were calculated: the first one for individuals exposed only at the adult stage such as the founders fish and the second for individuals that were exposed during their entire life cycle. BPA impacts on growth were introduced by decreasing the energy conductance 𝑣̇ in the DEB model leading to a reduction of the amount of reserves mobilized as suggested by [START_REF] Sadoul | Prediction of long-term variation in offspring metabolism due to BPA in eggs in rainbow trout using the DEB model[END_REF].

Embryo and larvae mortality

First, we supposed that BPA had no impact on the number of produced eggs by the founder females. Indeed, BPA concentrations from 1.75 to 5.00 µg/L delayed ovulation in brown trout but had no effect on the quality of eggs (egg mass and fertility) [START_REF] Lahnsteiner | Effect of bisphenol A on maturation and quality of semen and eggs in the brown trout, Salmo trutta f. fario[END_REF].

Furthermore, in a long-term experiment (444 days), [START_REF] Staples | Early life-stage and multigeneration toxicity study with bisphenol A and fathead minnows (Pimephales promelas)[END_REF] showed that the number of eggs per female and per day was impacted at BPA concentrations higher than 640 µg/L for both F0 and F1 fathead minnow reproducing adults. In another long-term experiment (164 days), no statistically significant changes in number of eggs and spawns were found for adult fathead minnow at BPA concentrations ranging from 1 to 640 µg/L [START_REF] Mihaich | Adult fathead minnow, Pimephales promelas, partial life-cycle reproductive and gonadal histopathology study with bisphenol A[END_REF] although induction of vitellogenin was found for both males and females starting from 64 µg/L. Finally, [START_REF] Sohoni | Reproductive Effects of Long-Term Exposure to Bisphenol A in the Fathead Minnow (Pimephales promelas)[END_REF] found a statistical decrease of the total number of eggs spawned per female fathead minnows after 164 days of exposure for a BPA concentration of at least 1280 µg/L. Nevertheless, BPA can impact the percentage of viable eggs and larvae when eggs and embryos were directly exposed to BPA as shown for zebrafish [START_REF] Saili | Neurodevelopmental Low-dose Bisphenol A Exposure Leads to Early Life-stage Hyperactivity and Learning Deficits in Adult Zebrafish[END_REF]. Furthermore, the percentage of viable eggs and larvae can also be impacted by BPA when only adults were exposed (not the eggs directly) from 20 µg/L as shown for Gobiocypris rarus and zebrafish [START_REF] Fan | Maternal Bisphenol A exposure impaired endochondral ossification in craniofacial cartilage of rare minnow (Gobiocypris rarus) offspring[END_REF][START_REF] Guo | Parental co-exposure to bisphenol A and nano-TiO2 causes thyroid endocrine disruption and developmental neurotoxicity in zebrafish offspring[END_REF]. Finally, in two long-term experiments where fathead minnow reproducing adults and offsprings were both exposed to BPA, the hatching success was reduced at 160 µg/L [START_REF] Staples | Early life-stage and multigeneration toxicity study with bisphenol A and fathead minnows (Pimephales promelas)[END_REF]) and 640 µg/L [START_REF] Sohoni | Reproductive Effects of Long-Term Exposure to Bisphenol A in the Fathead Minnow (Pimephales promelas)[END_REF]. Here, only studies where adults and offsprings were both exposed to BPA, as in our mesocosms, were kept. We modeled this impact by decreasing the number of eggs and larvae leaving the nest.

Male reproductive behavior

One study on three-spined sticklebacks showed an inhibition of the spiggin production [START_REF] Jolly | Detection of the anti-androgenic effect of endocrine disrupting environmental contaminants using in vivo and in vitro assays in the three-spined stickleback[END_REF] for BPA concentrations ranging from 2.5 pg/L to 250 µg/L. Indeed, in natural conditions, this protein is only secreted during the breeding season by male sticklebacks to use it as a glue during nest building [START_REF] Wootton | A Functional Biology of Sticklebacks[END_REF]. As the DEB-IBM does not directly model the spiggin production, we extrapolated this BPA effect by reducing the length of the breeding period. Thus, the number of built nests will decrease as males will stop reproductive processes sooner.

Gonad malformations

In the mesocosm experiment, malformations of the gonads (decrease of the ratio between the length of the gonads and the standard fish length) were observed for both male and female sticklebacks (de [START_REF] De Kermoysan | Effects of bisphenol A on different trophic levels in a lotic experimental ecosystem[END_REF]. A decrease of the gonadosomatic index (GSI) due to BPA exposure for both males and females have also been found in two studies [START_REF] Sohoni | Reproductive Effects of Long-Term Exposure to Bisphenol A in the Fathead Minnow (Pimephales promelas)[END_REF][START_REF] Mihaich | Adult fathead minnow, Pimephales promelas, partial life-cycle reproductive and gonadal histopathology study with bisphenol A[END_REF]. We introduced this effect on the energy allocated to sexual maturity and reproduction (predicted by the DEB model). However, we did not integrate this effect for the founders as they already had mature gonads.

Calibration of the dose-response parameters

The stress level for each type of effect integrated in the DEB-IBM was determined with a Hill equation (Eq. 1).

𝑠 (𝑡) = C (𝑡) 𝑛 𝑡𝑜𝑥 𝐸𝐶 50 𝑛 𝑡𝑜𝑥 + C (t) 𝑛 𝑡𝑜𝑥 (Eq. 1)
With 𝐶(𝑡) the exposure concentration at time t, 𝐸𝐶 50 the concentration leading to 50 % of effect, 𝑛 𝑡𝑜𝑥 the Hill coefficient and 𝑠(𝑡) the stress level at time t. 𝑛 𝑡𝑜𝑥 and 𝐸𝐶 50 values are specific to the type of BPA effect integrated in the model.

First, these dose-response parameters were fitted on the literature data using the R package "drc" [START_REF] Ritz | Dose-Response Analysis Using R[END_REF].

Second, to assess if the main BPA impacts were integrated in the DEB-IBM, the doseresponse parameters were recalibrated on the observed experimental population data with an iterative process using a genetic algorithm available within the software BehaviorSearch. The distance to minimize was the weighted sum of least squares (WLS) between the predictions and the observations of the endpoints of the population (Eq. 1).

𝑊𝐿𝑆 = ∑ (𝑆𝑖𝑚 𝑖 -µ 𝑜𝑏𝑠 𝑖 )² 𝜎 𝑜𝑏𝑠 𝑖 ² 𝑛 𝑖=1 (Eq. 2)
With 𝑛 the number of observed outputs, 𝑆𝑖𝑚 𝑖 the value of the simulated output 𝑖, µ 𝑜𝑏𝑠 𝑖 the mean of the observed output 𝑖 and 𝜎 𝑜𝑏𝑠 𝑖 the standard deviation of the observed output 𝑖.

Sensitivity analyses of the dose-response parameters

Sensitivity analysis by using the variance-based Sobol' method [START_REF] Sobol | Estimating the approximation error when fixing unessential factors in global sensitivity analysis[END_REF][START_REF] Saltelli | Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index[END_REF] were performed on the dose-response parameters varying 10 % of their values for the highest exposure scenario. The Sobol's method identifies the contribution of each parameter (First order index) as well as the interactions between parameters (Total index). To assess the impacts of parameter variations over time of the experiment, three sensitivity analyses were made at different time of simulations: during the beginning of the breeding period (70 days), during the end of the breeding period (140 days) and at the end of the experiment (211 days).

The sensitivity analyses were made with the "sobol2007" function implemented in the R package "sensitivity" [START_REF] Iooss | Global Sensitivity Analysis of Model Outputs: R package 'sensitivity[END_REF]) which uses the Mauntz estimators recommended for small first-order and total indices.

Model simulations

Predictive capacity of the model using the dose-responses fitted on literature data

To assess the predictive capacity of the model, population endpoints were simulated with the DEB-IBM using the individual dose-responses fitted on literature data for each BPA treatment.

1000 simulations were performed to account for the stochasticity of the model: inter-mesocosm environmental variability, inter-individual performance variability, and stochastic processes (for example survival). Simulation results were then compared to the observed experimental data. The relative errors (RE) and normalized root-mean-square deviations (NRMSD) between the simulations and the observations in the artificial streams were calculated to conclude on the predictive capacity of the model.

Explanation of the observed effects at the population level

Before integrating BPA effects into the model, population endpoint data from each treatment were compared to the simulated or observed control endpoints. The LOECs were calculated in two different ways using the methodology described in [START_REF] David | Modelling population dynamics in mesocosms using an individual-based model coupled to a bioenergetics model[END_REF]. First, they were calculated using the observed control endpoints in the artificial streams as a reference and a Dunnett's post hoc test as in de [START_REF] De Kermoysan | Effects of bisphenol A on different trophic levels in a lotic experimental ecosystem[END_REF]. Second, they were estimated comparing the distribution of the control endpoints estimated by the DEB-IBM to the observations made in the exposed mesocosms using a Kolmogorov-Smirnov test. The level of significance for all tests was 5 %. To estimate the distribution of the control endpoints, 1000 simulations of the DEB-IBM were performed.

Then, in order to gain more knowledge on the observed BPA impacts, 1000 simulations of the DEB-IBM were made for each BPA treatment with the recalibrated dose-responses on the observed population data in the artificial streams. The simulation results were compared to the simulations using the dose-responses fitted on literature data to evaluate the added value of this recalibration by calculating the RE and NRMSD between the simulations and the observations in the mesocosms.

Finally, the contribution of either direct or indirect BPA effects was assessed with the DEB-IBM with the recalibrated dose-response curves. Four types of simulations were made for the 100 µg/L BPA treatment: (i) without any BPA effects, (ii) integrating only indirect effects of BPA, (iii) integrating only the direct effects of BPA on fish and finally, (iv) integrating both direct and indirect effects. The RE and NRMSD for each type of simulation were compared to conclude on the contribution of direct and indirect effects. Furthermore, the ratios between the root-mean-square deviations of the simulations with the different BPA effects and the simulations without BPA effects for each output were calculated to conclude on the added benefit of the models.

Populational dose-response curves

Population dose-response curves were extrapolated between the 0 and 100 µg/L BPA treatments. To do that, the food scenarios were linearly interpolated between the actual tested BPA concentrations. The chosen extrapolated concentrations were 0.5, 1.5 and 2 µg/L and from 5 to 70 µg/L with a step of 5 µg/L. 1000 model simulations were made for each extrapolated concentration. The population dose-response curves were calculated by linearly interpolated between the median of the simulations per concentrations. They were then compared to the ones calculated on the observed data using a Hill model available in the R package "drc" [START_REF] Ritz | Dose-Response Analysis Using R[END_REF]).

Results

BPA exposure and ecological scenarios

The mean concentration throughout the experiment for each treatment were 0.97, 4.53 and 76.6 µg/L respectively for the 1, 10 and 100 µg/L BPA treatments. Figure 1 presents the three BPA exposure scenarios for the different treatments. Globally, BPA concentrations in water were lower than expected for the 10 and 100 µg/L treatments (between 46.5 and 91.8 % of the expected concentration for the 100 µg/L treatment and between 19.5 and 80.0 % for the 10 µg/L treatment). For the 1 µg/L treatment, concentrations were either lower or higher than the expected values (between 39.7 and 148.9 %).

Concerning the food scenarios, Figure 2 presents the total energy calculated using all stickleback preys considered in the model and Figure S1 in SI presents the energy for each prey type. The peak of total energy (around 110 days of the experiment) of prey observed in the control and the 1 µg/L BPA treatment was lower in the 10 µg/L treatment and did not occurred in the 100 µg/L BPA treatment (Figure 2). As shown in Figure S1, this result was mainly due to the dynamics of the macroinvertebrates.

Integration of the direct BPA effects and sensitivity analyses

Figure 3 presents the dose-responses fitted on literature data or on population observed data in the artificial streams for the four integrated BPA impacts in the DEB-IBM and Table 1 presents the dose-response parameter values. Furthermore, the dose-response curve fitted on literature data for the number of produced eggs by females showed no effect of BPA at 100 µg/L (Figure S9). Thus our assumption of neglecting this effect sounds justified. Figure 4 shows the results of the Sobol sensitivity analysis of the different types of outputs: female frequency, frequency of mature males, mean length of founder females and juveniles. The results for the other outputs can be found in SI Figure S4. Overall, the Sobol's total indices were important meaning that the interactions between the toxic parameters were strong.

The EC50 parameters of the effect on the duration of the breeding period had strong effects on the female frequency (Figure 4a). However, for the frequency of mature males, the parameter that had the most influence was the EC50 of the effect on the growth of non-founder fish (Figure 4b). Similarly, the juvenile mean length was strongly impacted by a variation of the EC50 and the Hill coefficient parameter values relative to the duration of the breeding period (Figure 4d).

However, the parameters relative to growth of founders strongly impacted the mean length of founder females (Figure 4c). Thus, globally all toxic parameters had strong influence on all outputs, except the parameters relative to the effect on egg survival which had low impacts (Figures 4 andS4).

Furthermore, the influence of the parameters was dependent on the simulation time (Figure S5). For example, the parameters related to effects on growth strongly influenced the mean length of juveniles at the beginning and middle of the experiment (at 70 and 140 days) whereas those related to a reduction of the duration of the breeding period had more influence at the end of the experiment.

Model predictive capacity

The comparison of the model simulations using the dose-response curve fitted on literature data with the mesocosms observations for each BPA treatment can be found in Figure S2. The model gave relatively good predictions of the observed population endpoints (RE = 26 % and NRMSD = 0.31 for all endpoints in the 100 µg/L BPA treatment). Nevertheless, the decrease of the frequency of mature males in the 100 µg/L BPA treatment was not predicted.

Furthermore, the mean length of females and males born in the mesocosms were slightly overestimated. For the other endpoints, the model predicted well the population dynamics under BPA stress. Concerning the length frequencies, the model gave good predictions for the control and the 1 µg/L treatment (Figure S3a,b). However, the DEB-IBM overestimated the length frequencies of the fish for the first cohorts for the 10 µg/L and 100 µg/L treatments at the end of the experiment (Figure S3c,d).

Explanation of the observed effects at the population level

Table S3 presents the comparison of the LOECs using the observed control endpoints in the artificial streams or using the distribution of the control endpoints estimated by our model as references for the statistical tests. Both approaches identified a LOEC at 100 µg/L for the juvenile mean length and the CV of their lengths. The LOEC for the frequency of mature males (F.M.m) was decreased to 10 µg/L and LOECs equal to 10 or 100 µg/L were found for four other endpoints (L.F.C00, L.F.CXX, L.M.CXX and CV.F.CXX) using the distribution of the control endpoints estimated by our model. All other endpoints had a LOEC above 100 µg/L.

Consequently, our model showed that seven endpoints were found to have a LOEC at 10 or 100 µg/L compared to three using only the observed mesocosm data.

After recalibration of the toxic parameters (see Table 1), the dose-response curves for the growth and the effects on gonads were shifted to the right which implies an increase of the BPA stress levels at the studied concentrations (Figure 3). No major differences were found for the effects on eggs and larvae survival and spiggin inhibition. This recalibration improved all the predictions of the endpoints (RE = 14 % and NRMSD = 0.22 for all endpoints in the 100 µg/L BPA treatment) (Figure 5), and especially the one concerning the frequency of mature males.

Finally, the length frequencies were still well predicted for the control and 1 µg/L treatment (Figure 6 a,b) and were strongly improved for the 100 µg/L BPA treatment (Figure 6d).

However, predictions for the 10 µg/L BPA treatment could still be improved as fish lengths of the first cohort were still slightly overestimated (Figure 6c). Figure 8 presents the benefit of the model considering the indirect and/or direct BPA effects compared to the one without BPA effects. First, the average benefit on all the population endpoints is 20.7 % with the model using both type of BPA effects against 16.5 % and 6.1 % for the model using direct or indirect effects alone respectively. Indeed, the benefit of the model using indirect effects alone was lower than the one using direct effects alone for most endpoints.

For example, almost no indirect effects could explain the frequency of mature males (0.86 %)

whereas for other endpoints like the female frequency, the benefit of the model considering the indirect effect alone was not negligible (25.8 %). Finally, the benefit of the model, regardless of the BPA effects, were negative for two endpoints (total abundance and the CV of the length of mature males).

Population dose-response curves

Figure S6 in SI presents the extrapolated food scenarios and Figure 9 the extrapolated doseresponse curves at the population level with the DEB-IBM compared to the dose-responses calculated on the observed data for the number of female founders, the CV of the length of the males born in artificial streams and the frequency of mature males. Figure S7 presents the results for all population endpoints. First, the model showed no impacts on total abundance and the number of male and female founders at the end of the experiment whereas the doseresponses fitted on the observed data suggested an effect (Figure 9a and S7). On the contrary, no effect on the juvenile frequency was shown with the dose-responses fitted on the observed data whereas the model showed one (Figure S7). For the juvenile mean length and the frequency of mature males, the shape of the extrapolated dose-response curves with the DEB-IBM showed a less noticeable effect between the 10 and 100 µg/L treatments compared to the dose-responses calculated on the observed data (Figure 9b and S7). Furthermore, dose-response curves having a slope which changed direction within the range of the tested treatments (i.e. non-monotonic dose-responses) were found for the CV of the lengths of both males and females born in mesocosms as well as the CV of the lengths of the mature males (Figure 9c and S7). Finally, for the other endpoints, the extrapolated dose-responses with the DEB-IBM presented the same shape than the ones calculated on the observed data (Figure S7).

Discussion

Population modelling in risk assessment was suggested to be a powerful tool to gain more knowledge on the mechanisms of chemical toxicity on populations [START_REF] Forbes | Adding Value to Ecological Risk Assessment with Population Modeling[END_REF]) and to extrapolate from ecotoxicological observations to relevant ecological effects [START_REF] Forbes | The extrapolation problem and how population modeling can help[END_REF]). Hence, it would help to improve decision-making and ecological management [START_REF] Goussen | Integrated presentation of ecological risk from multiple stressors[END_REF]). Here, we proposed a DEB-IBM to predict BPA effects on the population dynamics of three-spined sticklebacks in artificial streams. This model was based on a DEB-IBM previously developed for predicting population dynamics of three-spined sticklebacks without toxicant stress [START_REF] David | Modelling population dynamics in mesocosms using an individual-based model coupled to a bioenergetics model[END_REF].

Integration of the direct BPA effects on sticklebacks

The development of the individual dose-responses describing BPA effects was not straightforward to handle due to the heterogeneity of the literature studies (different studied species and exposure design). The exposure route via food was neglected in our model, while some studies have showed an accumulation of BPA in invertebrates [START_REF] Nurulnadia | Accumulation of endocrine disrupting chemicals (EDCs) in the polychaete Paraprionospio sp. from the Yodo River mouth, Osaka Bay, Japan[END_REF] suggesting that this exposure route could be possible for sticklebacks. However, the only study on fish that tested the effects of contaminated food with BPA [START_REF] Drastichová | Effect of Exposure to Bisphenol A and 17ß-estradiol on the Sex Differentiation in Zebrafish (Danio rerio)[END_REF]) showed effects at higher concentrations than those studied in the mesocosm experiment, thus the results

were not considered.

From literature, four main BPA effects on fish appeared to be interesting to include in the model (effects on growth, gonad maturity, eggs and embryo survival and male reproductive behavior). However, the BPA effects on the reproductive behavior of males were more complex to incorporate as very scarce literature data was available on the mating behavior of fish. One study showed a decrease of male locomotion and behavior during courtship in zebrafish [START_REF] Li | Behavioural effect of low-dose BPA on male zebrafish: Tuning of male mating competition and female mating preference during courtship process[END_REF]) at a low-dose BPA concentration but, to our knowledge the measured concentrations in water were not available. We thus extrapolated this effect by using the study on three-spined sticklebacks which showed an inhibition of the spiggin, protein used for nest building [START_REF] Jolly | Detection of the anti-androgenic effect of endocrine disrupting environmental contaminants using in vivo and in vitro assays in the three-spined stickleback[END_REF]). In the model, spiggin was correlated to the breeding period (inhibition of spiggin implies a reduction in the duration of the breeding season) which would thus decrease the success rate of the reproductive processes of males (in particular a decrease of the number of nests) [START_REF] Mori | The Breeding System of the 3-Spined Stickleback, Gasterosteus-Aculeatus (Forma-Leiura) with Reference to Spatial and Temporal Patterns of Nesting Activity[END_REF]. Other DEB-IBM parameters could have been impacted instead of the duration of the breeding period, like the nest mortality or the probability of stopping the reproduction processes. However, impacting the breeding period allows to take in account different possible physiological and behavioral effects such as a reduction of the overall courtship and nest building. Furthermore, this parameter influenced the most the outputs of the model as shown in [START_REF] David | Modelling population dynamics in mesocosms using an individual-based model coupled to a bioenergetics model[END_REF].

Globally, the DEB-IBM with the toxic parameters fitted on literature data at the individual level gave overall good predictions of the population endpoints, especially for the mean length and CV of the lengths of juveniles as well as the female frequency. Hence, this supports the idea that individual-based modelling may be a sound approach for predicting toxicant effects at the population level from the toxicological data recorded at the individual level [START_REF] Galic | Potential application of population models in the European ecological risk assessment of chemicals II: Review of models and their potential to address environmental protection aims[END_REF], Hommen et al. 2010[START_REF] Forbes | A Framework for Predicting Impacts on Ecosystem Services From (Sub)Organismal Responses to Chemicals[END_REF].

Explanation of the observed effects at the population level

The sensitivity analyses made for the 100 µg/L BPA treatment revealed that the influence of the toxic parameters was strongly dependent on the studied outputs and the time of the simulation (Figures 4,S4 and S5). First, the mean length of juveniles was mostly impacted by the variations of the parameter values relative to a reduction of the duration of the breeding period whereas the male and female mean lengths were most sensitive to changes in the toxic parameters related to growth. Furthermore, the female and juvenile frequencies and the total fish abundance (Figures 4a and S4 a,b) were mostly impacted by a variation of the toxic parameters related to a reduction of the duration of the breeding season. Indeed, as the last reproductive events were inhibited, less juveniles were present in the artificial streams at the end of the experiment which directly impacted the juvenile frequency as well as the mean length of the juveniles. Shortened reproductive season was already shown to impact cohort abundances with another IBM for zebrafish [START_REF] Hazlerigg | Population relevance of toxicant mediated changes in sex ratio in fish: An assessment using an individual-based zebrafish (Danio rerio) model[END_REF]) and the duration of the breeding period was already one of the critical parameter for the DEB-IBM developed without toxicant stress [START_REF] David | Modelling population dynamics in mesocosms using an individual-based model coupled to a bioenergetics model[END_REF]. Finally, the sensitivity of the population endpoints was clearly dependent on the type of BPA effects. Thus, it was necessary to include multiple BPA effects to accurately predict all population endpoints which highlights the complexity of mechanisms of toxicity that could have occurred in the mesocosm experiment.

The DEB-IBM simulations of the variables of the stickleback populations for each BPA treatment with the recalibrated dose-response curves on population data were close to the observed endpoints. Thus, this suggests that all main BPA effects have been integrated

in the model and explained the population structure found in the 100 µg/L treatment. However, some predictions, such as the length frequency of the fish of the first cohort, could still be improved, especially for the 10 µg/L BPA treatment, suggesting that other effects could have happened in the artificial streams, and consequently decreased the growth rate of the first cohort. For example, BPA was shown to have complex neurobehavioral effects. Indeed, it altered the locomotion and can induce hypoactivity as shown on fish larvae at concentrations between 2 and 200 µg/L [START_REF] Inagaki | Low dose exposure to Bisphenol A alters development of gonadotropin-releasing hormone 3 neurons and larval locomotor behavior in Japanese Medaka[END_REF][START_REF] Guo | Parental co-exposure to bisphenol A and nano-TiO2 causes thyroid endocrine disruption and developmental neurotoxicity in zebrafish offspring[END_REF]. BPA was also shown to reduce the aggressive behavior of adult zebrafish as well as induce a weaker adaptability to new environments [START_REF] Wang | Chronic bisphenol A exposure alters behaviors of zebrafish (Danio rerio)[END_REF]. On the contrary, [START_REF] Saili | Neurodevelopmental Low-dose Bisphenol A Exposure Leads to Early Life-stage Hyperactivity and Learning Deficits in Adult Zebrafish[END_REF] showed that low-dose BPA exposure levels led to early zebrafish life-stage hyperactivity and adults exhibited learning deficits. Thus, BPA effects on fish locomotion and activity are not fully understood but could impact the predation processes and growth density-dependence driven by fish interactions.

Contribution of the indirect and direct BPA effects

The comparison of the DEB-IBM simulations integrating direct and/or indirect effects showed that we obtained simulations closer to the observed endpoints when both direct and indirect impacts were taking into account for the 100 µg/L BPA treatment (Figure 7). Indeed, as seen in Figure 2, total energy of stickleback preys was lower at this treatment. This result was not surprising as BPA impacts on invertebrates have been widely described (Flint et al. 2012, Canesi and[START_REF] Canesi | Environmental Effects of BPA: Focus on Aquatic Species[END_REF]. In particular for the stickleback main preys, BPA was shown to affect growth and feeding of Asellus aquaticus during acute and chronic exposures to high BPA concentrations [START_REF] Plahuta | Adverse effects of bisphenol A on water louse (Asellus aquaticus)[END_REF]. Reproduction and mortality of Gammarus fossarum was also impacted starting from 5 µg/L of BPA [START_REF] Ladewig | Effects of bisphenol A on Gammarus fossarum and Lumbriculus variegatus in artificial indoor streams[END_REF], Schirling et al. 2006).

Developmental inhibition induced by BPA exposure was observed at very low exposure concentrations (10 ng/L) on Chironomus riparius larvae [START_REF] Watts | Exposure to 17α-ethinylestradiol and bisphenol Aeffects on larval moulting and mouthpart structure of Chironomus riparius[END_REF]. Finally, zooplankton species such copepods were also impacted at low BPA concentrations at 1 and 10 µg/L [START_REF] Marcial | Estrogenic compounds affect development of harpacticoid copepod Tigriopus japonicus[END_REF].

Assessing the impacts of species loss in an ecosystem requires to understand the processes that occur across trophic levels in order to consider possible cascading effects [START_REF] Duffy | Guidance on tiered risk assessment for plant protection products for aquatic organisms in edge-of-field surface waters[END_REF][START_REF] O'connor | Distinguishing between direct and indirect effects of predators in complex ecosystems[END_REF]). Here, as explained in de Kermoysan et al. ( 2013), indirect effects could explain the shift of the population structure towards higher lengths at the 100 µg/L BPA treatment. The model predictions of the mean length of juveniles were for example improved when including only the food scenario for the 100 µg/L treatment (Figure 7). The growth of the individuals could be reduced as less food was available leading to a larger number of fish in the length class of [20; 26 mm[ that did not reach the adult stage. Another explanation is that density-dependent mechanisms for growth were lower as the decrease of the food resource indirectly also reduced fish reproduction. As seen in Figure 8, the benefit of the model including indirect effects alone was 17.7 % which was not negligible compared to the benefit including direct effects alone (31.5 %).

On the contrary, the decrease of the frequency of mature males is mainly explained by direct BPA effects (Figure 8). Indeed, the decrease of the food resource in the environment could have explained a decrease of the maturity as described in the DEB model [START_REF] Kooijman | Dynamic Energy Budget for metabolic organization[END_REF]) but the results showed that it was clearly not sufficient to explain the observed effect. Concerning total abundance (Figure 8), the benefits of the models including the different types of BPA effects were negative. However, the large coefficient of variation of the observations (26.7 % for the 100 µg/L treatment) induces great uncertainty and then, the benefit may probably be smaller than the noise. Indeed, as seen in Figure 5, it is not possible to conclude on the relevance of the model predictions for total abundance as the three replicates were in the predictive interval of the simulations for the 100 µg/L treatment.

Globally, the combination of both direct and indirect effects gave better predictions (Figures 7 and8), highlighting the complex interactions between the trophic levels. Thus, this BPA case study highlights that modelling food availability as well as the effects on the physiology of the organisms would help to provide relevant chemical exposure-response relationships at the population level which is necessary for ecological risk assessment [START_REF] Forbes | A Framework for Predicting Impacts on Ecosystem Services From (Sub)Organismal Responses to Chemicals[END_REF].

Extrapolation of the ecotoxicological effects at the population level

In risk assessment, the use of the DEB theory [START_REF] Kooijman | Dynamic Energy Budget for metabolic organization[END_REF] coupled to an IBM was suggested to be a convenient approach to extrapolate the effect of toxicants measured on the individual life history traits to effects at higher levels of biological organization [START_REF] Beaudouin | Individual-based model of Chironomus riparius population dynamics over several generations to explore adaptation following exposure to uranium-spiked sediments[END_REF][START_REF] Martin | Extrapolating ecotoxicological effects from individuals to populations: a generic approach based on Dynamic Energy Budget theory and individual-based modeling[END_REF].

However, as the literature data concerned different fish species exposed to BPA with diverse experimental designs, this study gives an example of different extrapolation issues that could happened. More specifically, extrapolating effects between species adds uncertainties in predictions. Here, the comparison of the calibrated dose-response curves on population responses in artificial streams and the dose-response curves fitted on literature data showed that the BPA effects on growth and gonad maturity were stronger than assumed by literature data.

Furthermore, the assumption of the extrapolations between literature data and the way of integrating the effect in the DEB-IBM could also induce a discrepancy between the stress level calculated on the dose-response curve and the population observed responses. For example, the effect of BPA on the male reproductive behavior was extrapolated from the response of spiggin, an androgen and antiandrogen biomarker [START_REF] Katsiadaki | The potential of the three-spined stickleback (Gasterosteus aculeatus L.) as a combined biomarker for oestrogens and androgens in European waters[END_REF]. Nevertheless, we found that the recalibration of the dose-response parameters did not strongly change the dose-response curve of the spiggin inhibition on sticklebacks suggesting that including this effect on the duration of the breeding season was a relatively acceptable extrapolation. Finally, to go further on the modelling of the exposure, the population model could be coupled to toxicokinetic models that provide a quantitative mechanistic framework to predict the time-course of the concentration of chemicals in an organism [START_REF] Grech | Generic physiologically-based toxicokinetic modelling for fish: Integration of environmental factors and species variability[END_REF].

Furthermore, adverse outcome pathways (AOPs) have emerged as a useful framework to assess impacts of chemicals at the molecular level. AOPs describe mechanisms linking molecular initiating events to outcomes at an individual level. Here the concept could be applied to the case of sticklebacks as spiggin synthesis in sticklebacks is assumed to be regulated by a renal androgen receptor. Indeed, the inhibition of the spiggin production could occur if antiandrogens occupy the androgen receptor without activating it. Thus, spiggin inhibition here is a good example of how AOPs could be used to identify mechanisms of toxicity and extrapolate them from the molecular to the population level with a modelling approach.

The predictive capacity of our model was assessed on one reproductive season and thus, should be also evaluated for more than one season. To do that, the ecological scenarios would need to be extrapolated. This could be done by including a specific model of the aquatic food webs describing the dynamics over time. Furthermore, BPA effects on reproduction can be supposed to be long lasting or even irreversible as seen for example with the alterations of the gonadotroph cells in zebrafish at BPA concentrations of 100 and 1000 µg/L [START_REF] Molina | Hypothalamic-pituitary-ovarian axis perturbation in the basis of bisphenol A (BPA) reproductive toxicity in female zebrafish (Danio rerio)[END_REF]). In a study on rainbow trout egg exposed to BPA, Sadoul et al. ( 2017) also showed longterm effects on growth and gene expressions related to metabolism. Consequently, to make long-term predictions of BPA effects, it would be highly interesting to consider long-lasting effects of BPA on stickleback growth as well as on gonad maturity. However, according to the parsimonious principle, as the end of the mesocosm experiment was really close to the end of the exposure, the BPA long-lasting effects were neglected in our model.

To go further, integrating landscape characteristics (e.g. habitat, vegetation patches) in addition to species interactions and environmental conditions (e.g., temperature, food availability) in the ecological scenarios could be interesting as they may also be impacted by toxicants [START_REF] Rico | Developing ecological scenarios for the prospective aquatic risk assessment of pesticides[END_REF]. Furthermore, other ecological responses, like, for example, a shift of diet to other available species in response to the absence of the preys [START_REF] Allen | Temporal Patterns in Diet and Rate of Food-Consumption of the 3-Spined Stickleback (Gasterosteus-Aculeatus L) in Llyn Frongoch, an Upland Welsh Lake[END_REF] was not considered although this could minimize the indirect BPA effects.

Added value of ecological modelling to estimate toxicant thresholds

We here analyzed the power of modelling to estimate toxicant thresholds. First, the model without including toxicant effects was used to assess the LOECs (methods presented in David et al. ( 2019)) and showed that 39 % of the monitored population endpoints had a LOEC at 10 or 100 µg/L (only 17 % based on simple statistical comparisons (de Kermoysan et al. 2013))

(see Table S3). Moreover, the minimum LOEC value was estimated at 100 µg/L in de [START_REF] De Kermoysan | Effects of bisphenol A on different trophic levels in a lotic experimental ecosystem[END_REF], whereas it was lowered to 10 µg/L using data analysis based on population dynamics modelling. Consequently, our approach improved environmental risk assessment by decreasing the risk of false negative results when analyzing data from artificial streams.

Furthermore, our study offers the opportunity to compare dose-responses at the population level extrapolated by the DEB-IBM to the dose-responses calculated on the observed data with a Hill model (Figure 9 and S7). We found that in most cases, the extrapolated dose-response curves found with the model were consistent with the ones fitted on the observed data.

However, the model helped to refine the shape of the dose-responses between the tested concentrations. For example, the dose-response estimated by the model showed a slower decrease of the frequency of mature males. Furthermore, the model showed non-monotonic dose-response relationships, as frequently found for endocrine disruptors [START_REF] Vandenberg | Regulatory decisions on endocrine disrupting chemicals should be based on the principles of endocrinology[END_REF]. This was specially the case for the CV of the lengths of females and males born in the mesocosms and the mature males. Such unconventional dose-responses represents a challenge to current risk assessment [START_REF] Lagarde | Non-monotonic dose-response relationships and endocrine disruptors: a qualitative method of assessment[END_REF]. Hence, population models could add value to ecological risk assessment, extrapolating dose-responses at population level, which is essential for improving decision-making and protection of ecosystems [START_REF] Forbes | The extrapolation problem and how population modeling can help[END_REF], Galic et al. 2010[START_REF] Forbes | Adding Value to Ecological Risk Assessment with Population Modeling[END_REF].

Conclusions

The predictive capacity of our DEB-IBM in context of risk assessment of chemicals was assessed on stickleback population data from a mesocosm experiment where the artificial streams were exposed to BPA. We showed that including BPA effects on stickleback gonad formation, growth, male reproductive behavior and eggs and larvae survival gave us a relevant predictive model of the effects at the population level. Recalibration of the dose-response curves of the BPA effects on fish on experimental population data showed that the main BPA effects were integrated in the model. Indirect effects were not negligible and participated to the stickleback population structure observed in the 100 µg/L BPA treatment. Therefore, ecological models, such as ours, have a strong potential for improving environmental risk assessment.

Figures and tables

Table 1. Description and parameter values of the dose-response relationships fitted on literature data or recalibrated on stickleback population observations in the mesocosms (artificial streams).

Figure 7

 7 Figure7shows the comparison between the model predictions including either direct and/or

Figure 1 .

 1 Figure 1. BPA exposure scenarios based on the measured concentrations made in the

Figure 2 .

 2 Figure 2. Food scenarios throughout time for each BPA treatment (green: 0 µg/L, yellow: 1

Figure 3 .

 3 Figure 3. Dose-responses of the effects of BPA on growth (a), egg and larvae mortality (b),

Figure 4 .

 4 Figure 4. Results of the Sobol sensitivity analysis for the female frequency (a), the frequency

Figure 5 .

 5 Figure 5. Simulation of the toxicity of BPA using calibrated dose-responses on

Figure 6 .

 6 Figure 6. Probabilistic distributions of the length frequency predicted by the model

Figure 7 .

 7 Figure 7. Predictions of the mean length of each category of individuals (a), the CV of their

Figure 8 .

 8 Figure 8. Added benefit of the models considering the indirect BPA effects alone (light

Figure 9 .

 9 Figure 9. Dose-response curves at the population level for the number of female founders

Acknowledgment

The authors wish to thank G. De Kermoysan who contributed to the success of the mesocosm experiments and C. Tebby for her suggestions of improvement. We also would like to thank the three anonymous reviewers for their helpful comments and suggestions of improvements. The study was supported by the French Ministry in charge of Ecology and Sustainable Development for this study, within 190 program and the European Regional Development Fund (ERDF).


Jérôme Daubech
email: jerome.daubech@ineris.fr

Emmanuel Leprette

Christophe Proust

Guillaume Lecocq

Further insight into the gas flame acceleration mechanisms in pipes. Part I: experimental work

Keywords:

 published the conclusions of a rather large experimental work devoted to the gas flame acceleration down a long pipe. It was concluded that the flame propagation could be represented by a constantly accelerating piston. The acceleration parameter seems to be primarily linked to the expansion velocity of the burnt product. Other parameters seemed of secondary importance questioning in particular the respective roles of the turbulence of the flame and of the instabilities.

), varying the diameter of the pipe (150 and 250 mm) and the reactivity of the mixtures (methane-air and hydrogen air at various equivalence ratios). The smooth pipe is transparent enabling a direct visualization of the flame during the flame propagation and a refined resolution of the flame trajectory (in the steel pipes standard flame sensors were used). The pressure was measured at various locations but also the flow velocities in the boundary to try and detect any turbulence development. Only homogeneous and quiescent mixtures were studied and the flame was propagated from a closed ignition end toward the other open end. The results of the parametric study are presented in this paper.

Introduction

Coal mining industry remains one of the most high-risk industries in the world. Gas explosion accidents often happen following a leak of methane and ignition the flammable cloud formed in the mine infrastructures. The firedamp explosions generate important damages on critical equipment, injuries and dramatical deaths of people working in mines. As an example, there were 24 events related to methane hazard (explosions and ignitions) in Polish coal mines along the last 12 years. They caused the death of 62 miners and 75 were seriously injured. Losses in mining infrastructure and equipment were also significant. To improve the safety in mines and underground infrastructures, the European EXPRO project (prediction and mitigation of methane explosion effects for improved protection of mines infrastructures and critical equipment) was performed by the GIG institute, EMAGPL, KOMPWE in Poland, FSB in Spain and INERIS in France. One of the purpose of the project is to develop knowledge about the explosion process to determine the influencing parameters. [START_REF] Proust | Gas flame acceleration in long ducts[END_REF] published an interesting review of explosion process and flame acceleration in tunnel and pipes. In this work, he exposed different mechanisms for flame acceleration:

The most widely accepted mechanism so far is the continuous increase of the turbulence of the reactive mixture induced by the expanding burnt products pushing the reactants ahead [START_REF] Borghi | Turbulent combustion modelling[END_REF][START_REF] Clarke | Fast flames[END_REF]. Due to friction at the wall, turbulence would be generated in proportion of the mean flow velocity. The burning velocity would increase inducing an increase of the expansion velocity of the burnt products, hence of the velocity of the reactants.; [START_REF] Deshaies | Flame-speed sensitivity to temperature changes and the deflagration to detonation transition[END_REF] explain thanks to a theoretical development that the gradual acceleration of the flow ahead of the flame due to burnt product expansion is produced by a series of compression waves. The temperature of the reactants ahead of the flame increases accordingly as well as the burning velocity. The flame then selfaccelerates;

A flame acceleration mechanism can also result from the flame front structuration generated by the pressure waves and their reflections on the extremities [START_REF] Daubech | Contribution à l'étude de l'effet de l'hétérogénéité d'un prémélange gazeux sur la propagation d'une flamme dans un tube clos[END_REF]. Proust observed from his experiments that the flame propagation could be represented by a constantly accelerating piston. The acceleration parameter seems to be primarily linked to the expansion velocity of the burnt product. Other parameters seemed of secondary importance questioning in particular, the respective roles of the turbulence of the flame and of the instabilities.

The experimental results presented in this paper is a follow up of the examination of flame acceleration mechanisms. This work focuses on the influence of diameter and the roughness of straight long pipes closed at one end and opened at the other hand in the case of stoichiometric methane explosion.

Experimental set-up

The experimental device is a 24 m long pipe filled with a flammable mixture (Figure 1). The flammable mixture is prepared in a 2 m 3 spherical vessel using the partial pressure method and is injected in the tube. The initial pressure in the 2 m cube vessel is about 5 bar. The tube is swept by this mixture about 5 times its volume.

Tubes are made of transparent PMMA or steel. Two diameters are tested: 150 and 250 mm (Figure 2). The tube is closed on the side from which the gas is injected and the mixture ignited. The other extremity is covered with a plastic sheet drilled with the small hole to avoid a pressure rise-up of the tube during gas injection.

The instrumentation comprises:

Three pressure sensors: o Near the ignition point, at 5.5 m and at 15.5 m from the closed end flame detection in steel pipes:

o Flame probes: photovoltaic cells which catch the light produced by the flame (at 0.5, 4.5, 10.5 and 15.5 m from the ignition source) o Ionization probes in the first meter (4.5, 45 and 90 cm from ignition source -Figure 3) flame detection in transparent pipes:

o Photron fast video camera Two Pitot probes (Proust, 2018) to measure flow velocity (Figure 3)

The flammable mixture is ignited by an electrical spark (100 mJ) located at the center of the closed end of the tube (Figure 1). The tests are reproducible as the Figure 4 shows. The estimated errors on overpressure between two tests in the same configuration are around 5 % when the mixture is provided by the same batch and 10 % when the mixture is provided by two different batches.

Analysis of a typical test

This section presents the typical result obtained with an explosion of a stoichiometric methane-air mixture in the transparent 250 mm tube. Figure 5 presents the pressure signals measured in the transparent 250 mm tube near the ignition source, at 5.5 m and at 15.5 m from the ignition source. Figure 6 presents pictures of flame propagation when the transparent tube is captured on film in the middle of pipe. The flame is a paraboloid of revolution during its propagation.

Figure 7 presents the flame trajectory and the flame velocity. The flame trajectory is obtained thanks to a fast camera movie. A first pressure rise is observed up to 120 mbar at 75 ms. This first pressure peak P1 accurately corresponds to the first flame development of the flame from the ignition point to the side walls of the tube [START_REF] Kerampran | Etude des mécanismes d'accélération des flammes se propageant depuis l'extrémité fermée vers l'extrémité ouverte de tubes horizontaux de longueur variable, thèse de docteur de l[END_REF]. This first elongation could be the establishment of Darrieus-Landau instability [START_REF] Markstein | Non-steady flame propagation[END_REF]. Then the flame stops and the pressure decreases. Afterwards, the flame is submitted to strong accelerations, decelerations and stops due to the interactions with the acoustic waves in the tube. Indeed, the typical acoustic frequency of the tube can be estimated with the formula f = n.c/4L where c is the sound velocity estimated around 340 m/s in stoichiometric methane air mixture, L is the length of the tube and n=(3p+1) where p is an integer starting at zero and represents the number of vibration harmonic. The frequency of the second harmonics is around10 Hz and the frequency of velocity oscillation is about 11 Hz. These interactions are also visible when the flame velocity and flow velocity (measured by pitot probes in the middle of a section of tube at 10 and 18 m -Figure 9) are compared showing a good matching of flame velocity and flow velocity oscillations. The flame velocity shows strong oscillations at a frequency increasing as function of the flame position in the pipe (Figure 8). The blue line in the Figure 7 seems to show that the flame progresses at a mean velocity which increases. Thus, the flame seems to have, on average, a uniformly accelerated movement as suggested earlier [START_REF] Proust | Gas flame acceleration in long ducts[END_REF]. Note also that the amplitude of the oscillations of the flame velocity increase during the flame propagation. The value of this average acceleration is about 130 m/s 2 .

Nature of the flow

Specific measures of flow velocity are realized with pitot probes at 10,5 m from the ignition source one on the axis of the pipe and the second near the wall at 0.5 cm from the wall. The principle of the Pitot probes (Figure 3) is to create a break point where the flow line has a null velocity. The dynamic pressure is also null at this point, the total pressure is measured to which the "static" pressure measured by a pressure point downstream of the flow is subtracted from.

The Pitot probes linked with a differential pressure transducer, allows to measure the difference of the dynamic pressure and then to deduce, thanks to the Bernouilli's principle, the velocity of the flow and the fluctuations in the flow. Figure 10 and Figure 11 presents flow velocity measured in the 250 mm plastic and steel pipes and in the 150 mm steel pipe. This analysis strongly suggests that the flow in front of a flame propagating in a pipe is very different from a standard established boundary flow: a flat velocity profile, with very limited boundary, and as a consequence hardly turbulent seems to be produced like a "plug flow".

Parametric investigations

The influence of diameter, roughness of pipes and initial conditions was investigated.

Influence of the pipe diameter

The impact of the tube diameter was experimentally studied in the 24 m long transparent tube with diameters of 150 and 250 mm for stoichiometric methane/air mixture (Figure 14) The common point is the presence of the first peak of overpressure (P1) but the rest of the signal differ significantly with a second significant bump in the smaller tube, much less obvious in the larger pipe. It is known that the first peak corresponds to the first flame elongation, the overpressure in the 150 mm tube (around 130 mbar) is on the same order of magnitude than in the 250 mm tube (around 120 mbar). This value seems to be independent from the pipe diameter, which is consistent with Proust results [START_REF] Proust | Gas flame acceleration in long ducts[END_REF]. From the first peak P1, the flow velocity produced by the flame can be deduced: P1=ρ.c.Uflame [START_REF] Kerampran | Etude des mécanismes d'accélération des flammes se propageant depuis l'extrémité fermée vers l'extrémité ouverte de tubes horizontaux de longueur variable, thèse de docteur de l[END_REF] ρ : density of gases, Uflame : flame velocity, c : sound speed) . A maximum value of 25 m/s is calculated. Note that in the present configuration, this estimation can be compared to the flow velocity measured at 10 m from the ignition. Further, if the burnt product can be considered to be at rest behind the flame front, the flow velocity should be close to the flame velocity. In transparent tubes, the flame velocity can be estimated visually. When all these informations are confronted together a consistent result is obtained (Figure 15) demonstrating the flame and flow dynamics are correctly captured. The maximum flame velocity, 25 m/s, is much large than the expansion velocity E.Slad (E : expansion ratio of burnt gases, Slad : laminar flame speed). The ratio between both parameters is a flame folding factor (flame area to cross section) which amounts 10, consistently with previous findings [START_REF] Daubech | Contribution à l'étude de l'effet de l'hétérogénéité d'un prémélange gazeux sur la propagation d'une flamme dans un tube clos[END_REF].

Influence of pipe material

The influence of pipe material is studied by comparing the pressure signals obtained with the 250 mm plastic transparent tube and with the 250 mm steel tube (Figure 16) still for stoichiometric methane-air mixture. The transparent plastic pipe is considered as a perfectly smooth pipe. The roughness of metallic pipe is due to the oxidation of the wall. It can be estimated that the roughness is 100 μm. The global trend of the pressure signal is conserved. The frequencies of the oscillations are similar but their amplitudes differ. It is twice larger in the plastic pipe (signals at 15 m). Up to know there is no clear explanation to these differences.

Influence if the initial conditions

The impact of the initial conditions was experimentally investigated studied in the 24 m long metallic tube with diameter of 150 mm for stoichiometric methane/air mixture (Figure 17). The flammable mixture is ignited with an electrical spark (energy : 100 mJ) and a pyrotechnical match (energy : 60 J). The spark ignition is punctual whereas the ignition with a pyrotechnical match is more diffuse since several hot particles are emitted

In contradiction with the usual opinion that a more powerful ignition source would lead to a more severe explosion, the contrary is obtained. A possible explanation would be that multiple ignitions, like with the pyrotechnical match, would inhibit the first flame elongation phase from which the subsequent flowfield in the pipe depends on. The flame reaches earlier the edge than in punctual ignition. A further illustration of the potential influence of the first phase of the explosion (first pressure pulse) is provide in Figure 18 showing the pressure signals for stiochiometric CH4/air mixtures burning in 150 mm plastic and steel pipes (Figure 18). The slope and timing of the first pressure pulse is the same in both situations but P1 is higher in the steel pipe. This higher value could be explained by the influence of sensor ports located about 1 m from the ignition source (Figure 2). It is believed that he explosion of the small gas pockets trapped into the sensor ports could destabilize the flame, increasing its area and thus the first peak value.

Discussion

The present results and the parametric investigations confirm the previous findings [START_REF] Proust | Gas flame acceleration in long ducts[END_REF]. If the mechanisms explaining the flame propagation regimes still need to be clarified, it seems nevertheless that the first pressure pulse plays a significant role. The latter is clearly produced by the initial flame ball development and subsequent elongation. The subsequent flame propagation, after this initial elongation phase, is a constantly accelerating flame resembling a piston in terms of pressure generation. Over this general trend pressure oscillation are superposed due the acoustic resonance of the pipe.

The first elongation last from the spark ignition time to when the flame reaches the pipe lateral walls. As suggested earlier [START_REF] Proust | Gas flame acceleration in long ducts[END_REF][START_REF] Kerampran | Etude des mécanismes d'accélération des flammes se propageant depuis l'extrémité fermée vers l'extrémité ouverte de tubes horizontaux de longueur variable, thèse de docteur de l[END_REF], and confirmed by the present measurements, the magnitude of the overpressure is proportional to the material velocity of the flow, the latter parameter being very close to the flame velocity. The maximum flame velocity is proportional to the expansion velocity of the burnt gases represented by the parameter E.Slad. The coefficient of proportionality is simply the maximum elongation of the flame Af/A (or "aspect ratio"). The maximum flame elongation is about 10, which is close to the value measured on videos by Kerampran.

After the first peak, a self-accelerated motion of flame occurs despite the oscillation. A reasonable correlation between the expansion velocity of the burnt gases represented by the parameter E.Slad and some value of flame acceleration can be found [START_REF] Proust | Gas flame acceleration in long ducts[END_REF]. This analysis of flow velocity in explosion tests in plastic and metallic pipes shows that the flow created in front of the flame seems to be not turbulent and quite flat. The flow in front of the flame may be considered as a pressure zone (pressure plug) whose length depends on the first flame development and which propagates in the pipe.

The self-accelerated motion can be estimated using Vf = af.t where Vf is the flame velocity, af is the flame acceleration and t is the time. For stiochiometric methane/air mixture, this value is about 130 m/s2.

The last point deals with the acoustic oscillations after the first overpressure peak. It's possible to estimate the resonance frequencies of constant section pipe although it contains gases of different nature [START_REF] Kerampran | Etude des mécanismes d'accélération des flammes se propageant depuis l'extrémité fermée vers l'extrémité ouverte de tubes horizontaux de longueur variable, thèse de docteur de l[END_REF]. The tube with length L, closed at one end and open at the other, contains burnt gases up to the xf position. In this region, the acoustic waves propagate at ab velocity. In the other part of the tube, it contains unburned gases where acoustic waves propagates at au. The period T and the frequencies f of third harmonics oscillations are given by:

Figure 1 :

 1 Figure 1 : Overall device

Figure 3 :

 3 Figure 3 : Ionization probe (a), Pitot probe (b) and ignition source (c)

Figure 4 :

 4 Figure 4: Overpressure signals for a test in 150 mm steel pipe for 10 % CH4-air mixtures provided by the same batch (a) and for a test in 250 mm plastic pipe for 10 % CH4-air mixtures provided by two different batch (b)

Figure 5 :Figure 6 :Figure 7 :

 567 Figure 5 : Pressure signals measured in the transparent 250 mm tube near the ignition source, at 5.5 m and at 15.5 m from ignition source-stoichiometric methane air flame at ambient conditions

Figure 8 :

 8 Figure 8 : Frequency of oscillation versus distance from ignition source -stoichiometric methane air flame at ambient conditions (250 mm tube)

Figure 10 :

 10 Figure 10 : Velocity in plastic and steel pipes in the middle of a section and near the edge at 10 m from ignition source -stoichiometric methane air mixture at ambient conditions (24 m long)

Figure 12 :

 12 Figure 12 : Flow velocity measured in a 300 mm steel pipes (20 m long) in which the flow is established using a fan

Figure 14 :

 14 Figure 14: Impact of the tube diameter in the 24 m long transparent tube with diameters of 150 and 250 mm for stoichiometric methane/air mixture

Figure 15 :

 15 Figure 15 : Comparison between the flame velocity measured in transparent tube for which the first flame development is similar, the flow velocity measured by pitot probe at 10 m and the flame velocity deduced from the overpressure-stoichiometric methane air mixture at ambient conditions (250 mm pipe -24 m long)

Figure 16 :

 16 Figure 16 : Comparison between the pressure signals in 250 mm plastic transparent tube and 250 mm steel tube near the ignition, at 10 m and 15 m to the ignition source-stoichiometric methane/air mixture.

Figure 17 :

 17 Figure 17 : Pressure signals near the ignition source for stoichiometric methane/air mixture ignited with spark and pyrotechnical match -24 m long 150 mm steel pipe

Figure 18 :

 18 Figure 18 : Comparison between the pressure signals in 150 mm plastic transparent tube and 150 mm steel tube near the ignition, at 10 m and 15 m to the ignition source -10% methane/air mixture

  Figure19presents a comparison between the estimated frequencies and the frequencies deduced from the flame velocity.

Figure 19 :

 19 Figure 19 : Comparison between the estimated frequencies and the frequencies deduced from the flame velocity

Figure 20

 20 Figure 20 presents a comparison with experimental and estimated flame velocity. The first phase of flame development is completely characterized. The second phase is the coupling of the self-accelerated motion of flame whose value of acceleration is defined by E.Slad and the acoustic oscillations of pipe. The value of acceleration for stoichiometric methane/air mixture is 130 m/s 2 .

Figure 20 :

 20 Figure 20 : Comparison with experimental and estimated flame velocityThe estimation of overpressure thanks to the flame velocity is presented in Figure21. The model gives a reasonable agreement for the first acoustic oscillations, but overestimated the value of pressure. The mechanism of pressure production must be different with acoustic oscillations.

Figure 21 :

 21 Figure 21 : Estimation of overpressure thanks to the flame velocity

Acknowledgements

The authors gratefully acknowledge the financial contribution from the Research Fund for Coal and Steel (RFCS) for the EXPRO project.