Bayesian inference of a dynamical model evaluating Deltamethrin effect on Daphnia survival

—The toxicokinetic and toxicodynamic models (TK-TD) are very well-known for their ability, at both the individual and the population level, to accurately describe life cycles such as the growth, reproduction and survival of sentinel organisms under the inﬂuence of an ecological biomarker. Being dynamics, the consistent inference of life history and environmental traits parameters that engender them is sometimes very complex numerically, especially as these parameters vary from one individual to another. In this paper, we estimate the parameters of a survival model TK-TD already applied and validated by the implementation of the R package GUTS (the General Uniﬁed Threshold Model of Survival) by another coding applied to another very recent implementation of Bayesian inference with the R package deBInfer in order to evaluate the survival effects of our ecotoxicological biomarker called Deltamethrin on our Daphnia sample. The study allowed us to evaluate from a population point of view especially the threshold concentration not to be exceeded to observe a survival effect commonly known NEC (No effect Concentration) and possibly determine the correlations between different vari-ables of life history and the environment traits.

Abstract-The toxicokinetic and toxicodynamic models (TK-TD) are very well-known for their ability, at both the individual and the population level, to accurately describe life cycles such as the growth, reproduction and survival of sentinel organisms under the influence of an ecological biomarker. Being dynamics, the consistent inference of life history and environmental traits parameters that engender them is sometimes very complex numerically, especially as these parameters vary from one individual to another. In this paper, we estimate the parameters of a survival model TK-TD already applied and validated by the implementation of the R package GUTS (the General Unified Threshold Model of Survival) by another coding applied to another very recent implementation of Bayesian inference with the R package deBInfer in order to evaluate the survival effects of our ecotoxicological biomarker called Deltamethrin on our Daphnia sample. The study allowed us to evaluate from a population point of view especially the threshold concentration not to be exceeded to observe a survival effect commonly known NEC (No effect Concentration) and possibly determine the correlations between different variables of life history and the environment traits.
The simulation of the temporal evolution of processes leading to toxic effects on organisms is the major role of the use of toxicokinetictoxicodynamic models (TK-TD models) [17]. There is a diversity of TK-TD models for modeling seemingly simple survival according to the underlying assumptions (individual tolerance or stochastic death, speed of toxicodynamic damage recovery, threshold distribution). The General Unified Threshold Model for Survival (GUTS) is the more general survival TK-TD model from which a wide range of existing models can be inferred as special cases [17]. It has special cases of very appropriate model that can be adjusted to the survival data. As a result, it is actively contributing to increasing its application in ecotoxicology research as well as in the assessment of environmental risks related to chemicals.
However, it is known that in toxicokinetics and pharmacokinetics the evolution of xenobiotics (toxic or therapeutic) in a living organism is qualitative and quantitative. By means of a realistic description (ie anatomical, physiological and biochemical) of the absorption processes (inhalation, skin contact, ingestion or intravenous injection), distribution, metabolism and excretion (ADME process), the mechanistic models , which will result, allow the understanding and the simulation of this evolution of the dose of a substance in the various organs and fluids of the body [9]. The action of the organism on the substance defines the toxicokinetics (TK) whereas the opposite effect translates the toxicodynamics (TD). The equations that govern them are differential equations.
To answer why some individuals survive after exposure of chemicals while others die, Ashauer and al., 2015 [2] established the General Unified Threshold Model of Survival (GUTS), a mathematical relationship. In GUTS, there is two assumptions: the threshold of tolerance is individually distributed and that its overcoming causes sudden death among the individuals of a population and the existence of a certain threshold, above which death occurs stochastically, which all people share. As a result, GUTS appeared to be a promising development in the analysis of traditional survival curves and dose-response models.
Recently, Roman Ashauer and al., 2017 [3] treated the paradigm "dose is poison". They illustrated that it is not only the dose that makes the poison but also the sequence of exposure taking into account the toxicokinetic recovery assumptions (the lack of effect that once a chemical is removed from organism) and toxicodynamic recovery (the neglect of the other homeostasis recovery process may be rapid or slow depending on the chemical). To do this, they tested four toxic substances acting on different targets (diazinon, propiconazole, 4,6-dinitro-o-cresol, 4-nitrobenzyl chloride) on the freshwater crustacean Gammarus pulex.
In this study, special consideration is given to the application of Bayesian inference to the evaluation of the effects of Deltamethrin (a pesticide) on a toxicokinetic and toxicodynamic (TK-TD) survival model. Bayesian inference can be a very sophisticated tool for survival data analysis. It is well known for its ability to process data of any sample size, especially small samples as opposed to conventional methods.
Many statistical methods are currently too complex to be fitted using classical statistical methods, but they can be fitted using Bayesian computational methods [14], [23], [28]. However, it may be reassuring that, in many cases, Bayesian inference gives answers that numerically closely match those obtained by classical methods.
In this article, it is mainly to use, from another angle, a new approach to very recent Bayesian implementation [30] allowing the inference of parameters of the model TK-TD GUTS applicable to the adjustment of our survival data collected at the Interdisciplinary Laboratory of Continental Environments (LIEC). It is a very rigorous methodology insofar as it makes it possible to detect the different relations that can exist between the observable quantities of the unobservable quantities, the states and the parameters of the model. Simple to implement, it requires a differential equation TK-TD or DEBtox model, experimental data for the calculation of the likelihood on these data and a prior distribution assumption. A Markov Chain Monte Carlo procedure (MCMC) describes these inputs to estimate the posterior distributions of the parameters and any derived error variables, including model trajectories. This approach is designed with a MCMC diagnosis of inference, the visualization of posterior distributions of the parameters and trajectories of the model used. This manuscript assesses the long-term survival effects of a toxic substance (a pesticide) called Deltamethrin via the use of the highly reputable GUTS model for assessing the survival of living organisms under stressors such as toxic or pharmaceuticals. The plan adopted for the organization of this article is as follows: in the second section (II), we explain the experimental protocol established in the laboratory and present the model TK-TD GUTS used to translate our experimental protocol. In the third section (III), we discuss the results of the Bayesian analysis. We end in section (IV) with a conclusion and discussion.

A. Organism test
One of the three most widely used biological models for the ecotoxicological risk assessment of toxic substances, Daphnia is a major invertebrate of freshwater aquatic ecosystems. The experiments were conducted with clone A of Daphnia magna Straus 1820 (identified by Professor Calow, Uni-versity of Sheffield, United Kingdom). They are more than 40 years old at LIEC (University of Lorraine, France) [38]. Parthenogenetic cultures were carried out in 1L aquaria with LCV medium: a mixture (20/80) of LefevreCzarda (LC) medium and French mineral water called Volvic (V). This medium is supplemented with i) Ca and Mg in order to obtain a total hardness of 250 mg.L −1 and a Ca/Mg molar ratio of 4/1, and ii) a mixture of vitamins (0.1 mL.L −1 ) containing thiamine HCl (750 mg.L −1 ), vitamin B12 (10 mg.L −1 ) and biotin (7.5 mg.L −1 ). Parthenogenetic cultures of daphnids were maintained under a temperature of 20 • C, a photoperiod of 16 − 8 h lightdark and at a density of 40 organism per liter of culture medium. The Daphnia medium was renewed at least three times weekly and daphnids were fed with a mixture of three algal species (5×10 6 Pseudokirchneriella subcapitata, 2.5 × 10 6 Desmodesmus subspicatus, and 2.5 × 10 6 Chlorella vulgaris/Daphnia/day). These algae were also continuously cultivated in the laboratory using a nutrient LC medium.

B. Test chemical
Intensely used in agriculture, Deltamethrin is a class II pyrethroid insecticide that is harmful to freshwater ecosystems, especially the cladoceran Daphnia magna (Straus 1820) [37], [38]. The Deltamethrin (C 22 H 19 Br 2 NO 3 ) used in the experiments is the technical active substance of the formulation DECIS EC25 (25 g.L 1 ) commercialized by Bayer (Germany). Stock solutions were prepared by dissolving the toxicant in acetone immediately prior to each experiment.

C. Data sample
The experimental protocol was carried out during 21 days of observation. Without the control, five different doses of Deltamethrin (9, 20, 40, 80 and 160 ng.L-1, respectively) were administered to Daphnia magna, with a replicate of 10 for each dose submitted. The survivor count has allowed us to summarize our data sample in the table I.

D. Model Used
GUTS is part of mathematical modeling to quantify the temporal evolution of the survival of an organism population, statistically speaking. It is highly reputed for its ability to assess a population survival effects due to a chemical stressor presence (toxicity in other words) responsible for the individuals mortality in this population. Indeed, the toxicokinetic model criterion is explained by the fact that the ingested chemicals will affect a target site within the body before exerting a toxic effect thus causing damage over time. All TK-TD models including a damage state use either the assumption of individual tolerance or SD hypothesis (ie the existence of a single threshold not to be exceeded for all individuals). The modeling assumptions are not the same, it is obviously clear that the results and interpretations that will follow will differ thereafter. Let us not forget that the term "hazard" and specific terms of parametrization of the different models (such as killing rate, recovery rate constant or elimination rate constant) will be misinterpreted in both cases [17]. But GUTS was designed to overcome these confusions because playing a unifying role that merges different concepts of existing models. GUTS is a synthesis of all these models by mixing the aforementioned hypotheses. More complete documentation of GUTS formulation hypotheses can be found in [17]. For all these reasons, we take the GUTS model to adopt it to our survival data study. As in [1], the GUTS model considered is as follows: (1).
Where C(t) represent the toxic dose subjected linearly causing the time course of damage D(t).
The dominant rate constant denoted k e (in days −1 units) models the slowest process inducing the recovery of the exposed organism. In fact, the more slow the recovery in the individual, the more vulnerable he is to the damage. Note that in the body, there are systematically compensation mechanisms and damage repair. The assumption made in this GUTS model is that damage noted D(t) ( damage ) is considered to be the same for all individuals while knowing that once we exceed a certain threshold. The death considered at individual level as a stochastic event will occur and whose probability increases linearly with the damage. At the population level, this threshold is assumed to vary stochastically over the whole population. The hazard rate h z (t) (days −1 ) for individual with threshold z or NEC (No-Effect Concentration) in equation (2) below represents the "instantaneous probability to die" at individual level. The NEC define the concentration threshold not to be exceeded in the body, an amount that we would like to estimate on average. Once it is reached, it affects the health of the living organism.
where the proportionality constant k k (in ng.L −1 .days −1 units) is well known called killing rate and h b (in days −1 units) is the background mortality rate, that is, the control mortality rate, which is assumed to be constant over time [days −1 ]. The equation (3) expresses the probability S(t) that an individual of the population considered will survive until time t conditionally at the threshold z or NEC.

E. Statistical method
In contrast to visual estimation methods, which are often considered biased and not robust,  Bayesian statistics using kinetic data have been very successful over the last two decades [9]. For all these reasons, we use in this paper the Bayesian approach often considered from a practical point of view as a descriptive statistical analysis technique among the others [22]. In Bayesian statistics, any unknown entity is considered as a random variable, in particular parameters of the model used. An assumption of a prior distribution, assigned to each parameter to be estimated, is necessary before the experimental data analysis. Via the famous Bayes theorem, these prior information will be updated with the experimental data in order to retrieve posterior information. Only the Bayesian approach allows to integrate the knowledge that one has of a system by taking advantage of the experimental information [22]. It is a conjunction of the information provided by the probabilistic model by a prior distribution and experimental data. The R package used for our model parameters inferring is deBInfer [30]. We use the R package deSolve [20], [39] as underlined in [30] for the resolution of the implemented TK-TD model. To extrapolate likelihood on our experimental data, we use the Poisson log-likelihood function as defined in the equation (4). The log-likelihood of the data given the parameters, underlying model, and initial conditions is then a sum over the n observations at 1 The Gamma distribution 2 The Beta distribution 3 The Log-normal distribution each time point in t : Here we use small corrections (ec i ) i=1,··· ,6 that are needed because of the differential equations solutions can equal zero, whereas the parameter lambda of the poison likelihood must be strictly positive. We infer them later as suggested in [18], [20]. We set 20, 000 iterations for the MCMC procedure, cnt = 500 worth only 1231.06 seconds of execution with an Intel (R) Core (TM) i3-2350M CPU processor running at 2.30 GHz. The prior distributions assumptions as well as the parameters measures units are presented in the table II.

III. RESULTS AND DISCUSSION
The inference results are presented in tables III and IV. They were obtained using the major functions ode() of the R package deSolve [20] and de_mcmc() of the R package deBInfer [30]. Tables III and IV respectively give the empirical mean and standard deviation for each variable, plus standard error of the mean and the quantiles for each variable.
The threshold concentration above which there are effects on the survival of our test species (Daphnia magna) commonly called NEC is estimated cap 6.042 ± 2.418 ng.L −1 . It is similar to that estimated in one of our studies on the risk assessment of Deltamethrin on growth and reproduction treated separately [10]. This result is  figure 1. In this image, some chain trajectories are reasonable and consistent over time in that their posterior distributions are unimodal, sometimes resembling that of a normal distribution. We can cite for example the parameters NEC and the small corrections ec i=1,··· ,6 . Their prior distributions were those of a log-normal distribution. Unlike the distributions of the k e , k k and h b unimodal parameters, but suspect because they include a large number of outliers. This aberration would confirm the inference complexity of these types of studies. Let's not forget that these are constant rate. The study results are very consistent overall. The figure 2 perfectly shows a lack of detected correlation between parameters. The highest correlation value is 0.36 between k k and NEC, the two most important parameters in GUTS [17]. For proof purposes, we remove a burnin period of 1500 samples and examine parameter correlations in the figure 2 and overlap between prior and posterior densities. The figure 2 reflects the correlation lack between the different parameters of our dynamics evaluating Daphnia survival in the presence of our Deltamethrin stressor.
From the posterior, we simulate 500 trajectories of our TK-TD model while calculating at 95%HDI (Highest Posterior Density Intervals) for the de-  terministic part of the model. HDI sets (intervals) contain all values of the parameter θ such that the posterior density f θ |y is larger than some constant c α , where c α ensures that the coverage probability will be 1 − α. For each exposure concentration, figure 3 shows in the same graph, the experimental data and the model output describing the dynamics of alive Daphnia magna number during the 21 days of experience. It confirms that the inference procedure actually retrieves the model to our data. In addition, the fitted curves are obtained with small estimation errors, see figure 3. Post hoc trajectories adjust very well our observational data for different pesticide doses.

IV. CONCLUSION
This paper is very instructive in that it adapts the GUTS model to our survival data collected at LIEC through a more recent implementation of Bayesian inference (the R library deBInfer). Thus we ignored the use of GUTS (R package GUTS) implementation. With this new R library, it is easy to encode any toxicico-kinetic and toxicodynamic dynamics (TK-TD) then infer the parameters that compose it. Once differential en-coding is complete, the R package deBInfer has a function named de_mcmc() where is integrated that of ode() function of the R package deSolve specially designed for system differential solving such that ordinary, partial or delay differential equations. These last facilitate access to a lot of users types whether they are specialists in the field or not. Most of the life phenomena are modeled using Ordinary Differential Equations (EDO), Partial Differential Equations (PDE), or the Delay-Differential Equations (DDE). As a result, this R package deBInfer facilitates the transition from determinism to stochastic. As part of our study, it allowed us to consistently address our survival analysis with the GUTS TK-TD model use. It really facilitated the manipulation and inference of the parameters of a mechanistic model to describe the bioaccumulation kinetics and dynamics of survival effects in a contaminated environment of the pesticide Deltamethrin.

ACKNOWLEDGMENTS
This study was funded by the French cooperation and the African Center of Excellence in Mathematics, Computer Sciences and TIC (CEA-