# An Individual-Based Model of Zebrafish Population Dynamics Accounting for Energy Dynamics

Rémy Beaudouin, Benoit Goussen, Benjamin Piccini, Starrlight Augustine,

James Devillers, François Brion, Alexandre R.R. Pery

► To cite this version:

Rémy Beaudouin, Benoit Goussen, Benjamin Piccini, Starrlight Augustine, James Devillers, et al. An Individual-Based Model of Zebrafish Population Dynamics Accounting for Energy Dynamics. PLoS ONE, 2015, 10 (5), pp.e0125841. 10.1371/journal.pone.0125841. hal-01690258

HAL Id: hal-01690258

https://hal.science/hal-01690258

Submitted on 22 Jan 2018

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 Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

RESEARCH ARTICLE

# An Individual-Based Model of Zebrafish Population Dynamics Accounting for Energy Dynamics

Rémy Beaudouin¹*, Benoit Goussen¹, Benjamin Piccini², Starrlight Augustine³, James Devillers⁴, François Brion², Alexandre R. R. Péry¹,⁵

1. Unité Modèles pour l'Ecotoxicologie et la Toxicologie (METO), Institut National de l'Environnement Industriel et des Risques (INERIS), Verneuil en Halatte, France,
2. Unité Ecotoxicologie in vitro et in vivo (ECOT), Institut National de l'Environnement Industriel et des Risques (INERIS), Verneuil en Halatte, France,
3. Centre for Ocean Life, National Institute of Aquatic Resources, Technical University of Denmark, Charlottenlund, Denmark,
4. CTIS, Rillieux La Pape, France,
5. AgroParisTech, Paris, France

* remy.beaudouin@ineris.fr

## Abstract

Developing population dynamics models for zebrafish is crucial in order to extrapolate from toxicity data measured at the organism level to biological levels relevant to support and enhance ecological risk assessment. To achieve this, a dynamic energy budget for individual zebrafish (DEB model) was coupled to an individual based model of zebrafish population dynamics (IBM model). Next, we fitted the DEB model to new experimental data on zebrafish growth and reproduction thus improving existing models. We further analysed the DEB-model and DEB-IBM using a sensitivity analysis. Finally, the predictions of the DEB-IBM were compared to existing observations on natural zebrafish populations and the predicted population dynamics are realistic. While our zebrafish DEB-IBM model can still be improved by acquiring new experimental data on the most uncertain processes (e.g. survival or feeding), it can already serve to predict the impact of compounds at the population level.

## Introduction

Data used to estimate the likelihood of adverse ecological effects typically include responses of survival, growth, or reproduction of individuals measured after a specific exposure duration under constant laboratory conditions and in absence of ecological stress (e.g. predation and competition) [1]. These organism-level endpoints are far from the ecological features that the process aims to protect. Indeed, ecological risk assessment should protect the long-term persistence of populations of species in space and time under naturally varying field conditions and in the presence of other stressors (e.g. food limitation). However, except the ecotoxicological data provided by mesocosm experiments and a few field studies [2–5], data on impacts of chemical substances on populations or higher biological levels are very sparse. In this context, population models can play an important role in bridging the gap between what is measured (organism-level endpoints) and what needs to be protected (population-level endpoints) [6].

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

OPEN ACCESS

**Citation:** Beaudouin R, Goussen B, Piccini B, Augustine S, Devillers J, Brion F, et al. (2015) An Individual-Based Model of Zebrafish Population Dynamics Accounting for Energy Dynamics. PLoS ONE 10(5): e0125841. doi:10.1371/journal.pone.0125841

**Academic Editor:** Christopher Joseph Salice, Towson University, UNITED STATES

**Received:** October 24, 2014

**Accepted:** March 26, 2015

**Published:** May 4, 2015

**Copyright:** © 2015 Beaudouin et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the paper and its Supporting Information files.

**Funding:** This study was funded by the French Ministry in charge of Ecology and Sustainable Development, within the frameworks of the National Research Programme on Endocrine Disruption (PNRPE) and of the Programme 190.

**Competing Interests:** The authors have declared that no competing interests exist.
Zebrafish (*Danio rerio*) offers many practical and methodological advantages that make this small fish species an attractive vertebrate model now used worldwide in a variety of biological disciplines ranging from basic developmental biology to applied toxicology [7]. In the last few years, the popularity of the zebrafish model for investigating chemical ecotoxicity increased notably to address the issues posed by Endocrine Disrupting Chemicals (EDCs) and their reproductive and developmental effects in vertebrate species. This lead to significant knowledge and methodological advances to assess modes of action and effects of EDCs [8–11].

While these experimental approaches can provide relevant information regarding the capacity of chemicals to disrupt key endocrine-regulated physiological processes at organism level, they are limited to address impact at population level. Hence, the development of population dynamics model for zebrafish appears crucial in order to extrapolate from toxicity data measured on organism to biological levels relevant to support and enhance ecological risk assessment [6]. Individual-based models (IBMs) are suitable population models to increase the relevance of ecotoxicity tests since they incorporate available mechanistic knowledge on the links between responses at the individual level and responses at the population level [12]. Predictions for stressed populations could be achieved provided that dose response relationships are known for the key parameters [13–15]. Moreover, toxicokinetic/toxicodynamic (TK/TD) models, able to link exposure concentration and effects at target organ/system level, can be easily integrated into IBMs and then, the processes and the corresponding effects that lead to toxicity within an organism can be dynamically simulated [16]. A physiologically-based toxicokinetic model [7] and IBM [17] have been recently developed for the zebrafish. However, the IBM developed by Hazlerigg et al. [17] describes life cycle processes without integrated mathematical framework and without integrating some of the main ecological factors driving the population dynamics: photoperiod, water temperature, and food availability [18, 19]. Hence, a more mechanistical approach would offer advantages to improve the realism of such an IBM and to allow, in the future, the integration of toxic effects on physiological processes.

The consequence of the high level of complexity and precision of IBMs is that these models need a large amount of data on the physiological processes involved at the organism-level [20]. The Dynamic Energy-Budget (DEB) theory [21] is a unified mathematical framework composed of models describing mechanistically the acquisition and use of energy to relate growth, maintenance, development, and reproduction to environmental parameters such as food availability and temperature. Then, nesting a DEB model within a population dynamics IBM can be useful when existing models and/or data on organism level physiology are available [2, 22]. A zebrafish DEB model has been developed by Augustine et al. [23] and provides part of the physiological information essential to develop an IBM.

This study aimed to develop a model of the zebrafish population dynamics by coupling the DEB model with an IBM, hereafter referred to as the DEB-IBM. To achieve this, new experimental data on zebrafish growth and reproduction were produced and alternative hypotheses to the existing zebrafish DEB model [23] were used. Then a DEB-IBM integrating the main ecological factors (photoperiod, temperature, and food dynamics) was developed. To perform predictions for stressed populations, the first step was to carefully check the prediction accuracy and reliability of the model under control conditions (e.g. fish population dynamic in the absence of exposure [2]). To this end, we performed a thorough sensitivity analysis of our DEB-IBM and compared the model predictions of the population dynamic in the absence of exposure to field data collected by Spence et al. [24] and Hazlerigg et al. [17].
# Materials and Methods

## New experimental data production

First, length-mass allometry of juveniles was investigated to test the anisomorphic growth hypothesis (change in shape during growth *i.e.*, surface area grows proportional to volume) proposed by Augustine et al. [23]. Other experiments were designed to produce data to adjust the DEB model for the juveniles (lengths monitored at different feeding regimes), females (fecundity monitored at two different temperatures) and males (female and male lengths monitored at different feeding regimes). All experiments were approved by a local ethics committee: Institut National de l'Environnement Industriel et des Risques (INERIS) "ECOT_11_050_RBU" and "ECOT_11_053_RBu". Fish were killed by an overdose of MS222.

**Zebrafish origin, maintenance and culture.** Wild type larvae and adult zebrafish (AB strain) originated from our breeding unit (INERIS, Verneuil-en-Halatte, France). Adult zebrafish were maintained in 3.5 L aquaria in a recirculation system (Zebtec, Techniplast, France) at a 14:10 light: dark photoperiod, 27.1 ± 2.0°C. They were allowed to reproduce (2 males for 1 female) in 50L tanks (14:10 light: dark 29°C). Fertilized eggs were collected, disinfected 5 min in water supplemented with 0.1% of commercial bleach (2.6% of sodium hypochlorite). Zebrafish were maintained under semi-static conditions until 21 days and then transferred into 3.5 L aquaria of the recirculation system. For zebrafish culture, fish were fed with commercial diets of different qualities and size to cover all stages of development of fish (Fig. A in S1 Text).

**Juvenile allometry.** To assess length-mass allometry of juveniles, standard length (measured from the tip of the nose to the base of the caudal fin; mm) and wet mass (mg) of 60 juveniles aged from 10 to 50 days post fecundation (dpf) were measured with a precision of 0.1 mm and 0.1 mg, respectively. The relationship between the logarithms of individual length and body masses was then assessed using a linear model.

**Growth data.** Juveniles' lengths were monitored from 49 to 107 dpf. The juveniles were randomly assigned to three feeding treatments differing in feeding frequencies. Each feeding treatment was replicated three times with $n = 15$ fish per replicate (15 fish × 3 replicates × 3 treatments = 135 fish in 9 groups/tanks). Fish were fed twice a day with dehydrated food at 6% of their median mass according to the following feeding treatments: (i) each day (Food 1), (ii) two out of three days (Food 0.75), and (iii) every two days (Food 0.5). Fish mass was derived from the length by length-mass allometry relationship and updated at each length measurement. The fish were measured at 49, 56, 63, 70, 77, 91, 104, and 107 dpf.

The measurements of the standard length (mm) were measured by placing the fish onto a Petri dish containing water (less than 5 mm deep). The Petri dish was placed on graph paper and photographed with a digital camera. The digital file was then used to measure the length of the fish to the nearest 0.1 mm using Image J software [25].

Another experiment (following the same protocol) was held with fish groups composed of 30 individuals per group and only two feeding treatments (three replicates per feeding treatment, six fish groups). These individuals were fed following two feeding treatments (6% of their median mass): (i) each day (Food 1), (ii) every two days (Food 0.5). At the end of the experiment, these individuals were sexed by observing their gonad under a binocular microscope.

**Reproduction data.** Two groups of eight females of similar length (from 35 to 43 mm in standard length) were randomly selected. The first group was maintained at 29°C and reproductive outputs were monitored daily for 21 days. The second group was maintained at 26°C and reproductive outputs were monitored daily for 11 days. Each female was held with two males (selected randomly among male of similar length) in a 5 L aquaria supplied continuously with dechlorinated tap water at a flow rate of 5 L·h⁻¹. Each evening after the last feeding event, a spawning trap (glass square-shaped beakers with a 2 mm mesh net on the top) was placed in
each aquarium. In the morning, the spawning traps were collected and the eggs were counted. Thirty eggs were randomly selected, cleaned, and incubated until 24 hours post-fertilization (1 dpf) at 28°C in autoclaved and aerated aquarium water. After the incubation period, unfertilized eggs were counted. For each group, the spawning at the day 0 was not counted, because this spawning (>500 eggs) corresponded to the eggs cumulated for several days. During the experiment, fish were fed four times a day with SDS 400 (SDS diet) and live brine shrimp *Artemia salina* at a rate of 30% of the fish biomass per day.

## Dynamic Energy Budget model

DEB theory [21] is based on a mathematical description of the uptake and use of energy within an organism, to describe mechanistically the energy flux to physiological process. Moreover, in toxicology, the analysis of (eco)toxicological data through models based on DEB theory (i.e. DEBTox models) is relevant to mechanistically assess the effect of toxic compounds [26, 27]. According to the DEB theory, energy is taken up from food, assimilated, and stored into reserves. This energy is then dispatched between three main processes: (i) maintenance, (ii) growth, and (iii) maturation/reproduction. The general DEB framework assumes that individuals' growth follows a von Bertalanffy growth curve provided feeding is ad libitum or at constant density [21]. However, as presented by _bookmark2Augustine et al. [23] the zebrafish growth curve could be a sigmoid. The authors used several dataset from Eaton and Farley [28], Bagatoo et al. [29], Schilling [30], Lawrence et al. [31], Best et al. [32], and Gòmez-Requeni et al. [33] to assess their growth curve hypothesis and to estimate parameters of their DEB model. Augustine et al. [23] put forth the hypothesis that the metabolism of larvae accelerates after birth until juvenile stage, i.e. fish of about 10 mm total body length and with adult fins and pigments [34]. The hypothesis of metabolic acceleration comes with concomitant anisomorphic growth [21]. The metabolic acceleration concept is assumed to impact both the growth curve and the incubation time. Here we propose alternative hypotheses based on feeding limitation and environmental factors. Indeed, as demonstrated on other organisms, a sigmoid growth curve can be explained by a limitation either in the uptake [35, 36] or in the quality of the food [36]. We thus adapted the set of equations from Augustine et al. [23] in order to take into account this new hypothesis. The biological interpretation (and values) of the DEB parameters are listed in Table 1 using standard DEB notations [21]. The compound parameters, which are function of primary parameters, are presented in the supporting information (S2 Text). The primary parameters used in the DEB model were corrected depending on the temperature using a temperature correction function based on the Arrhenius equation (detailed in the supporting information). Thereby, parameters $v$ (mm d$^{-1}$) the energy conductance, $\{\dot{p}_{Am}\}$ (J d$^{-1}$ mm$^{-2}$) the maximum surface area specific assimilation rate, and $[\dot{p}_{M}]$ (J d$^{-1}$ mm$^{-3}$) the volume specific somatic maintenance costs were corrected. The DEB model including food limitation and fixed size at puberty reads:

$$S_f(l) = \alpha \left[ 1 - \left( 1 + \frac{l^3}{l^3} \right)^{-1} \right] \quad (1)$$

$$\frac{de}{dt} = \frac{\dot{k}_{M}g}{l} [(1 - s_f)f - e] \quad (2)$$

$$\frac{dl}{dt} = r_B(e - l) \quad \text{with} \quad r_B = \frac{\dot{k}_{M}g}{3(e + g)} \quad (3)$$
<table><caption>Table 1. Parameter abbreviations, values, descriptions and units of the DEB model.</caption><thead><tr><th>Abbreviation</th><th>Value</th><th>Description</th><th>Unit</th><th>Source</th></tr></thead><tbody><tr><td>T<sub>A</sub></td><td>3000</td><td>Arrhenius temperature</td><td>K</td><td>[23]</td></tr><tr><td>T<sub>R</sub></td><td>293</td><td>Reference temperature</td><td>K</td><td>[23]</td></tr><tr><td>δ</td><td>0.20</td><td>Shape coefficient ; V<sup>1/3</sup> = δ.L</td><td>-</td><td>Fitted</td></tr><tr><td>{ṗ<sub>Am</sub>}</td><td>4.72</td><td>Maximum area specific assimilation rate</td><td>J.d<sup>-1</sup>.mm<sup>-2</sup></td><td>Fitted</td></tr><tr><td>v̇</td><td>0.60</td><td>Energy conductance</td><td>mm.d<sup>-1</sup></td><td>Fitted</td></tr><tr><td>κ</td><td>0.70</td><td>Fraction of energy to growth/somatic</td><td>-</td><td>Fitted</td></tr><tr><td>α</td><td>0.84</td><td>Fraction of food accessible at first feeding</td><td>-</td><td>Fitted</td></tr><tr><td>l<sub>b</sub></td><td>0.079</td><td>Scaled length at the first feeding</td><td>-</td><td>Fitted</td></tr><tr><td>l<sub>p</sub></td><td>0.58</td><td>Scaled length at puberty</td><td>-</td><td>Fixed</td></tr><tr><td>l<sub>f</sub></td><td>0.163</td><td>Length at half maximal food assimilation</td><td>-</td><td>Fitted</td></tr><tr><td>[ṗ<sub>M</sub>]</td><td>0.44</td><td>Volume somatic maintenance costs</td><td>J.d<sup>-1</sup>.mm<sup>-3</sup></td><td>Fitted</td></tr><tr><td>[E<sub>G</sub>]</td><td>2.35</td><td>Cost of synthesis of a unit of structure</td><td>J.mm<sup>-3</sup></td><td>Fitted</td></tr><tr><td>E<sub>0</sub></td><td>1.25</td><td>Initial amount of energy in an egg</td><td>J</td><td>Fitted</td></tr><tr><td>L<sub>0</sub></td><td>0.25</td><td>Initial length</td><td>mm</td><td>Fixed</td></tr><tr><td>R<sub>M</sub></td><td>406</td><td>Maximum number of egg per day</td><td>d<sup>-1</sup></td><td>Fitted</td></tr><tr><td>f<sub>lim</sub></td><td>0.92</td><td>Male appetite factor limitation</td><td>-</td><td>Fitted</td></tr></tbody></table>

doi:10.1371/journal.pone.0125841.t001

$$ \frac{dR}{dt} = \frac{R_M}{1 - l_p^3} \left( \frac{g + l}{g + e} e^{l^2} - l_p^3 \right) \quad (4) $$

If $l < l_b$ then $f = 0$ and if $l < l_p$ then $\frac{dR}{dt} = 0$ $sf$ is the size-dependant energy ingestion function, $f$ (-) is the ratio of the actual ingestion rate divided by the maximal ingestion rate for a given body size, $e$ (-) is the scaled reserve density, and $r_B$ (d⁻¹) is the von Bertalanffy growth rate. It should be noticed that the birth corresponds to the opening of the mouth. Thus, $L_b$ (mm) corresponds to the length at the opening of the mouth. $L_f$ (mm) is the length at which the ingestion rate is half the maximum ingestion rate and $L_p$ (mm) is the length at puberty. All the physical lengths (i.e. $L$, $L_b$, $L_f$ and $L_p$) were scaled by the maximal physical length $L_{inf}$ resulting in the scaled lengths noted $l$, $l_b$, $l_p$ and $l_f$, respectively. Maximal physical length ($L_{inf}$) is a function of the DEB maximal length and a shape coefficient (Table A in S3 Text). $R$ represents the cumulative number of eggs produced, $R_M$ (d⁻¹) the maximum reproduction rate, and $g$ (-) the energy investment ratio. The DEB model was also adapted for males by assuming that after puberty their food intake is modified by an appetite factor ($f_{lim}$-). This factor was calibrated on the experimental data (S1 Text). Thereby, after puberty, the parameter $f$ was modified using the following equation:

$$ f_{male} = f_{lim} \times f \quad (5) $$

The DEB model was calibrated using data from our experiments (Fig. B and C in S1 Text) in addition to a subset of the data used by Augustine et al. [23]. More information on the calibration process is provided in supporting information (S2 Text).

## Individual-based model of zebrafish population dynamics

The individual based-model (IBM) description follows the ODD (Overview, Design concepts, Details) protocol [37]. Our IBM was developed from the DEB model presented in the previous section. Mating behaviour and survival sub-model was adapted from the model proposed by Hazlerigg et al. [17]. The detailed description of the Design concepts, initialization, Input data
and submodel are provided in supporting information (S3 Text). The model was implemented under Netlogo 5.0.4 [38] and the code of the model is provided in supporting information (S4 Text, S1 File and S2 File).

**Purpose.** The zebrafish IBM was developed to extrapolate from individual data measured on organisms under laboratory conditions to biological endpoints convenient to support ecological risk assessment of toxic compounds. To this end, the zebrafish model was developed for a population living in well-suited habitat, i.e. the abiotic factors (e.g. temperature, photoperiod) were modelled to represent the natural habitat of this species which are found in North-Eastern India, Bangladesh and Nepal. Field populations were observed near Kolkata [39], thus we selected as reference climatology that of Kolkata, i.e. a monsoon tropical climate. The zebrafish population that are modelled here inhabits a « median » habitat of the ones described in [24, 39], i.e. an isolated well-vegetated pond.

**Entities, state variables, and scales.** The model includes two categories of agents: patches and fish agents. Three different patches were modelled: vegetation cover, breeding grounds, and open water patches. Fish agents were sub-divided into four sub-stages: eggs-larvae, juveniles, females, and males. Eggs and larvae (before initiation of feeding) have very similar attributes and processes, and so were therefore grouped into one stage. Eggs-larvae change into juveniles at the onset of exogenous feeding and become females or males at sexual maturity. All agents are characterised by the state variables age (dpf), generation (0 for the initial fish), and survival rate (day⁻¹). An eggs-larvae agent represents an entire brood and is characterised by one state variable: the number of surviving eggs. The juvenile, female, and male agents have additional state variables: length (mm), wet body weight (mg), sex, scaled energy density (dimensionless), degree of satiety (dimensionless), and a variable accounting for the variability of the performance between individuals (dimensionless). Female and male agents have each another state variable: the number of eggs (dimensionless) and a state variable characterising the absence/presence and position of the mating territory, respectively.

The fish agents were spatialized (i.e. agent located in geometrical space) in a two dimensional grid of 900 patches (30 × 30 patches of 20 × 20 cm) representing a typical pond of 36 m². The pond water depth was 0.5 m. Each patch (20 cm × 20 cm × 50 cm) represents the real size of observed zebrafish territories [40]. The time step was 1 day for the IBM (1/10 day for DEB model sub-part).

**Process overview and scheduling.** The schedule is schematized in Fig 1. At each time step, first the environment global variables, i.e. water temperature, photoperiod, population variables (e.g. fish biomass), and the food dynamic model are updated. Second, fish move each in turn and in random order: juveniles move to find a vegetated patch, males try to establish a territory or move randomly, and females try to find a free males with a territory or move randomly. Then, the fish agents challenge their survival. If alive, fish agents grow. After that, if the conditions are met, eggs-larvae agents can hatch, juvenile agents can undergo puberty, and female agents can spawn. At the end of the time step, the simulation results are backed up.

## Statistical methods

**Sensitivity analysis.** The sensitivity analysis orders the inputs by importance, identifying the main contributors to the variation in the model outcomes. Two different methods were used for the DEB model and IBM analysis to keep a reasonable time of computing. Sensitivity analysis of the DEB model was performed using the variance-based Sobol method [41, 42]. This sensitivity analysis method is a global and model independent method (nonlinear and non-monotonic functions and models can be handled) and is based on variance decomposition (quantifying the amount of variance that each parameter contributes to the unconditional
<figure><img src="image_4.png" /><figcaption>**Fig 1. Schema of the schedule of the model of zebrafish population dynamics.** Schema of food sub-model was issued from Li and Yakupitiyage [61] and schema of the dynamic energy budget (DEB) model was issued from Kooijman [21]. An arrow before an action denotes that action was realized in random order by fish agents one after another in swift succession.</figcaption></figure>

doi:10.1371/journal.pone.0125841.g001

variance of the model output). This approach allows identification of the contribution of each parameter (First order Sobol' sensitivity index) and the interactions between parameters (Total Sobol' sensitivity index). Sensitivity analysis of the IBM was performed using Morris method [43]. This method identifies the few important factors at a limited cost of simulations. It also gives information on the global contribution of each parameter and the linearity of the effects or the interaction with other parameters. Details of these analyses are provided in (S5 Text).

**Model parameterisation.** DEB model parameters calibrations were performed using GNU MCSim software [44], which is a program for the statistical analysis of Bayesian hierarchical models by Markov Chain Monte Carlo (MCMC). Details are provided in the supporting information (S2 Text). The parameter values of the IBM were picked-up from the literature (Tables 2 and 3) and three parameters were fitted (no *a priori* information): πₐ, πc, and Fᵢₙ (survival parameters and food input during the monsoon; Table 3). Two different datasets were used to calibrate and assess predictability of the DEB-IBM model: data collected by Hazlerigg et al. [17] and data collected by Spence et al. [24], respectively. These authors have measured length on a random sample of fish captured from a sub-part of the monitored populations (field sampling methods were ineffective in catching smaller individuals). These data were compared with the fish length distribution predicted by the model for the sub-part of the population considered (after the simulation of three years to stabilise the population; 1000 simulations). Hence, we assume that the probability to be sampled for a given fish was related to the frequency of the fish length in the sub-population. Calibration was done using a genetic algorithm provided by the open-source software BehaviorSearch [45] with a generational population-model and a mutation-rate, crossover-rate, population-size, and tournament-size equal to 0.01, 0.7, 50 and 3, respectively. The distance was the sum of squares of the difference between the length distribution predicted by the model (after 1110 days, fish with a length > 18.4 mm) and the length distributions observed by Hazlerigg et al. [17] (length distributions had similar bin of 1 mm).

**Statistical analysis.** Data were analysed using ANOVA or RM-ANOVA. When not normal or not homoscedastic, data were transformed using the Box Cox method [46]. Scheffe's
<table><caption>Table 2. Parameter abbreviations, values, descriptions and units of the food sub-model.</caption><thead><tr><th>Abbreviation</th><th>Value</th><th>Description</th><th>Unit</th><th>Source</th></tr></thead><tbody><tr><td>G2Kcal</td><td>1.88</td><td>Coefficient of energy by mass of zebrafish</td><td>Kcal.g<sup>-1</sup></td><td>[23]</td></tr><tr><td>s</td><td>21.08</td><td>Proportionality of food nutrient quantity to fish biomass</td><td>-</td><td>[61]</td></tr><tr><td>h<sub>n</sub></td><td>0.2</td><td>Half saturation nitrogen (N)</td><td>mg.L<sup>-1</sup></td><td>[61]</td></tr><tr><td>h<sub>p</sub></td><td>0.02</td><td>Half saturation phosphorus (P)</td><td>mg.L<sup>-1</sup></td><td>[61]</td></tr><tr><td>K<sub>n</sub></td><td>0.01</td><td>N-fixation coefficient of phytoplankton</td><td>g.kcal<sup>-1</sup>.day<sup>-1</sup></td><td>[61]</td></tr><tr><td>K<sub>fn</sub></td><td>0.017</td><td>N content of fish</td><td>g.kcal<sup>-1</sup></td><td>[61]</td></tr><tr><td>K<sub>an</sub></td><td>0.0224</td><td>N content of phytoplankton</td><td>g.kcal<sup>-1</sup></td><td>[61]</td></tr><tr><td>K<sub>hn</sub></td><td>0.0192</td><td>N content of heterotrophic components</td><td>g.kcal<sup>-1</sup></td><td>[61]</td></tr><tr><td>K<sub>sn</sub></td><td>0.003</td><td>Release coefficient of N in sediment</td><td>day<sup>-1</sup></td><td>[61]</td></tr><tr><td>K<sub>nl</sub></td><td>0.17</td><td>coefficient of inorganic N loss to air</td><td>day<sup>-1</sup></td><td>[61]</td></tr><tr><td>K<sub>ap</sub></td><td>0.001</td><td>P content of phytoplankton</td><td>g.kcal<sup>-1</sup></td><td>[61]</td></tr><tr><td>K<sub>hp</sub></td><td>0.001</td><td>P content of heterotrophic components</td><td>g.kcal<sup>-1</sup></td><td>[61]</td></tr><tr><td>K<sub>pr</sub></td><td>0.0006</td><td>Release coefficient of P in sediment</td><td>m.day<sup>-1</sup></td><td>[61]</td></tr><tr><td>K<sub>ps</sub></td><td>0.28</td><td>Coefficient of inorganic P sedimentation to sediment</td><td>m.day<sup>-1</sup></td><td>[61]</td></tr><tr><td>K<sub>s</sub></td><td>0.14</td><td>Coefficient of heterotrophic food sedimentation</td><td>day<sup>-1</sup></td><td>[61]</td></tr><tr><td>K<sub>d</sub></td><td>0.12</td><td>Coefficient of heterotrophic food decomposition</td><td>day<sup>-1</sup></td><td>[61]</td></tr><tr><td>λ<sub>max</sub></td><td>1.6</td><td>Maximum growth coefficient for phytoplankton growth</td><td>day<sup>-1</sup></td><td>[61]</td></tr><tr><td>K<sub>r</sub></td><td>0.1</td><td>Coefficient of phytoplankton respiration</td><td>day<sup>-1</sup></td><td>[61]</td></tr><tr><td>K<sub>ml</sub></td><td>0.6</td><td>Coefficient of autotrophic food entering heterotrophic food pool</td><td>day<sup>-1</sup></td><td>[61]</td></tr><tr><td>Ir</td><td>6.547</td><td>Reference solar radiation</td><td>10<sup>6</sup> cal.m<sup>-2</sup>.day<sup>-1</sup></td><td>[61]</td></tr><tr><td>Ir<sub>a</sub></td><td>0.000017</td><td>Light extinct coefficients for autotrophic components</td><td>pond.Kcal<sup>-1</sup></td><td>[61]</td></tr><tr><td>Ir<sub>b</sub></td><td>0.000015</td><td>Light extinct coefficients for heterotrophic components</td><td>pond.Kcal<sup>-1</sup></td><td>[61]</td></tr><tr><td>T<sub>opta</sub></td><td>30</td><td>Optimal temperature for phytoplankton growth</td><td>°C</td><td>[61]</td></tr><tr><td>K<sub>T1</sub></td><td>0.004</td><td>Effects of temperature below Topta on growth</td><td>°C<sup>-2</sup></td><td>[61]</td></tr><tr><td>K<sub>T2</sub></td><td>0.008</td><td>Effects of temperature above Topta on growth</td><td>°C<sup>-2</sup></td><td>[61]</td></tr></tbody></table>

doi:10.1371/journal.pone.0125841.t002

post hoc test was used to compare two groups with a significant limit of p < 0.05. Statistical analyses were performed using R statistical environment [47].

# Results

## New experimental data

All the newly generated data on juvenile allometry, growth, and reproduction are presented in S1 Text. A linear model was adjusted between the mass logarithm and the length logarithm. The slope and intercept of the linear model between the mass logarithm and the length logarithm were 3.22 with a 95% confidence interval of (2.97; 3.46) and -5.24 with a 95% confidence interval of (-5.83; -4.65), respectively. Hence, the slope of the linear model did not significantly differ from 3 (organism mass is proportional to the cubic length in case of isomorphic allometry; Fig. F in S1 Text). Significant effects of feeding regime and time were observed in the group of 15 and 30 fish. For the second group (30 fish), a significant sexual dimorphism on the length was observed at 107 dpf.

Inter-spawn intervals were 1.59 ± 1.10 days and 1.23 ± 0.76 days for the females held at 26°C and 29°C, respectively (Fig. D in S1 Text). The interval was equal to 24h for 71% and 85% of the intervals of the females monitored at 26°C and 29°C, respectively. The temperature did not impact the mean number of eggs per spawning (263 ± 104.2 eggs). A significant effect of female/male couple on the percentage of viable eggs at 24h post fecundation was observed. The
<table><caption>Table 3. Parameter abbreviations, values, descriptions and units of the IBM.</caption><thead><tr><td>Abbreviation</td><td>Value</td><td>Description</td><td>Unit</td><td>Reference</td></tr></thead><tbody><tr><td>W<sub>d</sub></td><td>0.5</td><td>Water depth of the pond</td><td>m</td><td>[24, 39]</td></tr><tr><td>W<sub>v</sub></td><td>18</td><td>Water volume of the pond</td><td>m<sup>3</sup></td><td>[24, 39]</td></tr><tr><td>Np<sub>b</sub></td><td>207</td><td>Number of patches with breeding grounds</td><td>-</td><td>[17]</td></tr><tr><td>Np<sub>v</sub></td><td>207</td><td>Number of patches with vegetation cover</td><td>-</td><td>[17]</td></tr><tr><td>Nb<sub>j</sub></td><td>300</td><td>Number of juveniles in the initial population</td><td>-</td><td>[17]</td></tr><tr><td>Nb<sub>m</sub></td><td>35</td><td>Number of males in the initial population</td><td>-</td><td>[17]</td></tr><tr><td>Nb<sub>f</sub></td><td>35</td><td>Number of females in the initial population</td><td>-</td><td>[17]</td></tr><tr><td>End</td><td>1110</td><td>Simulation duration</td><td>d</td><td>This study</td></tr><tr><td>σ</td><td>0.235</td><td>Inter-individual variability of energy acquisition</td><td>-</td><td>This study</td></tr><tr><td>P.τ</td><td>12</td><td>Photoperiod threshold to the onset/offset of the reproduction</td><td>h</td><td>This study</td></tr><tr><td>T.τ</td><td>22.5</td><td>Temperature threshold to the onset/offset of the reproduction</td><td>°C</td><td>[62]</td></tr><tr><td>H<sub>0.5</sub></td><td>24</td><td>Female density inducing 50% reduction of hte hacthing rate</td><td>g.mm<sup>-3</sup></td><td>[63]</td></tr><tr><td>H<sub>max</sub></td><td>0.89</td><td>Optimal hatching rate</td><td>-</td><td>[64–66]</td></tr><tr><td>R.τ</td><td>263</td><td>Number limit of eggs to spawn</td><td>-</td><td>This study</td></tr><tr><td>SR<sub>μ</sub></td><td>50</td><td>Sex-ratio genetic mean</td><td>-</td><td>[67]</td></tr><tr><td>SR<sub>σ</sub></td><td>23.1</td><td>Sex-ratio genetic variability</td><td>-</td><td>[67]</td></tr><tr><td>SR<sub>a</sub></td><td>-0.0496</td><td>Slope of the effect of temperature on sex-ratio</td><td>°C<sup>-1</sup></td><td>[68]</td></tr><tr><td>SR<sub>b</sub></td><td>27.9</td><td>Effect of temperature on sex-ratio</td><td>-</td><td>[68]</td></tr><tr><td>π<sub>a</sub></td><td>0.0292</td><td>Natural mortality probability</td><td>d<sup>-1</sup></td><td>Fitted on population data</td></tr><tr><td>π<sub>b</sub></td><td>- 0.3820</td><td>Allometric scaling factor</td><td>-</td><td>[17]</td></tr><tr><td>π<sub>c</sub></td><td>0.9576</td><td>Density-independent mortality constant</td><td>d<sup>-1</sup></td><td>[17]</td></tr><tr><td>π<sub>d</sub></td><td>0.0089</td><td>Density-dependent mortality constant</td><td>-</td><td>Fitted on population data</td></tr><tr><td>π<sub>p</sub></td><td>0.0250</td><td>Daily egg predation probability</td><td>d<sup>-1</sup></td><td>[17]</td></tr><tr><td>π<sub>e</sub></td><td>2.839e<sup>-06</sup></td><td>Effect of age on mortality</td><td>d<sup>-1</sup></td><td>[69]</td></tr><tr><td>π<sub>g</sub></td><td>550</td><td>Age threshold of mortality due to aging</td><td>d</td><td>[69]</td></tr><tr><td>F<sub>m</sub></td><td>0.91</td><td>Male appetite modified by male puberty</td><td>-</td><td>This study</td></tr><tr><td>F<sub>inputs</sub></td><td>0.0142</td><td>inputs during the monsoon period</td><td>-</td><td>Fitted on population data</td></tr><tr><td>H<sub>a</sub></td><td>60.9</td><td>age at birth in degree.days</td><td>°C. d<sup>-1</sup></td><td>[34]</td></tr><tr><td>H<sub>b</sub></td><td>10.3</td><td>Threshold of the degree.days</td><td>°C</td><td>[34]</td></tr><tr><td>W<sub>a</sub></td><td>3.205</td><td>Slope of the relationship log(W)/Log(L)</td><td>-</td><td>This study</td></tr><tr><td>W<sub>b</sub></td><td>-5.193</td><td>intercept of the relationship relation log(W)/Log(L)</td><td>-</td><td>This study</td></tr></tbody></table>

doi:10.1371/journal.pone.0125841.t003

median percentage of viable eggs at 24 hours post fecundation per females varied from 40 to 100% (Fig. E in S1 Text).

## DEB and IBM sensitivity analysis

Sensitivity analysis performed on the DEB model (Fig 2 and S5 Text) showed that the main contributors to the DEB model output variations for the growth part of the model were the experimental and reference temperature, the shape coefficient, the maximum surface area specific assimilation rate, the actual ingestion rate divided by the maximal ingestion rate for a body size, the volume specific somatic maintenance costs, and the fraction of energy mobilised from the reserves which is allocated to growth and somatic maintenance for the growth part of the model. The experimental temperature, the reference temperature, the actual ingestion rate divided by the maximal ingestion rate for a body size, and the maximal reproduction rate were the main contributors to the reproduction part of the DEB model. As expected, these
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Fig 2.** Mean of the DEB model sensitivity analysis results for length equations (A) and reproduction equations (B). Parameters are ordered according to the Sobol's first order indices (light gray). Dark gray indices are the Sobol's total indices. First order (Si) and total Sobol' sensitivity indices (STi) were estimated at 25, 50, 75, 100, 200, and 400 dpf for length prediction and at 100, 200, and 400 dpf for the reproduction prediction.</figcaption></figure>

doi:10.1371/journal.pone.0125841.g002

parameters of the DEB model were also in the most influent parameters in the DEB-IBM sensitivity analysis (Fig 3).

The DEB-IBM sensitivity analysis showed that all the important parameters presented both high mean and standard deviation of the Morris' elementary effects (Fig 3 and Fig. B and C in S5 Text). This means that these parameters present non-linear effect or strong interaction with other parameters. Among the most influent parameters, a great part of them are linked to ecological factors. For instance, the photoperiod parameter (P.τ), parameters driving effects of temperature (T<sub>r</sub>, T.τ), and several parameters of the food sub-model (h<sub>p</sub>, K<sub>pr</sub>, and K<sub>s</sub>) have a high impact on the model outputs (Fig 3). In particular, the photoperiod parameter (P.τ) presents a large μ* and σ for total abundance and juvenile frequency at 1095 days (May, just after the beginning of the reproduction). This denotes a strong regulation of the recruitment of juveniles by photoperiod, more than by the water temperature in our model. As expected, the parameters driving the mortality were important (π<sub>c</sub>, π<sub>b</sub>, π<sub>d</sub> and π<sub>a</sub>; Fig 3). These parameters impact mainly the total abundance, the male and female frequency and the mean length (Fig. B in S5 Text). Interestingly, the initial parameters of the simulations (Nb<sub>j</sub>, Nb<sub>F</sub>, Nb<sub>M</sub>, and End) have a very low impact on the global index (Fig 3).

## DEB and IBM model calibration

The DEB calibration was performed using three independent MCMC chains. The Gelman and Rubin convergence statistics were inferior to 1.2 for all parameters [48, 49]. This indicates that
<figure><img src="image_7.png" /><figcaption>**Fig 3.** Morris' global index (A), mean of the elementary effects (B) and standard deviation of the elementary effects (C) determined for the IBM of the zebrafish population dynamics. Presented indices were the mean of the indices determined for the total number of fish, the frequency of adult/juvenile, the mean length of adults and juveniles at 1095, 1156, 1217, 1277, 1339, 1400 days from the beginning of the simulations. Sensitivity analysis indices were calculated on the mean of 30 model repetitions for each parameter combination. "Dummy" is a parameter with no influence on the model outputs. The other parameters are presented Tables 1, 2 and 3.</figcaption></figure>

doi:10.1371/journal.pone.0125841.g003

the three chains converged to the same solution. Our DEB model provided a relevant fit of both the growth and the reproduction data for several datasets at several temperatures from both new experiment and the literature (Figs 4 and 5). The *f* parameter, linked to environmental feeding conditions was fitted for each experiment. Prior and posterior distributions of the parameters are presented in Table A in S2 Text except for the *f* parameter which is presented with the relevant figures.

To fit the three IBM parameters (survival parameters and food input during the monsoon), three independent fitting to the data collected by Hazlerigg et al. [17] were done and converged to the same solution. The adjusted values of $\pi_a, \pi_c$, and $F_{in}$ are presented in Table 3. Fig 6 presents the model prediction compared to the data collected by Hazlerigg et al. [17]. The model provided an accurate prediction of these data, outside of the frequencies of the larger individuals (from 29 to 33 mm) in one population monitored (i.e., this population presented a higher calibration distance with the best parameterisation). In this population, larger individuals were more frequent than other length class, contrary to the two other observed populations and the simulated populations. The data reported by Spence et al. [24] were compared to model predictions without any parameter fitting (Fig 7A and 7B). Observation done by Spence et al. [24] are
<figure><img src="image_8.png" /><figcaption>**Fig 4.** DEB model simulations against observed length for fish from various literature experiments. All lengths are presented as standard length (mm). Circles and crosses represent observations, lines represent model predictions. Panel A represents dataset from Lawrence et al. [31], experimental temperature ($T_{exp}$) of 28.5°C, estimated ingestion level ($f$) of 0.54 (circles) and $f = 0.80$ (crosses). Panel B represents dataset from Gómez-Requeni et al. [33], $T_{exp}$ = 28°C, $f = 0.96$. Panel C represents dataset from Schilling [30], $T_{exp}$ = 28°C, $f = 0.52$. Panel D represents dataset from Bagatto et al. [29], $T_{exp}$ = 25°C, $f = 0.47$. Panel E represents dataset from Eaton and Farley [28], $T_{exp}$ = 25.5°C, $f = 0.86$. Panel F represents dataset from Best et al. [32], $T_{exp}$ = 25°C, $f = 0.87$</figcaption></figure>

doi:10.1371/journal.pone.0125841.g004

censored due to the field sampling methods. The length above which the data were censored is supposed to be 15 mm (fifth percentile of the lengths observed by Spence et al. [24]). Thereby, if individuals smaller than 15 mm were excluded from the model predictions (Fig 7B), the evolution through the year of the length distribution as predicted by our model is in good agreement with the length distribution of the sub-sample of the population monitored (Fig 7B). The proportion of individuals excluded from the simulated populations varies from a minimum of 17% during non-breeding period to a maximum of 91% during the breeding period (Fig 8B).

All variables of the model presented an annual cycle with a clear difference between the non-breeding period and the period favourable for reproduction (monsoon period) (Fig 8). The juvenile biomass represented a small part of the biomass (panel A), whereas in abundance, the juveniles represented a large part of the population during the breeding-period (panel B). Fish biomass was larger than autotrophic and heterotrophic biomass (panel D).

## Discussion

One of the main challenges to ecotoxicology and ecological risk assessment of chemicals is to develop tools and strategies to characterize environmental hazard and risk assessment of chemicals for aquatic environment. To this end, we developed a model of the zebrafish population dynamics that combined DEB and IBM models to extrapolate from ecotoxicological data measured on organisms in laboratory conditions to biological endpoints convenient to support ecological risk assessment of toxic compounds. For that purpose, new experimental data on zebrafish growth and reproduction were produced and used to propose alternative hypotheses
<figure><img src="image_9.png" /><figcaption>**Fig 5.** DEB model simulations compared to the experimental data produced in this study. All lengths are presented as standard length (mm). Points represent observations (mean ± SD for length data), lines represent model predictions. (A) Length data, experimental temperature ($T_{exp}$) of 27°C, estimated ingestion level ($f$) of 0.70 (circles), 0.81 (squares), and 0.93 (triangles). (B) Median of the cumulated number of eggs, $T_{exp} = 26°C, f = 0.48$ (points) and $T_{exp} = 29°C, f = 0.93$ (triangles); Error bars represent the first and third quartiles (C) Model predictions against observed length data for males (black circles and solid line) and females (grey squares and dashed line). $T_{exp} = 27°C, f = 0.76$. (D) Model predictions against observed length data for males (black circles and solid line) and females (grey squares and dashed line). $T_{exp} = 27°C, f = 0.99$.</figcaption></figure>

doi:10.1371/journal.pone.0125841.g005

to the zebrafish DEB model published by Augustine et al. [23]. Then a DEB-IBM couple model integrating the main ecological factors (photoperiod, temperature, and food dynamics) was developed.

Our data on the relationship between masses and lengths of juveniles (from the hatching length to near 15 mm) do not support the hypothesis of anisomorphic growth. Hence, we proposed a zebrafish DEB model based on an alternative hypothesis: a limitation of the food/energy intake linked to the length at the hatching. Such food limitation intake modification relative to standard DEB models has already been used for other organisms [35, 50].

Two other modifications of the model were done to reduce the computation costs of the whole DEB-IBM model. The maturity variables were removed by assuming a constant length at birth and puberty. For the length at birth, as the sensitivity analysis showed a very low impact of the size at birth on the population dynamics, the possible deviations due to this simplification compared to full DEB model are clearly insignificant. For the length at puberty, DEB model predicted that this physiological parameter could vary of about 10% for large variation of food availability [23]. Nevertheless, this limited variation due to food, in addition of the overall small impact of this parameter on the population dynamics showed by the sensitivity analysis (see Fig 5), suggest that this simplification is not of high impact on the population dynamic in our DEB-IBM model.

The second modification concerned the incubation time. An alternative hypothesis was proposed to predict the incubation time, the embryo dynamics have been excluded by using an equation based on the Kimmel et al. [34] empiric formula. According to the data collected by
<figure><img src="image_10.png" /><figcaption>**Fig 6. Probabilistic distributions of the length frequency predicted by the model length compared to frequency distributions observed in [17].**
Circles, crosses and triangles represent the length frequency distribution of the three observed populations. Full lines represent the median length frequency distribution of the simulated populations. Colour level represents the frequency of simulated populations (n = 1,000) having a given percentage of individuals for a given class length. Frequency inferior to < 1e-05 was represented in white. The length class was one millimetre. The frequencies in the populations of the length class frequencies was calculated using class of 0.01.</figcaption></figure>

doi:10.1371/journal.pone.0125841.g006

Kimmel et al. [34], for the zebrafish, the incubation time is mostly driven by the local temperature.

The incubation time and maturity variables simplifications of the model were done to reduce the computation costs of the whole DEB-IBM model while still fits experimental data very well. As the model predictions were accurate, it can be assumed that the modifications present relevant alternative hypothesis to the previous published zebrafish DEB model [23].

Another modification of the DEB food intake was introduced to take into account the sexual dimorphism. Indeed, an appetite factor lowered the male's food intake. This modification was introduced in order to take into account the male zebrafish behaviour modification after the puberty and then the lower male maximal length [24]. Indeed, our data suggested that the energy demand of the males is lower than that of females (males eat less than females), as the males were generally smaller even if both sexes were fed in excess. This was also observed in several other teleost fish [51, 52]. The appetite factor modification is not so far from the type A metabolic acceleration described by Kooijman [53]. The author describes a shift in the surface area-specific assimilation capacity between males and females of the same species at some stage of the development. Some other studies addressing DEB modelling for both males and females in marine species have opted to calibrate several DEB model parameters independently on both males and females [54] or to focus on independent $κ$ parameter for both sexes [55] to explain the sexual dimorphism.
<figure><img src="image_11.png" /><figcaption>**Fig 7.** Zebrafish length distributions predicted compared to the length observed by Spence et al. [24]. (A) Predicted length distributions of all fish. (B) Predicted length distributions of fish > 15 mm. Black full lines represent the median length, purple dotted lines represent the first and third quartile of the fish length, and blue lines represent the limits of 95% of the fish length. Bar plot represent the distribution of the length of the fish sampled by Spence et al. [24] (sample of 120 fish).</figcaption></figure>

doi:10.1371/journal.pone.0125841.g007

The length at puberty was not fitted in our model as we considered that available data did not allow a relevant fit of this parameter. Therefore, this parameter was fixed according to literature.

The sensitivity analysis performed on the zebrafish DEB model showed that experimental conditions are the main contributor to the model output variation. Hence, a particular care should be paid to experimental conditions such as feeding and temperature across the whole life cycle of the organisms in order to fully analyse experimental data with modelling tools, and in particular, ecotoxicological data.

Few data were found to assess the predictability of our model. In addition, these data are censored due to the field sampling methods. Hazlerigg et al. [17] noted that field sampling methods were ineffective in catching smaller individuals. Indeed, less than 5% of fish sampled by Spence et al. [24] measured less than 15 mm. Hence, to compare the predicted length distribution to these reported by Spence et al. [24], the fish smaller than 15 mm have to be excluded. Consequently, length distributions were estimated on a sub-sample of the entire population and at some time points, a large proportion of the population was not observed. These two points resulted in a relative high uncertainty on the length class frequencies estimated and then, could explain some divergences between the model predictions and the observations. In particular, the relative high frequency of the larger individuals (from 29 to 33 mm) in one population observed in [17] contradicted usual fish population length distributions. Indeed, usually, older (larger) individuals are less frequent than younger fish [13].

Annual variation of the different state variables that were predicted by the food dynamic sub-model was classically observed in aquatic systems with similar characteristics to the one we modelled in this study [56–58]. Indeed, fish have a strong impact on the auto and heterotrophic compartments (top-down control). Over the entire year, in our food dynamic sub-model, balance between the nutriment input and output is supposed null to ensure sustainability of
<figure><img src="image_12.png" /><figcaption>**Fig 8.** Zebrafish population dynamics predicted through three years. (A) zebrafish biomass, (B) frequencies of the different fish stages (excluding eggs-larvae), (C) mean length of the different zebrafish stages, and (D) food dynamic state variables.</figcaption></figure>

doi:10.1371/journal.pone.0125841.g008

the ecosystem. Hence, the nutrient loss is compensated by inputs during the monsoon period, which is described as a high food availability period [56]. Actually, during the monsoon the extent of freshwater habitats is maximal with wide variation [56] and runoff conveys a large amount of heterotrophic particles and nutrients from terrestrial sources to aquatic ecosystem [59].

When DEB-IBM model predictions were compared to the observations, the predicted population dynamics seems likely and coherent with ecological knowledge on fish population dynamics. Nevertheless, some processes could be upgraded to improve the mechanistic realism of our model. First, it should be noticed that a lack of precise environmental information regarding the photoperiod and daily water temperature in the natural range of the zebrafish imposed to re-build scenario for these two inputs (linear interpolation or model prediction). However, as sensitivity analysis showed, the parameters which are related to these factors are highly influential on the model outputs (DEB model and IBM). Therefore, more precise information on
these environmental factors and on their impacts on the zebrafish population dynamics should
be acquired.

Contrary to the model developed by Hazlerigg et al. [17], we did not integrate the reduction of mortality rate of juveniles according to their location. Actually, there is no carrying capacity for the vegetated patch, so all juveniles stay in vegetated patches and which led to a global modification of the relationship between length and survival probability. Thus, it does not seem necessary to take into account the impact of the location on the juvenile mortality. Survival process is globally the most uncertain in our model. Indeed, very few data are available to model and calibrate this process. However, the parameters driving the mortality were identified as important in the sensitivity analysis, except the parameters driving the senescence. This is in agreement with the observations in wild populations in which senescence was not observed [56]. It is likely that fish die in natural populations before the senescence process occurs. Globally, more detailed data on the entire wild population structure and dynamics should allow us to more precisely model the survival probability. The relative feeding level is determined assuming that fish feeding has a paralleling mechanism, i.e. the food resource is shared equally to all the fish (proportionally to their need). Nevertheless, the food resources are probably different between juveniles and adults, and dominant individuals attempt to monopolise a food source [56]. In addition, as no mortality would occur in adult female zebrafish after 21 days of starvation [60] and as any other data are available, no explicit modelling of the starvation effect on mortality was integrated. This assumption is realistic, as many fish species evolved the capability to endure prolonged food shortages due to seasonal change in food availability [19]. However, further experiments are needed to understand the mechanism of the food partition and starvation, and to improve this part of our model.

Finally, the annual cycle of zebrafish growth rates should be investigated. Currently, the energy fluxes modelled by the DEB model are identical between the non-breeding period (winter) and the period with active reproduction (spring and summer) and energy use for reproduction is supposed dissipated in winter. No observation on this process was found for the zebrafish. In other fish species, it was observed that non-reproducing females grew more than the productively active females [19]. New data on growth obtained during photoperiodic offset of the reproduction should allow us to propose more precise modelling of the non-breeding period.

Conclusion

A DEB-IBM model has been successfully developed for the zebrafish taking into account the main environmental parameters for the fish population dynamics. The DEB-IBM model was calibrated, deeply analysed, and the predictions were compared to the few published zebrafish population observations. Globally, when the model predictions were compared to the observations, the predicted population dynamics seems highly probable. While, our zebrafish DEB-IBM model can still be improved by acquiring new experimental data on the most uncertain processes (e.g. survival or feeding), it can already serve to predict impact of compounds at the population level and serve as a basis for future works aiming at assessing the ecological impact of chemicals acting on the endocrine system of fish. Hence, the next step will be to compare our model predictions on the impacts of some compounds at the population level with the observed data to test its predictability performances.

**Supporting Information**

S1 File. Model input file. Input file of temperature and photoperiod scenario. This file is called by the model as "InputDataFromMay.txt".
(TXT)
**S2 File. ZIP file archive of the Netlogo model file.**

(ZIP)

**S1 Text. New experimental data.**

(DOCX)

**S2 Text. DEB model calibration.**

(DOCX)

**S3 Text. Models details.**

(DOCX)

**S4 Text. Netlogo code of the zebrafish IBM.**

(DOCX)

**S5 Text. Sensitivity analysis.**

(DOCX)

**Acknowledgments**

Authors would like to thank Enrico Mombelli who contributed to improve the quality of this manuscript. The authors are also very grateful to the three reviewers who greatly helped to improve the content and the clarity of this paper.

**Author Contributions**

Conceived and designed the experiments: RB BP FB ARRP. Performed the experiments: RB BP FB ARRP. Analyzed the data: RB BG SA ARRP. Contributed reagents/materials/analysis tools: RB BG SA JD ARRP. Wrote the paper: RB BG SA BP JD FB ARRP.

**References**

1. Forbes VE, Calow P, Grimm V, Hayashi TI, Jager T. Katholm A, et al. Adding Value to Ecological Risk Assessment with Population Modeling. Human and Ecological Risk Assessment. 2011; 17(2):287-99.

2. Beaudouin R, Dias V, Bonzom JM, Pery ARR. Individual-based model of *Chironomus riparius* population dynamics over several generations to explore adaptation following exposure to uranium-spiked sediments. Ecotoxicology. 2012; 21(4):1225-39. doi: 10.1007/s10646-012-0877-4 PMID: 22396021

3. de Kermoysan G, Joachim S, Baudoin P, Lonjaret M, Tebby C, Lesaulnier F, et al. Effects of bisphenol A on different trophic levels in a lotic experimental ecosystem. Aquatic Toxicology. 2012; 144-145:186-98. http://dx.doi.org/10.1016/j.aquatox.2013.09.034. doi: 10.1016/j.aquatox.2012.11.009 PMID: 23220413

4. Kidd KA, Blanchfield PJ, Mills KH, Palace VP, Evans RE, Lazorchak JM, et al. Collapse of a fish population after exposure to a synthetic estrogen. Proceedings of the National Academy of Sciences of the United States of America. 2007; 104:8897-901. doi: 10.1073/pnas.0609568104 PMID: 17517636

5. Iwasaki Y, Kagaya T, Miyamoto K, Matsuda H, Sakakibara M. Effect of zinc on diversity of riverine ben-thic macroinvertebrates: estimation of safe concentrations from field data. Environmental toxicology and chemistry / SETAC. 2011; 30(10):2237-43. doi: 10.1002/etc.612 PubMed PMID: 21721036.

6. Forbes VE, Calow P. Sibly RM. The extrapolation problem and how population modeling can help. Environmental Toxicology and Chemistry. 2008; 27(10):1987-94. PMID: 19108041

7. Péry ARR, Devillers J, Brochot C, Mombelli E, Palluel O, Piccini B, et al. A Physiologically Based Toxicokinetic Model for the Zebrafish Danio rerio. Environmental Science & Technology. 2014; 48(1):781-90. doi: 10.1021/es404301q

8. Menuet A, Pellegrini E, Anglade I, Blaise O, Laudet V, Kah O, et al. Molecular characterization of three estrogen receptor forms in zebrafish: binding characteristics, transactivation properties, and tissue distributions. Biology of reproduction. 2002; 66(6):1881-92. Epub 2002/05/22. PubMed PMID: 12021076.

9. Brion F, Le Page Y, Piccini B, Cardoso O, Tong SK, Chung BC, et al. Screening Estrogenic Activities of Chemicals or Mixtures In Vivo Using Transgenic (cyp19a1b-GFP) Zebrafish Embryos. Plos One. 2012; 7(5):10. doi: 10.1371/journal.pone.0036069 PubMed PMID: WOS:000305335000023.
10. Vosges M, Le Page Y, Chung BC, Combarnous Y, Porcher JM, Kah O, et al. 17 alpha-Ethinylestradiol disrupts the ontogeny of the forebrain GnRH system and the expression of brain aromatase during early development of zebrafish. Aquatic Toxicology. 2010; 99(4):479-91. doi: 10.1016/j.aquatox.2010.06.009 PubMed PMID: WOS:000282621000006.

11. Cosnefroy A, Brion F, Maillot-Marechal E, Porcher JM, Pakdel F, Balaguer P, et al. Selective Activation of Zebrafish Estrogen Receptor Subtypes by Chemicals by Using Stable Reporter Gene Assay Developed in a Zebrafish Liver Cell Line. Toxicol Sci. 2012; 125(2):439-49. doi: 10.1093/toxsci/kfr297 PubMed PMID: WOS:0002993346000012.

12. Devillers J, Devillers H, Decourtye A, Fourrier J, Aupinel P, Fortini D. Agent-based modeling of the long term effects of pyriproxyfen on honey bee population. In: Devillers J, editor. In Silico Bees. Boca Raton: CRC Press; 2014. p. 179-208.

13. Beaudouin R, Ginot V, Monod G. Improving mesocosm data analysis through individual-based modeling of control population dynamics: a case study with mosquitofish (Gambusia holbrooki). Ecotoxicology. 2012; 21(1):155-64. doi: 10.1007/s10646-011-0775-1 PMID: 21877229

14. 2009. Development and validation of an individual based Daphnia magna population model: The influence of crowding on population dynamics; p. 310-29.

15. Gabsi F, Hammers-Wirtz M, Grimm V, Schäffer A, Preuss TG. Coupling different mechanistic effect models for capturing individual- and population-level effects of chemicals: Lessons from a case where standard risk assessment failed. Ecological Modelling. 2014; 280:18-29.

16. Ashauer R, Agatz A, Albert C, Ducrot V, Galic N, Hendriks J, et al. Toxicokinetic toxicodynamic modeling of quantal and graded sublethal endpoints: A brief discussion of concepts. Environmental Toxicology and Chemistry. 2011; 30(11):2519-24. http://dx.doi.org/10.1002/etc.639. doi: 10.1002/etc.639 PMID: 21805502

17. Hazlerigg CRE, Tyler CR, Lorenzen K, Wheeler JR, Thorbek P. Population relevance of toxicant mediated changes in sex ratio in fish: An assessment using an individual-based zebrafish (Danio rerio) model. Ecological Modelling. 2014; 280:76-88. http://dx.doi.org/10.1016/j.ecolmodel.2013.12.016.

18. Godin J-G. Behavioural Ecology of Teleost Fishes. Godin J-G, editor: Oxford University Press; 1997. 384 p.

19. Wootton RJ. Ecology of Teleost Fishes: Kluwer Academic Publishers, Dordecht, The Netherlands.; 1998. 386 p.

20. Grimm V, Railsback SF. Individual-based Modeling and Ecology. Princeton University Press ed. Princeton, New Jersey 2005. 428 p.

21. Kooijman SALM. Dynamic Energy Budget Theory for Metabolic Organisation: Cambridge University Press; 2010. 532- p. doi: 10.1111/j.1469-185X.2008.00053.x PMID: 19016672

22. Martin B, Jager T, Nisbet R, Preuss T, Hammers-Wirtz M, Grimm V. Extrapolating ecotoxicological effects from individuals to populations: a generic approach based on Dynamic Energy Budget theory and individual-based modeling. Ecotoxicology. 2013; 22:574-83. doi: 10.1007/s10646-013-1049-x PMID: 23430409

23. Augustine S, Gagnaire B, Floriani M, Adam-Guillermin C, Kooijman SALM. Developmental energetics of zebrafish, Danio rerio. Comparative Biochemistry and Physiology—A Molecular and Integrative Physiology. 2011; 159(3):275-83. http://dx.doi.org/10.1016/j.cbpa.2011.03.016. doi: 10.1016/j.cbpa.2011.03.016 PMID: 21440658

24. Spence R, Fatema MK, Ellis S, Ahmed ZF, Smith C. Diet, growth and recruitment of wild zebrafish in Bangladesh. Journal of Fish Biology. 2007; 71(1):304-9. doi: 10.1111/j.1095-8649.2007.01492.x PubMed PMID: WOS:000247441800025.

25. Rasband WS. ImageJ, U. S. National Institutes of Health, Bethesda, Maryland, USA. 2012.

26. Billoir E, Delignette-Muller ML, Péry ARR, Geffard O, Charles S. Statistical cautions when estimating DEBtox parameters. Journal of Theoretical Biology. 2008; 254:55-64. doi: 10.1016/j.jtbi.2008.05.006 PMID: 18571678

27. Kooijman SALM, Bedaux JJM. Analysis of toxicity tests on Daphnia survival and reproduction. Water Research. 1996; 30:1711-23. doi: 10.1016/0043-1354(96)00054-1

28. Eaton RC, Farley RD. Growth and the Reduction of Depensation of Zebrafish, Brachydanio rerio, Reared in the Laboratory. Copeia. 1974; 1:204-9.

29. Bagatto B, Pelster B, Burggren WW. Growth and metabolism of larval zebrafish: effects of swim training. Journal of Experimental Biology. 2001; 204:4335-43. PMID: 11815657

30. Schilling TF. The morphology of larval and adult zebrafish. In: Nüsslein-Volhard C, Dahm R, editors. Zebrafish: A Practical Guide: Oxford University Press Inc; 2002. p. 59-83.

31. 2008. Rapid growth and out-crossing promote female development in zebrafish (Danio rerio); p. 239-46.
32. Best J, Adatto I, Cockington J, James A, Lawrence C. A Novel Method for Rearing First-Feeding Larval Zebrafish: Polyculture with Type L Saltwater Rotifers (*Brachionus plicatilis*). Zebrafish. 2010; 7(3):289-95. http://dx.doi.org/10.1089/zeb.2010.0667. doi: 10.1089/zeb.2010.0667 PMID: 20936984

33. 2010. A reference growth curve for nutritional experiments in zebrafish (Danio rerio) and changes in whole body proteome during development; p. 1199-215.

34. Kimmel CB, Ballard WW, Kimmel SR, Ullmann B, Schilling TF. Stages of Embryonic Development of the Zebrafish. Developmental dynamics. 1995; 203:255-310.

35. Goussen B, Beaudouin R, Dutilleul M, Buisset-Goussen A, Bonzom J- M, Péry ARR. Energy-based modelling to assess effects of chemicals on Caenorhabditis elegans: A case study on uranium. Chemo-sphere. 2014; 120:507-14. doi: 10.1016/j.chemosphere.2014.09.006 PMID: 25278179

36. Zimmer EI, Jager T, Ducrot V, Lagadic L, Kooijman SALM. Juvenile food limitation in standardized tests: a warning to ecotoxicologists. Ecotoxicology. 2012; 21(8):2195-204. Epub 2012/07/31. doi: 10.1007/s10646-012-0973-5 PubMed Central PMCID: PMCPmc3475972. PMID: 22843241

37. Grimm V, Berger U, DeAngelis DL, Polhill JG, Giske J, Railsback SF. The ODD protocol A review and first update. Ecological Modelling. 2010; 221(23):2760-8.

38. Wilensky U. NetLogo. http://ccl.northwestern.edu/netlogo/. Center for Connected Learning and Computer-Based Modeling, Northwestern University. Evanston, IL. 1999.

39. Engeszer RE, Patterson LB, Rao AA, Parichy DM. Zebrafish in the Wild: A Review of Natural History and New Notes from the Field. Zebrafish. 2007; 4(1):21–U126. doi: 10.1089/zeb.2006.9997 PubMed PMID: WOS:000207986600003.

40. Spence R, Smith C. Mating preference of female zebrafish, Danio rerio, in relation to male dominance. Behavioral Ecology. 2006; 17(5):779–83. doi: 10.1093/beheco/arl016

41. Saltelli A, Annoni P, Azzini I, Campolongo F, Ratto M, Tarantola S. Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index2010. 259–70 p.

42. 2007. Estimating the approximation errors when fixing unessential factors in global sensitivity analysis; p. 957–60.

43. Morris MD. Factorial sampling plans for preliminary computational experiments. Technometrics. 1991; 33:161-74.

44. Bois FY, Maszle D. MCSim: a simulation program. Journal of Statistical Software. 1997; 2(9):http://www.stat.ucla.edu/journals/jss/v02/i9 (software available online at www.toxi.ineris.fr).

45. Stonedahl F, Wilensky U. Finding Forms of Flocking: Evolutionary Search in ABM Parameter-Spaces. Proceedings of the MABS workshop at the Ninth International Conference on Autonomous Agents and Multi-Agent Systems; Toronto, Canada.2010.

46. Box GE, Cox DR. An analysis of transformations. Journal of the Royal Statistical Society: Series B. 1964; 26(2):211-52.

47. R Core Team R. R: A Language and Environment for Statistical Computing. 2013.

48. Gelman A, Rubin DB. Inference from Iterative Simulation Using Multiple Sequences. Statistical Science. 1992; 7(4):457-72.

49. Brooks SP, Gelman A. General Methods for Monitoring Convergence of Iterative Simulations. Journal of Computational and Graphical Statistics. 1998; 7(4):434-55.

50. Jager T, ÁLVarez OA, Kammenga JE, Kooijman SALM. Modelling nematode life cycles using dynamic energy budgets. Functional Ecology. 2005; 19:136–44. doi: 10.1111/j.0269-8463.2005.00941.x

51. Beaudouin R, Ginot V, Monod G. Growth characteristics of eastern mosquitofish Gambusia holbrooki in a northern habitat (Brittany, France). Journal of Fish Biology. 2008; 73(10):2468–84. doi: 10.1111/j.1095-8649.2008.02101.x PubMed PMID: ISI:000261786800014.

52. Wootton RJ. Energy allocation in the threespine sticklaback. In: Bell MA, Foster SA, editors. The evolutionary biology of the threespine stickleback. New York: Oxford university press Inc.; 1994. p. 571.

53. Kooijman SALM. Metabolic acceleration in animal ontogeny: an evolutionary perspective. Journal of Sea Research. 2014; 94:128-37.

54. Freitas V, Cardoso JFMF, Lika K, Peck MA, Campos J, Kooijman SALM, et al. Temperature tolerance and energetics: a dynamic energy budget-based comparison of North Atlantic marine species. Phil Trans R Soc B. 2010; 365(1557):3553–65. doi: 10.1098/rstb.2010.0049 PMID: 20921053

55. Bodiguel X, Maury O, Mellon-Duval C, Roupsard F, Le Guellec AM, Loizeau V. A dynamic and mecha-nistic model of PCB bioaccumulation in the European hake (Merluccius merluccius). Journal of Sea Re-search. 2009; 62(2-3):124-34.

56. Spence R, Gerlach G, Lawrence C, Smith C. The behaviour and ecology of the zebrafish, Danio rerio. Biol Rev Camb Philos Soc. 83. England2008. p. 13-34. PMID: 18093234
57. Bala N, Hasan MR. Seasonal fluctuations in water levels and water quality in Oxbow Lakes in relation in fish yields and social conflict: ICLARM, Conf. Proc.; 1999. 280 p.

58. Affan A, Jewel AS, Haque M, Khan S, Lee JB. Seasonal Cycle of Phytoplankton in Aquaculture Ponds in Bangladesh. Algae. 2005; 20:43-52.

59. Rice PJ, Horgan BP. Evaluation of nitrogen and phosphorus transport with runoff from fairway turf managed with hollow tine core cultivation and verticutting. Science of The Total Environment. 2013; 456-457(1):61-8. http://dx.doi.org/10.1016/j.scitotenv.2013.02.051.

60. Drew RE, Rodnick KJ, Settles M, Wacyk J, Churchill E, Powell MS, et al. Effect of starvation on transcriptomes of brain and liver in adult female. Physiological genomics. 2008; 35(3):283-95. doi: 10.1152/physiolgenomics.90213.2008 PubMed PMID: 18728227.

61. Li L, Yakupitiyage A. A model for food nutrient dynamics of semi-intensive pond fish culture. Aquacultural Engineering. 2003; 27(1):9-38. doi: 10.1016/s0144-8609(02)00037-7 PubMed PMID: WOS:000180234500002.

62. Hisaoka KK, Firlit CF. Ovarian cycle and egg production in the zebrafish, *Brachydanio rerio*. Copeia 1962:302-5.

63. Spence R, Smith C. Male territoriality mediates density and sex ratio effects on oviposition in the zebrafish, Danio rerio. Animal Behaviour. 2005; 69:1317-23. http://dx.doi.org/10.1016/j.anbehav.2004.10.010.

64. Markovich ML, Rizzuto NV, Brown PB. Diet Affects Spawning in Zebrafish. Zebrafish. 2007; 4(1):67-74.

65. Balasubramani A, Pandian T. Endosulfan suppresses growth and reproduction in zebrafish. Current Science. 2008; 94(7):10.

66. Uusi-Heikkila S, Wolter C, Meinelt T, Arlinghaus R. Size-dependent reproductive success of wild zebrafish *Danio rerio* in the laboratory. Journal of Fish Biology. 2010; 77:552-69. doi: 10.1111/j.1095-8649.2010.02698.x PMID: 20701640

67. Liew WC, Bartfai R, Lim Z, Sreenivasan R, Siegfried KR, Orban L. Polygenic Sex Determination System in Zebrafish. PLoS ONE. 2012; 7(4):e34397. doi: 10.1371/journal.pone.0034397 PMID: 22506019

68. Sfakianakis DG, Leris I, Mylonas CC, Kentouri M. Temperature during early life determines sex in zebrafish, Danio rerio (Hamilton, 1822). Journal of Biological Research-Thessaloniki. 2012; 17:68-73. PubMed PMID: WOS:000303663000008.

69. Gerhard GS, Kauffman EJ, Wang X, Stewart R, Moore JL, Kasales CJ, et al. Life spans and senescent phenotypes in two strains of Zebrafish (Danio rerio). Experimental Gerontology. 2002; 37:1055-68. http://dx.doi.org/10.1016/S0531-5565(02)00088-8. PMID: 12213556