Article pubs.acs.org/est
Constructing Time-Resolved Species Sensitivity Distributions Using a Hierarchical Toxico-Dynamic Model Guillaume Kon Kam King,*,†,‡ Marie Laure Delignette-Muller,*,†,‡,§ Ben J. Kefford,∥ Christophe Piscart,⊥ and Sandrine Charles†,‡,# †
Université de Lyon, F69000 Lyon, France Université de Lyon, Université Lyon 1, CNRS, UMR 5558, Laboratoire de Biométrie et Biologie Evolutive, Bâtiment Grégor Mendel, 43 Boulevard du 11 novembre 1918, F69622 CEDEX Villeurbanne, Villeurbanne, France § VetAgro Sup Campus Vétérinaire de Lyon, F69280 Marcy l’Etoile, France ∥ Institute for Applied Ecology, University of Canberra, Canberra, ACT 2601, Australia ⊥ UMR CNRS 6553, Ecobio, Université de Rennes 1, 263 Avenue du Général Leclerc, 35042 Rennes, France # Institut Universitaire de France, 103 boulevard Saint-Michel, 75005 Paris, France ‡
S Supporting Information *
ABSTRACT: Classical species sensitivity distribution (SSD) is used to assess the threat to ecological communities posed by a contaminant and derive a safe concentration. It suffers from several well-documented weaknesses regarding its ecological realism and statistical soundness. Criticism includes that SSD does not take time-dependence of the data into account, that safe concentrations obtained from SSD might not be entirely protective of the target communities, and that there are issues of statistical representativity and of uncertainty propagation from the experimental data. We present a hierarchical toxicodynamic (TD) model to simultaneously address these weaknesses: TD models incorporate time-dependence and allow improvement of the ecological relevance of safe concentrations, whereas the hierarchical approach affords appropriate propagation of uncertainty from the original data. We develop this model on a published data set containing the salinity tolerance over 72 h of 217 macroinvertebrate taxa, obtained through rapid toxicity testing (RTT). The shrinkage properties of the hierarchical model prove particularly adequate for modeling inhomogeneous RTT data. Taking into account the large variability in the species response, the model fits the whole data set well. Moreover, the model predicts a time-independent safe concentration below that obtained with classical SSD at 72 h, demonstrating under-protectiveness of the classical approach.
1. INTRODUCTION Critical for managing chemical stressors in the environment is the ability to assess the risk of particular chemicals and complex effluents to ecological communities. Approaches to assessing the risk for a biotic community of multiple species exposed to a stressor vary in complexity, data requirements, and ecological realism. Applying assessment factors to the acute toxicity value of one sensitive species gives a first approximate idea of an arbitrary safe concentration, whereas a more accurate characterization requires extensive testing of many species as well as studies in model ecosystems (i.e., mesocosms). Species sensitivity distribution (SSD) stands between the two extremes of assessment factors and model ecosystems. It relies on sensitivity data from a few tested species, which are assumed to be a sample of a community of species.1 From this sample, an estimation of the distribution of the sensitivities within the community is derived. This approach is an improvement on safety factors because it removes some arbitrariness in the © XXXX American Chemical Society
choice of the factor and of the sensitive species. SSD remains a manageable approach for many contaminants because it does not require data collected in complex model ecosystems with interacting species; characterizing the response of species tested independently is enough. For these reasons, SSD has become the recommended method for many environmental regulatory bodies such as ECHA,2 ANZEC,3 and the U.S. EPA.4 Classically, SSD takes as input toxicity thresholds, also called critical effect concentrations (CECs). Those are typically an x% effective concentration (ECx) (i.e., the concentration for which a species suffers an x% reduction compared to the control response). The output of the SSD is a hazardous concentration Received: April 28, 2015 Revised: September 18, 2015 Accepted: September 25, 2015
A
DOI: 10.1021/acs.est.5b02142 Environ. Sci. Technol. XXXX, XXX, XXX−XXX
Article
Environmental Science & Technology for p% of the community (HCp), the concentration below the CEC for 1 − p% of the species in the community. The SSD approach is not exempt from criticism.1 Among those are issues of both ecological realism and statistical concerns. The first criticism regarding ecological realism is that SSD does not account for time.5 The CECs are usually determined for one exposure time, which may be short relative to the life cycle of the species. SSD can only be relevant for exposures of that particular time-scale. A second remark concerning ecological realism bears on the difficulty to interpret the HCp. Whether an HCp truly protects 1 − p% of the community crucially depends on the type of CEC. An HC5 based on EC50’s is a concentration for which 95% of the species are affected at a level between 0 and 50%. At such an HC5, a community might suffer strong adverse effects if many species are affected at 45%. Several authors have advocated that SSD should be based on no observed effect concentrations (NOECs) to address this issue. However, the use of this type of CEC has been disparaged extensively because they are based on a wrong interpretation of statistical tests (i.e., no statistically significant effect does not mean no effect), are strongly dependent on the experimental setting, and favor poor resolution on the concentration scale.6−9 Criticisms on statistical grounds include: the sample of species tested to infer the sensitivity distribution of the community is not representative of known communities and not random,10 sample sizes are usually small, and uncertainty on the estimated CECs is generally not taken into account.11−13 Various approaches have been proposed to improve SSD. Time dependence has been addressed by extrapolating longterm effects from short-term exposure through the use of acuteto-chronic (ACR) transformation14−16 or by specifically modeling the toxicity over time with a dynamic energy budget ecotoxicological model (DEBtox).17,18 The potential implications of taking time dependence into account have also been investigated, showing that if species become more sensitive as exposure increases and shifts the mean of the SSD accordingly then it is also critical to consider the evolution of the variance of the SSD to predict the evolution of the HC5.5 A solution proposed for the difficulty of interpreting the HC5 was to use no effect concentration (NEC) models.9,10,19,20 Contrary to HC5 based on a distribution of ECx, HC5 based on a distribution of NEC can be considered as the threshold below which 95% of the species suffer no direct effect for the end point measured. NECs are estimated as parameters of a threshold model and not from statistical tests, so they do not suffer from the same defects as NOECs. They represent a threshold below which the concentration has no effect on the end point considered. Rapid toxicity testing21 (RTT) has been proposed as a way to overcome some of the statistical issues, namely, small sample size, randomness, and representativity. RTT aims to approximate the sensitivity of a large number of species in order to get a representative sample of the community. The focus of RTT is on the estimation of the variability within the community rather than on the precise estimation of the sensitivity of a few species. Yet, RTT raises the issue of including uncertainty on the tested species sensitivity with increased sharpness. Rapid test data have a large uncertainty that is not included in classical SSD. There have been proposals to include that uncertainty by describing species sensitivity as censored data,12,22−25 using a frequentist or a Bayesian approach.13,25,26
There is currently no method which addresses all these issues simultaneously. We propose in this paper a new hierarchical toxico-dynamic (TD) model for SSD that provides a solution to all the aforementioned problems. We model mortality as a stochastic process resulting from the damage induced by the contaminant to an organism. Contaminant concentration and damage are linked by a phenomenological two-compartment model describing all the biochemical and physiological processes leading to mortality.19 Using this TD model, we model explicitly the time dependence of individual species sensitivity using parameters that do not depend on time themselves, arriving at a time independent SSD. The TD model is parametrized using a no effect concentration (NEC). Constructing an SSD using the NEC as end point allows for a convenient interpretation of the HC5. Additionally, the hierarchical approach allows a consistent management of the statistical uncertainty that is propagated from the original raw data to the HC5. The hierarchical model incorporates all the information available on the tested species regardless of the quality of the data, making it easy to include rare species for which it would not be possible to fit a concentration−response model.
2. METHODS 2.1. Data. The data set used for this study was published and described in detail.24 It contains salinity tolerance of 217 macroinvertebrates taxa from southern Murray Darling Basin, Victoria, Australia,27 and Britanny, France.24 The data consist of the mortality of stream invertebrates exposed to artificial seawater with observations of survival after 24, 48, 72, and sometimes 96 h (more details in Supporting Information section S1). We used the raw data (i.e., the number of survivors at each time point) rather than a summary such as the lethal concentration for 50% of the test population (LC50). Data were collected according to the RTT framework28 to approximate the sensitivity of large numbers of field-collected taxa. The data set is strongly inhomogeneous in the number of concentrations tested per taxa (1−98 measured concentrations), in the number of organisms tested per concentration (1−17) and in the number of replicates per treatment (1−3). Notably, there are some common species for which it is possible to estimate an LC50 precisely and rare species for which it is only possible to estimate a censored value for the LC50. The data set includes organisms that disappeared because they were eaten, completed the aquatic phase of their life-cycle, or were otherwise lost and thus could not be followed up. (The latter are hereafter referred to as “lost to follow-up” organisms.) 2.2. Model. 2.2.1. TD Model. The temporal description of the survival data is based on a toxico-kinetic toxico-dynamic (TKTD) model presented in Jager et al.,19 and following their classification, it is a stochastic death model using external concentration as the dose metric. A TKTD model represents the various processes involved in toxicity from the penetration of the contaminant into the organism and its transformations (TK) to its toxic effects on survival (TD). The model presented here aggregates everything related to toxico-kinetics into a damage compartment for reasons explained below. Hence, it is not really a TKTD model but rather a TD model. It contains a mechanistic description of toxic effects on survival over time with four parameters: a NEC threshold, the control mortality (m0), a parameter controlling the delay between exposure and effect (kr), and a parameter controlling how the concentration of toxicant decreases the probability to survive (ks). (A table B
DOI: 10.1021/acs.est.5b02142 Environ. Sci. Technol. XXXX, XXX, XXX−XXX
Article
Environmental Science & Technology with all the symbols is provided in section S2 of the Supporting Information.) The parameters are time-invariant, which makes it possible to predict the mortality outside of the time span covered by the data without the uncertainty exploding and to compute a time-invariant HC5. In particular, it is possible to compute an LCx for any time. As noted by Smit et al.,17 in this model any LCx will decrease asymptotically in time toward the NEC threshold. One traditional presentation uses a onecompartment TKTD model to link the internal concentration to the external concentration (e.g., Kooijman and Bedaux).18 Jager et al.19 report that a simple one-compartment model linking external and internal concentration usually captures the kinetics of whole-body residues well. However, our survival data contains no information on the initial internal concentration, and it cannot be assumed to be negligible because salts are crucial components of organisms. This prevents the direct use of the model described by Jager and colleagues.19,29−31 As explained in section S3 of the Supporting Information, we resorted to modeling damage to the organism instead of internal concentration. Here, damage describes the state of the organism resulting from all the biochemical and physiological processes that are involved in toxicity.19 The mortality rate results from the state of damage of the organism. It is possible to measure the biochemical and physiological parameters and determine how they change with salinity exposure. However, it is more difficult (albeit possible in principle) to establish what specifically is damage, and it could be specific to the taxa considered. If it is demonstrated that the time course of damage reflects the time course of a specific descriptor of a known mode of action, then damage could be interpreted more specifically.19 An advantage of the damage is that it is a concept sufficiently general to allow a single model to account empirically for a variety of processes. We assumed a linear dependence of the hazard rate on the damage, which yielded the following survival function for one species at time t and constant concentration Cw. (See Supporting Information, section S2 for the derivation.) S(t ) = e
w r r ⎡ ⎤ −m0t − ks⎢⎣(Cw − NEC)(t − t NEC) + C r (e−k t − e−k t NEC)⎥⎦ k
Figure 1. Plot of the survival function (eq 1) versus time for various concentrations. (a) Model without kr, (i.e., kr → ∞ ⇒ tNEC → 0). (b) Model with a small value of kr.
After this delay, the effect of salinity increases with time and saturates, which results in a change of concavity in the curves (Figure 1b). More details about the interpretation of the parameters of the TD model are provided in section S4 of the Supporting Information. The parameters characterizing the tolerance of a species to salinity are kr, ks, and NEC. The assumptions in this model are (1) There is no intraspecific variability in sensitivity to the contaminant. (2) There is a species-specific damage threshold below which there is no hazard due to the contaminant. (This assumption leads to the existence of a species-specific NEC, see Jaynes32 for a discussion of this assumption.) (3) The hazard rate for a species is proportional to the internal damage above the threshold plus a constant. This means that it has a timedependent component due to the contaminant plus a constant component due to other causes. This implies that there is no significant hormesis or essentiality. (4) Damage evolves following a one-compartment model, it increases proportionally to the external concentration and decreases proportionally to the damage level. (5) The mechanistic TKTD processes that cause acute and chronic toxicity are the same, and the greater sensitivity often observed with expanded exposure is the result of a buildup of damage. 2.2.2. Conditional Binomial Error Model. Jager et al.19 used a multinomial formulation of the error model on the number of survivors. The multinomial and the conditional binomial error models33 are two mathematically equivalent formulations of the same error model for survival without lost to follow-up organisms (section S5 of the Supporting Information). The presence of lost to follow-up organisms precludes the use of the multinomial model formulation and suggests instead modeling the number of survivors at time conditionally to the number of organisms alive at time tk−1 and that have not disappeared. Replacing the multinomial error model in the presence of lost to follow-up organisms makes the assumption that the mechanism by which they are lost is not related to toxicity
(1)
where S(t) is the probability for an organism to survive until time t, NEC (concentration) is the no-effect concentration, m0 (time−1) is the control mortality per time unit, kr (time−1) is the damage recovery rate and ks ((concentration × time)−1) is the killing rate as defined by Kooijman.18 (See also section S2 of the Supporting Information.) kr induces a delay between exposure and effect. ks controls the magnitude of the effect of salinity on survival. tNEC (time) is the time at which an organism with negligible initial damage starts to feel the effect of salinity and is defined by t NEC = −
1 ⎛⎜ NEC ⎞⎟ ln 1 − kr ⎝ Cw ⎠
(2)
For Cw < tNEC or t < tNEC, the term by which ks is multiplied in eq 1 is set to 0. In our paper, ln x stands for the natural logarithm of x, and log x stands for the base 10 logarithm of x. When the external concentration is below the NEC, only the parameter m0 has an influence on survival (Figure 1 for the lowest concentration). When the external concentration is above the NEC, there is a delay (tNEC) before the internal damage reaches a level affecting the survival of the organism. This delay depends on the parameters kr, NEC, and the external concentration (eq 2). C
DOI: 10.1021/acs.est.5b02142 Environ. Sci. Technol. XXXX, XXX, XXX−XXX
Article
Environmental Science & Technology (and does not influence the estimation of the parameters). The probability for species exposed to concentration Cwj to survive until tk after having lived until tk−1 is Si , j(tk|tk − 1) =
Si , j(tk) Si , j(tk − 1)
(3)
Cwj
where the dependence on was omitted to simplify the notation. We modeled the observed number of survivors Ni,j,k with the following binomial distribution: Ni , j , k ∼ B(Nip, j , k , Si , j(tk|tk − 1))
(4)
Npi,j,k
where is the number of organisms alive at the previous time and that have not disappeared since. Note that Npi,j,k ≠ Ni,j,k−1 because some organisms might have disappeared between tk−1 and tk. 2.2.3. Hierarchical SSD Model. A classical SSD is the distribution of a single summary of species tolerance to a stressor (the CEC). This CEC is usually a function of the parameters of the concentration−response model, such as the LC50, which discards a lot of the information available in the raw data. Specifically, the CEC does not contain any uncertainty from the fit of the model; it is ignored in the computation of the HC5.13 In the case of RTT of rare species, this uncertainty can be very large, and it is not ideal to ignore it. Instead, we chose to model the multivariate distribution of all the parameters of the concentration−response model and their correlations to make the most of the raw data. This way, the uncertainty from the fit of the model can be propagated at the level of the estimation of the community parameters. The shrinkage property (pooling of information across species shrinks together the model parameter estimates) of a hierarchical SSD model is a convenient way of dealing with heterogeneous data, allowing each species to impact the estimation of the community parameters at the extent of the information available on that species. The two levels of hierarchy are the species and the community (Figure 2), where the species play the role of the random effects in a mixed model. Because the log-normal distribution is the distribution most frequently used in classical SSD,34 we followed that custom and assumed a log-normal distribution for each of the parameters. We also model a potential correlation among the parameters with the nondiagonal elements of the correlation matrix. μ and Σ are defined in terms of position, scale, and correlation parameters for each of the four dimensions of the multivariate normal distribution: μ = (μ log k r , μ log ks , μ log m0 , μ log NEC )
(5)
∑ =ξ Λξ
(6)
Figure 2. Probabilistic directed acyclic graphical model of the hierarchical model54 (also called Bayesian network). Ellipses represent variables; rectangles represent covariables. Solid lines represent stochastic links; dotted lines represent deterministic links. To avoid repetition, similar subunits are summarized by plates: the inner set of plates denotes the different times, the intermediate set of plates denotes the different concentrations, and the outer set of plates denotes the different species. Because the graph is directed and acyclic, any two variables are conditionally independent given the value of their parents. log x is used as the logarithm base 10 of x.
The links between the nodes on Figure 2 are θj ∼ 5m(μ , Σ) (stochastic), Si,j(tk|tk−1) = (Si,j(tk)/Si,j(tk|tk−1)) (deterministic), and Ni , j , k ∼ )(Nip, j , k , Si , j(tk|tk − 1)) (stochastic). 2.3. Fit of the Model. The model was implemented using the software Stan35 via the R interface RStan.36 As JAGS37 or WinBUGS,38 Stan performs Bayesian inference, but it uses Hamiltonian Monte Carlo with the No-U-turn sampler39 instead of Markov chain Monte Carlo Gibbs sampling. Stan was chosen for convenience in writing the threshold model and for the availability of various methods to define the prior on the correlation matrix. The code for fitting the model, the full data set, and a file with all the estimated TK parameters are provided in section S8 of the Supporting Information. We used 4 chains with 5000 warm-up iterations and 5000 post-warm-up draws per chain. We chose vague uniform priors for the scale parameters of the normal distributions σlog k r , σlog ks , σlog m0 , σlog NEC ∼
