New Insights to Compare and Choose TKTD Models for Survival

Subscriber access provided by READING UNIV

Article

New insights to compare and choose TKTD models for survival based on an inter-laboratory study for Lymnaea stagnalis exposed to Cd Virgile Baudrot, Sara Preux, Virginie Ducrot, Alain Pavé, and Sandrine Charles Environ. Sci. Technol., Just Accepted Manuscript • DOI: 10.1021/acs.est.7b05464 • Publication Date (Web): 03 Jan 2018 Downloaded from http://pubs.acs.org on January 3, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Environmental Science & Technology is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 30

Environmental Science & Technology

1

New insights to compare and choose TKTD models

2

for survival based on an inter-laboratory study for

3

Lymnaea stagnalis exposed to Cd

4

Virgile BAUDROT‡,a, Sara PREUX‡,a,b, Virginie DUCROTc, Alain PAVEa, Sandrine CHARLESa*

5

a

6

Évolutive, F-69100 Villeurbanne, France

7

b

8

Fédérale de Lausanne EPFL, Lausanne, Switzerland

9

c

10

Univ Lyon, Université Lyon 1, UMR CNRS 5558, Laboratoire de Biométrie et Biologie

School of Architecture, Civil and Environmental Engineering ENAC, École Polytechnique

Bayer AG, CropScience Division, Environmental Safety, Monheim, Germany

TOC-Art

11 12 13

ACS Paragon Plus Environment

1


Page 2 of 30

14

ABSTRACT

15

Toxicokinetic-toxicodynamic (TKTD) models, as the General Unified Threshold model of

16

Survival (GUTS), provide a consistent process-based framework compared to classical dose-

17

response models to analyze both time and concentration-dependent data sets. However, the

18

extent to which GUTS models (Stochastic Death (SD) and Individual Tolerance (IT)) lead to a

19

better fitting than classical dose-response model at a given target time (TT) has poorly been

20

investigated. Our paper highlights that GUTS estimates are generally more conservative and

21

have a reduced uncertainty through smaller credible intervals for the studied data sets than

22

classical TT approaches. Also, GUTS models enable estimating any x% lethal concentration at

23

any time (LCx,t), and provide biological information on the internal processes occurring during

24

the experiments. While both GUTS-SD and GUTS-IT models outcompete classical TT

25

approaches, choosing one preferentially to the other is still challenging. Indeed, the estimates of

26

survival rate over time and LCx,t are very close between both models, but our study also points

27

out that the joint posterior distributions of SD model parameters are sometimes bimodal, while

28

two parameters of the IT model seems strongly correlated. Therefore, the selection between these

29

two models has to be supported by the experimental design and the biological objectives, and

30

this paper provides some insights to drive this choice.

31

KEYWORDS

32

GUTS models, classical dose-response models, constant exposure, aquatic toxicity, x% lethal

33

concentration

34


2

Page 3 of 30


35

INTRODUCTION

36

New chemicals are created every day to be used among others as industrial products, pesticides,

37

pharmaceuticals and so on, for which safety testing may be required in most parts of the world

38

[1]. In the European Union, the REACH (Registration, Evaluation, Authorization and Restriction

39

of Chemicals) regulation has been adopted in 2007 in order to “improve the protection of human

40

health and the environment from the risks that can be posed by chemicals” [2]. It demands to

41

companies the identification and management of risks linked with every substance in the context

42

of their usages. Since 1981, the Organization for Economic Cooperation and Development

43

(OECD), an intergovernmental organization, develops the OECD Guidelines for the Testing of

44

Chemicals, which consist in internationally agreed methods for the testing of chemicals [1].

45

These guidelines are generally matching the expectations of the REACH and other regulations,

46

and are consequently updated and enhanced frequently. For instance, reproductive toxicity tests

47

methods on molluscs are being developed by the OECD since 2011 to support chemical risk

48

assessment [3,4]. Corresponding test guidelines have been published in 2016 [5–7], which

49

include statistical modelling approaches for data analysis.

50

These statistical approaches are derived from the 2006-published OECD guideline, which

51

provides practical guidance on how to analyse the data from toxicity tests [8]. It specifies that

52

dose-response modelling is well adapted to estimate the ECx of a chemical substance (the

53

effective concentration that causes x% of effect, being an LCx if the response is death), as well as

54

its uncertainty. When the dose-response modelling relates to a given exposure duration, we will

55

call this approach a target time (TT) analysis. Several dose-response models exist and are

56

frequently used to analyse collected data from toxicity tests [8].


3


Page 4 of 30

57

As TT analyses do not take into account both time and concentration simultaneously, using them

58

to model data sets depending on these two variables may raise some issues. For instance, in TT

59

analyses, only data at a given TT are used, so that parameters strongly depend on the chosen TT.

60

In addition, the large part of ignored observations may lead to a lower accuracy on parameter

61

estimates. Toxicokinetic-toxicodynamic (TKTD) models, which simulate the time-course of

62

processes leading to toxic effects, may then represent a promising alternative. TKTD models

63

consist of two parts. First, the toxicokinetic part (TK) translates the external concentration of the

64

chemical substance in the medium into an internal concentration within the exposed organism

65

over time. Then, the toxicodynamic part (TD) links this internal concentration to the effects on

66

life-history traits over time [9]. A very large number of TKTD models exist to describe different

67

effects of chemical substances on life history traits. The General Unified Threshold model for

68

Survival (GUTS) has been suggested by Jager et al. (2011), making GUTS a standardized

69

framework unifying TKTD models for survival [9]. Jager’s study showed that a lot of existing

70

survival models can be derived from GUTS as special cases with reduced versions of model

71

equations [9]. Among these special cases, Stochastic Death (SD) and Individual Tolerance (IT)

72

models are probably – nowadays – the most used ones. Stochastic Death models assume that all

73

individuals are identically sensitive to the tested chemical substance above a certain internal

74

threshold concentration and that mortality is a stochastic process once this threshold is reached.

75

On the contrary, Individual Tolerance models are based on the Critical Body Residues (CBR)

76

approach, which assumes that each organism dies as soon as its internal concentration reaches its

77

own internal threshold. The individual internal threshold concentration is then assumed to follow

78

a probability distribution among the organisms [9].


4

Page 5 of 30


79

In the present study, we revisited the data sets from the two ring-tests studied by Ducrot et al. in

80

2014 [3] and Charles et al. in 2016 [5]. These ring-tests were performed in the scope of the

81

consolidation of the OECD Test Guideline for the assessment of chemical effects on Lymnaea

82

stagnalis (the Great Pond Snail) reproduction [3,5,6]. The present study focused however only

83

on the survival of this hermaphrodite freshwater snail when exposed to cadmium (Cd). Although

84

reproduction is not under focus in this paper, the data from the ring-tests were chosen as a basis

85

for our work because they provide a unique collection of survival data over the long-term: this is

86

ideal for our purpose to compare different modelling tools in their ability to address the influence

87

of the interaction between time and concentration on survival during laboratory toxicity tests.

88

The first objective of our paper was indeed to model the time-course of survival of L. stagnalis

89

exposed to Cd, with both GUTS-SD and GUTS-IT approaches, and to compare the performance

90

of the two models. Our second objective was to compare the results of the TKTD approach with

91

the results of the more classical analysis at TT.

92

MATERIALS AND METHODS

93

Principle of the toxicity test

94

To evaluate the impact of Cd on L. stagnalis (simultaneous hermaphrodite), 6 replicates of 5

95

reproducing adults (that is 30 snails per treatment) were exposed to a range of fixed

96

concentrations of Cd, and 6 additional replicates of 5 snails used as controls. Prior to the tests,

97

snails were sampled from a parasite-free cohort, checked for identical size (27 ± 2 mm), shell

98

integrity and reproduction ability [3,5]. As regards Cd, from 5 to 7 concentrations per laboratory

99

were tested, ranging from 16 to 486 µg.L-1 (average measured concentrations). The

100

concentration, the temperature (20°C ± 1°C) and the dissolved oxygen concentration ([ ] > air

101

saturation value) were maintained constant over 56 days long, during which the reproduction and


5


Page 6 of 30

102

the mortality of snails were monitored twice a week [3]. Once counted and recorded, dead snails

103

were removed from the test vessels.

104

This experimental protocol was applied twice: by seven European laboratories between 2011 and

105

2013 (prevalidation ring-test) from which two have been removed because of failure in the

106

protocol (all individual dying whatever the concentration; or non-significant effects were

107

recorded on reproduction, the life-history trait of interest in the study) thus giving five data sets

108

(hereafter denoted Labs 01, 04, 07, 14 and 15) ; and by thirteen laboratories between 2013 and

109

2014 (validation ring-test), among which six (not involved in the prevalidation ring-test) were

110

testing Cd (hereafter denoted Labs 02, 05, 06, 10, 11 and 13). Consequently, we used eleven data

111

sets, five from the prevalidation ring-test, and six from the validation ring-test. The experience of

112

all these laboratories with mollusk testing varied a lot, from inexperienced to expert ones. For

113

more information on the experimental protocol, see [3,5].

114

Significance of the effects of cadmium on survival

115

As done for reproduction [5], we first checked if effects of Cd on survival were significant (in

116

comparison to controls) and concentration dependent. This was performed using the Jonckheere-

117

Terpstra hypothesis test, a nonparametric test determining the statistical significance of trend

118

between ranked variables [10]. This test was run under the R software [11], with the ‘clinfun’ R-

119

package and the ‘jonckheere.test’ function [12]. Control samples were used as reference and the

120

number of iterations was fixed to 106. Results from these analyses showed that the survival rate

121

significantly decreased with increasing Cd concentrations (Jonckheere-Terpstra p-values at day

122

56 < 0.05 as shown in Table S1 in Supporting Information) for all laboratories.

123

Modelling of survival data


6

Page 7 of 30


124

Principle of GUTS models

125

In the following, we detail the mathematical equations of the GUTS approach describing the

126

survival rate of organisms exposed to a constant concentration of Cd over time. During the

127

chronic toxicity tests we studied, concentrations ci ( varying between 1 and 5 or 1 and 7) of Cd

128

were indeed held constant and the mortality was recorded twice a week. Moreover, no

129

measurement of internal concentration was performed during the tests, meaning that the scaled

130

internal concentration is a latent variable as described by the toxicokinetics model part detailed

131

hereafter by equation (1).

132

If we note , the number of survivors at time (with < < ⋯ < and = 0) and

133

concentration , the data sets are composed of triplets = {( , , , )}, of observations. As

134

only the external concentration is available in our data sets, the TK part of GUTS translates

135

first the external concentration into a scaled internal concentration : ()

136

137

=

(

− ())

(1)

being the dominant rate constant, corresponding to the slowest compensating process

138

(damage repair or elimination of the toxicant) dominating the overall dynamics of toxicity [9].

139

The number of survivors , at time given the number of survivors , " at time " is

140

assumed to follow a conditional binomial distribution characterized by the number , " of

141

organisms. The conditional probability for an organism to survive between times " and is

142

given by [13,14]:

143

, ∼ $(, " ,

% (& )

), ≥ 1

% (&'( )

(2)


7


Page 8 of 30

144

+ ( ) is the probability to survive until time under external concentration (with + ( ) =

145

1). The expression of +() depends on the model and will consequently be detailed further for

146

GUTS-SD and GUTS-IT.

147

Specific assumptions of GUTS-SD

148

The GUTS Stochastic Death model (GUTS-SD) supposes, as mentioned in the introduction, that

149

all the organisms have the same internal threshold concentration (denoted , and also called NEC

150

for No-Effect-Concentration), and that, once exceeded, the instantaneous probability to die

151

increases linearly with the internal concentration [9]. At time in our toxicity tests, the

152

internal concentration of Cd is assumed to be null for every organism at every external

153

concentration c in the water, and so no death occurred. In addition, if is above ,, we need to

154

define time - from which the internal concentration, (- ), is equal to threshold ,, leading to

155

death events among the observed organisms [9]. Time - is then defined as dependent on

156

threshold ,:

-

- = − /0 (1 − )

157

(3)

.

158

Once time - is defined, it is possible to write the expression of +(), the probability to survive

159

until time under concentration , as in [9] by:

160

161 162

+() = 123 45 −ℎ(7)879

(4)

With ℎ being the instantaneous probability to die (also called hazard rate) defined as follows: ℎ(7) =

:;2 (0, (7)

− ,) + ℎ>

(5)


8

Page 9 of 30


163

with

164

ℎ> the background hazard rate (i.e., the instantaneous mortality rate in the water control).

165

Equation (5) translates the fact that the instantaneous probability to die increases linearly once

166

() is above ,.

167

Based on [13], if internal concentration is below z, or if time is below - , the survival

168

function is consequently expressed as:

being the killing rate, z the threshold concentration for survival of all the organisms and

+() = 123(−ℎ> )

169 170

171

(6)

Otherwise: +() = 123 (−ℎ> −

(

− ,)( − - ) − & (123(− .

)

− 123(−

- ))

(7)

172

Based on equations (6) and (7), it is possible to determine any ?@A at any time whatever x

173

(once the model parameters are estimated), as being the concentration for which

174 175 176

+() = (1 − 2/100)+( )

(8)

We set 2 = 50 to get the ?@D at any time . Specific assumptions of GUTS-IT

177

For the GUTS-IT model, we assume that the threshold concentration is distributed among the

178

organisms. Here, we assume the distribution to be log-logistic with a median E and a shape F.

179

The probability to survive until time and concentration c can then be expressed as follows:

180

181

+() = 123 (−ℎ> ) (1 −

H(('IJK('&. )) 'M G( ) L

)

(9)

As for GUTS-SD, this expression enables to determine any LCx at any time t whatever x.


9


182

Page 10 of 30

Implementation of the GUTS models

183

To implement the GUTS models (SD and IT), we chose the R-package ‘morse’ [15] that allows

184

the simultaneous use of JAGS and R software thanks to the R-package ‘rjags’ [16]. Package

185

‘morse’ consequently enables the implementation of Bayesian inference.

186

There are multiple ways to implement TKTD models, either using a frequentist approach, see for

187

instance Albert et al., 2016 [17]; or Bayesian inference based on Bayes rule, see for example

188

Delignette-Muller et al., 2017 [13] who illustrated a higher robustness of the Bayesian

189

framework compared to the frequentist one in the case of sparse data sets. In addition, the use of

190

Bayesian inference has already provided promising results in ecotoxicology [18-23].

191

Bayesian inference consists in fitting a probability model to a data set and results in probability

192

distributions of the model parameters as well as of unobserved quantities [24]. Bayesian

193

inference is based on the use of Monte Carlo Markov Chains (MCMC) simulations to infer the

194

parameter estimates of every implemented model. This statistical method uses observations to

195

update what is already known on the parameters (priors), in order to provide their joint posterior

196

distributions [24]. For each one of the three models and each data set, three independent MCMC

197

chains were run in parallel. The number of iterations and the thinning of these chains were

198

determined by using the Raftery and Lewis method [25] as implemented in the ‘rafter.diag’

199

function of the R-package ‘rjags’. We finally checked the chain convergence by using the

200

Gelman statistics [24].

201

For the GUTS-SD model, the choice of the four prior distributions for parameters ℎ> ,

202

z, is based on [13], which favours weakly informative priors in a log scale based on simple

203

assumptions such as the threshold for effects is expected to be in the range of the tested

,

and


10

Page 11 of 30


204

concentrations. So, whatever parameter N, a realistic range of values was defined between a

205

lower value N and a maximum value NOA . From these two extreme values, a log-normal

206

prior distribution was defined as follows:

207

/PQ N ∼ (

RST(U (VW )G RST(U (VWXJ ) RST(U (VWXJ )" RST(U (VW )

,

)

Y

(10)

208

For the GUTS-IT model, concerning the three parameters ℎ> ,

209

distributions were similar to the ones used for the GUTS-SD analysis (equation (10)). The prior

210

distribution of F is however defined differently, by a non-informative log-uniform distribution

211

between -2 and 2.

212

In order to avoid any confusion, the parameters defining each one of the two models are gathered

213

together in Table 1.

214

and E, the chosen prior

Goodness-of-fit

215

Goodness-of-fit of both GUTS models (SD and IT) was assessed by computing the Deviance

216

Information Criterion (DIC) and performing a cross-validation. The DIC is a classical Bayesian

217

measure of predictive accuracy based on posterior estimates. The cross-validation is more robust

218

[24] whereas less used because time-consuming. In cross-validation, models are fitted on a

219

"training" data set, and the resulting parameter estimates are used to calculate the predictions for

220

the other data sets (so-called “validation” data sets). Each of the eleven laboratories was

221

successively used for training and validation, as well as the pool of data from both the pre-

222

validation and the validation ring-tests, leading to a total of thirteen data sets for the cross-

223

validation. For each posterior estimate of the thirteen data sets, we also calculated the prediction

224

for the entire pool of pre-validation and validation data. Then, the goodness-of-fit of all kinds of

225

predictions was quantified from the percentage of observed data lying within the 95% predicted


11


Page 12 of 30

226

credible interval. A prediction is considered as good when around 95% of the observations lie

227

within the 95% predicted credible intervals. The results of the cross-validation are summarized in

228

a 14x13 matrix of these percentages.

229

Target time analysis

230

This classical analysis of survival data at a given target time (denoted TT, usually fixed at the

231

latest observed time point) is based on the log-logistic model [14]. In this modelling approach,

232

the data sets are composed of triplets of observations = {( , , , ,ZZ )} , where ci is the

233

concentration in the media ( varying between 1 and 5 or 1 and 7), , the initial number of

234

individuals at concentration ci and ,ZZ the number of survivors at [[ and concentration ci. This

235

model supposes that the deaths of two organisms are two independent events and that the

236

survival rate after days is a function \ of the concentration. Consequently, the number ,ZZ of

237

survivors at TT follows a binomial distribution [15]: ,ZZ ∼ $(, , \( ))

238 239

(11)

\() being the mean survival rate at TT and concentration c, defined as [14]: \(c) =

240

^ _ I

G4 9

(12)

241

where d, e and b are positive parameters corresponding to the survival rate in the control

242

(parameter d), the LCD (parameter e) and a value related to the effect intensity of the chemical

243

substance (parameter b).

244

By analogy with the GUTS models, and as provided within the R-package ‘morse’, prior

245

distributions of e and b were defined as:


12

Page 13 of 30


/PQ 1 ∼ 4

246

RST(U b ( )cG RST(U bOA ( )c RST(U bOA ( )c" RST(U b ( )c

and

247

,

/PQ d ∼ e(−2,2)

Y

9 (13)

(14)

248

As regards the prior distribution of parameter d, it is assumed to be uniform between 0 and 1.

249

RESULTS

250

GUTS-SD and GUTS-IT analyses

251

We compared the observed and the estimated number of survivors provided as medians and 95%

252

credible intervals for GUTS-SD and GUTS-IT models. This comparison is represented on

253

Figure 1 for Lab.14 and on Figure S1 in Supporting Information for the other laboratories of the

254

two ring-tests and both models. The results obtained for Lab.14 are indeed quite representative of

255

the overall results of all laboratories and of both models. This laboratory will consequently be

256

used to illustrate all outputs hereafter. Whatever the goodness-of-fit of those graphs, all of them

257

reveal very similar patterns for both SD and IT models. In addition, we compared the observed

258

and estimated values by plotting the Posterior Predictive Check (PPC) [15, 26], as shown on

259

Figure S2 in Supporting Information. Here again, GUTS-SD and GUTS-IT appear very similar,

260

as also pointed out by the average percentage of observed points in the 95% credible intervals of

261

the predictions which are 96.32% for both GUTS-SD and GUTS-IT models for Lab. 14. As

262

illustrated by the diagonal of Figure 2, those percentages remain almost the same whatever the

263

laboratory for both GUTS-SD and GUTS-IT models.

264

The goodness-of-fit measured with the Deviance Information Criterion (DIC) indicates that the

265

GUTS-SD model is slightly better (with a mean DIC at 199.94 and a standard deviation at 46.57)

266

than the GUTS-IT model (with a mean DIC at 211.54 and a standard deviation at 52.65). A

267

contrario, based on the cross-validation (Figure 2), pairwise comparisons between all


13


Page 14 of 30

268

laboratories gives mean percentages at 91.55% for the GUTS-SD model and at 92.08% for the

269

GUTS-IT model. This suggests a slightly better fit with the GUTS-IT model. The cross-

270

validation illustrates that predictions from both GUTS-SD and the GUTS-IT models based on

271

parameter estimates from any laboratory or any combination of the laboratories (read columns in

272

Figure 2) are reliable whatever the laboratory for which these predictions are made for.

273

Implementing MCMC simulations to infer all model parameter estimates finally generated a

274

large sample of the joint posterior distribution. The median value of the four parameter estimates

275

and their 95% credible intervals (defined as the range between the 2.5% and 97.5% quantiles) are

276

presented in Table S2 for all laboratories and both models. The joint posterior distribution can

277

also be projected in each plane of parameter pairs in order to visualize the correlation between

278

parameters of each laboratory as shown on Figure S3 in Supporting Information. Finally, we

279

noticed that all inference outputs for the GUTS-SD and the GUTS-IT models were very close

280

whatever the laboratory.

281

We compared common parameters of GUTS-SD and GUTS-IT (

282

laboratories (see Figure S5 in Supporting Information). For each laboratory, the dominant rate

283

constant (

284

sets (Labs 01, 04, 15 and 11). No typical pattern appears when comparing background mortality

285

ℎ> for both models. The no-effect concentration threshold (,) and the median of the threshold

286

distribution (E) show a high similarity between both models, except for Lab.10 that exhibits a

287

large 95% credible interval for , of GUTS-SD compared to E of GUTS-IT.

288

)

,

ℎ> and , with E) for all

is higher in GUTS-SD than in GUTS-IT, with a strong difference for some data

Comparison of GUTS models and target time analysis


14

Page 15 of 30


289

As explained in the Materials and Methods part, it is possible for both GUTS models to estimate

290

the LC50 at any exposure duration , as being the concentration for which the probability to

291

survive at time t is 0.5. LC50 values are also estimated with a TT analysis, as it is exactly

292

parameter 1 (equation (12)).

293

The LC50 estimates (median and 95% credible intervals) under the three approaches are

294

represented on Figure 3 for Lab.14 and on Figure S4 in Supporting Information for the other

295

laboratories (boxplots represent the LC50 at a given target time, and the continuous curves

296

correspond to the GUTS models). We can see on these figures that curves of the GUTS models

297

do not start at day 0: at the beginning of the experiment, there were no dead snails and

298

consequently it is not possible to estimate the LC50 at day 0. We notice again the similarity

299

between GUTS-SD and GUTS-IT models. The target time LC50 estimates are close to GUTS

300

ones, especially from a 35-days exposure duration. After 42 days of exposure, target time LC50

301

estimates do not change anymore, and the three types of estimates are very close both in median

302

and in uncertainty.

303

DISCUSSION

304

TKTD models are known to have several advantages over classical target time (TT) dose-

305

response models mainly due to their mechanistic derivation. In this paper, we highlight the

306

advantages of GUTS models to fit time- and concentration-dependent survival data in

307

ecotoxicology. We show that GUTS models provide smaller 95% credible intervals and more

308

conservative predictions of LCx than classical TT analysis. We also show the impossibility to

309

choose between GUTS-SD and GUTS-IT models based on either the estimated values of model

310

parameters, their joint posterior distributions, their correlations and the time-course of the LC50

311

over 56 days, or based on dedicated tools for goodness-of-fit like DIC and cross-validation.


15


312

Page 16 of 30

TKTD models added-value

313

The main difference between TKTD models and classical TT analyses is that TKTD models take

314

into account both time and concentration simultaneously [9]. Therefore, all the collected data are

315

used by TKTD models (and not only data at a given TT) [27], so that the dynamics of both the

316

exposure and the effects can be better comprehended over time. Importantly, TKTD models

317

enable extrapolating the results beyond the experimental conditions, as soon as the assumptions

318

on the dynamic processes underlying the toxic response are made (in the scope of this paper it

319

means choosing between GUTS-SD and GUTS-IT), as highlighted by Jager et al. [9]. For

320

example, we could have modelled the time-course of survival for longer time, or it would have

321

been possible to simulate time-varying concentration effects even though the concentration was

322

held constant in experiments used to estimate model parameters [27]. More generally,

323

environment displays plenty of situations that have not been tested in laboratories. In this

324

context, obtaining robust predictions of adverse effects is a major goal for regulatory risk

325

assessment. This goal is achievable with TKTD models.

326

Moreover, TKTD models add some mechanistic understanding in comparison with classical TT

327

models, since they are built upon general biological processes [9]. In addition, GUTS parameters

328

are time independent whereas they are time dependent in a TT model. It means that TKTD

329

models give information on the biological processes that take place during the experiment

330

through the biological significance and value of the estimated parameters [28]. Comparing and

331

interpreting the results of different experiments and exposure scenarios is consequently easier

332

and more meaningful with TKTD models [9]. For instance, we can see in Table S2 that the

333

background mortality is the highest for Lab.05 with both GUTS-SD and GUTS-IT approaches.

334

Having this knowledge on background mortality is helpful to disentangle between effects of the


16

Page 17 of 30


335

chemical substance and effects of the experimental conditions. In addition, dominant rate

336

constant

337

the underlying process) and the threshold concentration for effects (either represented by

338

parameters z or α) is related to the no effect concentration on survival. This information is

339

particularly interesting for the understanding of factors that affect sensitivity to a chemical

340

substance in a given species and when comparing the sensitivity between species [28]. Getting

341

similar threshold concentrations (, or E) with both models strengthens our confidence in the

342

reliability of these parameters estimates. In our study, the dominant rate parameter related to the

343

toxicokinetics part (

344

already observed in [29], so that we should cautiously consider the biological interpretation of its

345

parameter value which may depend on several life-history traits (e.g., organism size).

346

Accordingly, comparing both GUTS models has the potential to bring a better understanding of

347

various possible effects and recovery mechanisms [28]. Therefore, this has the potential to bring

348

more realism to the environmental risk assessment.

349

As regards the time-course of the LC50, we can see on Figure 2 and Figure S4 in Supporting

350

Information that GUTS-type approaches often provide more conservative predictions at

351

intermediate exposure durations (between 14 and 38 days depending on the laboratory) than the

352

TT analysis. Also, the 95% credible interval is often smaller with the GUTS models than with

353

the TT analysis. These are important advantages of the GUTS models over the classical analyses,

354

which appear particularly relevant in the field of environmental risk assessment. Indeed,

355

environmental risk assessment could benefit from GUTS-based LC50 values because their

356

estimation appears more precise due to the use of the entire data set. The possibility to compute

357

LC50 values with their credible interval (that is the range of uncertainty) at any time and

is related to the recovery of the organisms (its value enlightens about the speed of

)

was greater for GUTS-SD compared to GUTS-IT. Such a trend was


17


Page 18 of 30

358

concentration is also of particular interest, for example to make predictions in situations where

359

data are lacking, or where new chemical exposure profiles need to be explored.

360

GUTS-SD and GUTS-IT performance

361

We highlighted in the previous section the advantages of the TKTD approaches. In the scope of

362

this paper, we analysed the data with two of them: GUTS-SD and GUTS-IT. As illustrated, the

363

time-course of the LC50 median values and of their 95% credible intervals, as well as the time-

364

course of the predicted number of survivors and the posterior predictive checks are very similar

365

for both models. Results from the cross-validation also reinforce this statement. It is

366

consequently almost impossible to distinguish the two models based on a visual assessment of

367

the fit quality as shown on Figures 1, 2, S1, S2 and S3. These results corroborate a similar

368

statement by Ashauer et al., 2013, for pesticides and several species of fish [29]. The slight

369

differences between the two models are more visible if we look at the estimated parameters.

370

Table S2 and Figure S5 shows that the median of the background mortality ℎ> and of the

371

threshold concentration (represented either by , or E) estimates as well as their 95% credible

372

intervals are close when estimated with GUTS-SD and GUTS-IT (the difference is always below

373

a factor 1.75 for ℎ> and 1.40 for the no-effect-concentration). The median of dominant rate

374

constant

375

times higher for Lab.11). A sensitivity analysis was performed (see all Figures in the section

376

“sensitivity analysis” of Supporting Information) to better understand the relative influence of

377

each parameter on the shape of the predicted survival curve. In this purpose, three parameters

378

over four were fixed and the remaining parameter was varied within ± 50%. This sensitivity

379

analysis showed that

380

the most important influence on the shape on the survival rate over time under our constant

is however often much higher when estimated with GUTS-SD (until nearly 260

and , for GUTS-SD, or

and

E for GUTS-IT, are the parameters with


18

Page 19 of 30


381

exposure conditions. This result agrees with Ashauer et al., who also stressed that the sensitivity

382

of modelled survival towards , and

383

of the exposure patterns this sensitivity could be high [29]. This is because

384

GUTS-SD model are correlated for some laboratories, and

385

all the laboratories of the two ring-tests in the GUTS-IT model (as we can see on Figure S3 in

386

Supporting Information). The consequence of this strong correlation is that changing one of the

387

parameter values has an incidence on the other parameter value, as shown by the fact that some

388

tested couples of (

389

between

390

itself. This could lead us toward the use of the GUTS-SD approach as the preferred method.

391

Nevertheless, if we look at the correlation plot of the GUTS-SD analysis on Figure S3, we can

392

see that some posterior distributions (as for Lab.10) are bimodal. These bimodalities in the

393

posterior distributions appear independently of the number of iterations per MCMC chain. They

394

increase the uncertainty on outputs as the LC50 estimates. However, as we mentioned in the

395

previous “TKTD models added-value” part, the uncertainty with TKTD models we tested can be

396

quantified and still remains smaller in average than the uncertainty of the current classical TT

397

analysis.

398

According to Ashauer et al. [30], the precision and accuracy of some parameters of the GUTS

399

proper model (a combination of both GUTS-SD and GUTS-IT models), depend a lot on the

400

number of individuals per experiment. In their study, the precision and accuracy of the killing

401

rate

402

fixed to 100 instead of 20 in every sample. As there were 30 snails per concentration for each

403

laboratory in our data sets, we could assume that the uncertainty on the killing rate

and

,

differed a lot depending on the exposure pattern: for some

and

and , in the

α are strongly correlated for

E) values were not in the correlation plot line. The fact that the correlation

E is strong for all the laboratories is probably inherent to the GUTS-IT model

and of the median threshold , increased a lot when the initial number of individuals was

for GUTS-


19


Page 20 of 30

404

SD (especially high for laboratories in the validation ring-test) and of the median threshold (, in

405

the GUTS-SD model or α in the GUTS-IT model) could have been reduced by increasing the

406

number of individuals in the experiments.

407

To summarize, it is difficult to select one of the reduced GUTS models presented in this paper as

408

being overall better than the other, as Ashauer et al. [30] or Ducrot et al. [28] noticed already.

409

Both models give close results and their performances depend on the criteria we want to focus

410

on: the GUTS-IT model should allow to avoid bimodal posterior distributions, the GUTS-SD

411

model should allow to avoid strong correlations between parameters, and the use of a TKTD

412

approach provides the time-course of the LC50 with a reduced uncertainty by taking advantage of

413

the whole set of raw data. Biological arguments as the species and/or the chemical substance

414

could motivate the choice towards SD or IT. Nevertheless, based on our results for L. stagnalis

415

exposed to Cd, there was no relevant biological reason to decide. While this could be a matter of

416

concern from a statistical point of view, having two complementary models may be

417

advantageous for environmental risk assessment. Indeed, fitting is fast enough to be performed

418

with both models. Also, checking their goodness-of-fit can easily be handled based on tools as

419

proposed in the present paper: the joint posterior distribution of parameters, the posterior

420

predictive check, the deviance information criterion and the cross-validation. Then, if both SD

421

and IT models fill all goodness-of-fit criteria, as was observed for most of our fits, this testifies

422

the quality of the datasets, and therefore the implemented experimental protocols used to collect

423

them. On the other hand, if both models show bad fit results, the experimental design or the

424

sample size of the datasets could be questioned. In between, if one model appears better than the

425

other, the choice of the most appropriate one should depend on the question at hand.


20

Page 21 of 30


426

Based on the data sets we studied, our results highlight the added-value of TKTD models, which

427

take into account both time and concentration simultaneously, and which are based on the

428

internal processes that take place in response to the exposure to a chemical substance.

429

Considering our results, it was difficult to define one criteria which could help choosing between

430

GUTS-SD or GUTS-IT. Indeed, both models fitted pretty well the observations and each of them

431

had pros and cons as regards the parameter estimates. Consequently, picking one model more

432

than the other depends on the specific requirements and no general rule can be dictated.


21


Page 22 of 30

433

ASSOCIATED CONTENT

434

Supporting Information. The .PDF file provides all complementary fitting results for all

435

laboratories of the two ring-tests, as well as the sensitivity analysis. This information is available

436

free of charge via the Internet at http://pubs.acs.org. Raw data are available upon request from

437

the corresponding author.

438

AUTHOR INFORMATION

439

Corresponding Author

440

* tel.: +33 4 72 43 29 00, email: [email protected]

441

Orcid

442

Sandrine Charles: 0000-0003-4604-0166

443

Author Contributions

444

The manuscript was written through contributions of all authors. All authors have given approval

445

to the final version of the manuscript. ‡ These authors contributed equally.

446

ACKNOWLEDGMENT

447

The authors thank the French National Agency for Water and Aquatic Environments (ONEMA,

448

now denominated French Agency for Biodiversity), the Région Auvergne Rhône-Alpes, the

449

Agence Régionale de Santé Auvergne Rhône-Alpes, the Direction Régionale de

450

l'Environnement, de l'Aménagement et du Logement Auvergne Rhône-Alpes for the financial

451

support. We acknowledge Laurent Lagadic for his support with regard to data generation. We

452

also acknowledge the laboratories which participated in the OECD ring-tests for the development

453

of the OECD test guideline on reproductive toxicity to L. stagnalis for sharing their data i.e.


22

Page 23 of 30


454

AstraZeneca (Brixham, UK), BASF SE (Limburgerhof, DE), Bayer AG (Monheim am Rhein,

455

DE), CEFAS (Lowestoft and Weymouth, UK), FERA (York, UK), Ghent University (Ghent,

456

BE), Goethe University (Frankfurt am Main, DE), Ibacon GmbH (Rossdorf, DE), INRA

457

(Rennes, FR), Fraunhofer IME (Schmallenberg, DE), University of Aveiro (Aveiro, PT),

458

University of Liège (Liège, BE), University of Southern Denmark (Odense, DK), Swedish

459

University of Agricultural Sciences (Uppsala, SW), Texas Tech University, (Lubbock, TX,

460

USA), WIL Research (now denominated Charles River, Ashland, USA.

461

REFERENCES

462

(1) OECD. More about OECD Test Guidelines. URL

463

http://www.oecd.org/env/ehs/testing/more-about-oecd-test-guidelines.htm [cited 12

464

September 2017].

465

(2) European Chemicals Agency. Understanding REACH. URL

466

https://echa.europa.eu/regulations/reach/understanding-reach [cited 12 September 2017].

467

(3) Ducrot V, Askem C, Azam D, Brettschneider D, Brown R, Charles S, Coke M, Collinet M,

468

Delignette-Muller ML, Forfait-Dubuc C, Holbech H, Hutchinson T, Jach A, Kinnberg KL,

469

Lacoste C, Le Page G, Matthiessen P, Oehlmann J, Rice L, Roberts E, Ruppert K, Davis JE,

470

Veauvy C, Weltje L, Wortham R, Lagadic L. 2014. Development and validation of an

471

OECD reproductive toxicity test guideline with the pond snail Lymnaea stagnalis

472

(Mollusca, Gastropoda). Regulatory Toxicology and Pharmacology. 70:605–614.

473

(4) Sieratowicz A, Schulte-Oehlmann U, Wigh A, Oehlmann J. 2013. Effects of test media on

474

reproduction in Potamopyrgus antipodarum and of pre-exposure population densities on

475

sensitivity to cadmium in a reproduction test. Journal of Environmental Science and Health,

476

Part A – Toxic/Hazardous Substances and Environmental Engineering. 48:481–488.


23


Page 24 of 30

477

(5) Charles S, Ducrot V, Azam D, Benstead R, Brettschneider D, De Schamphelaere K, Filipe

478

Goncalves S, Green JW, Holbech H, Hutchinson TH, Faber D, Laranjeiro F, Matthiessen P,

479

Norrgren L, Oehlmann J, Reategui-Zirena E, Seeland-Fremer A, Teigeler M, Thome JP,

480

Tobor Kaplon M, Weltje L, Lagadic L. 2016. Optimizing the design of a reproduction

481

toxicity test with the pond snail Lymnaea stagnalis. Regulatory Toxicology and

482

Pharmacology. 81:47–56.

483 484 485 486 487 488

(6) OECD. 2016. Test No. 243: Lymnaea stagnalis Reproduction Test. OECD Guideline:1–31. doi:10.1787/9789264264335-en. (7) OECD. 2016. Test No. 242: Potamopyrgus antipodarum Reproduction Test. OECD Guideline: 1–23. doi:10.1787/9789264264311-en. (8) OECD. 2006. Current approaches in the statistical analysis of ecotoxicology data: a guidance to application. OECD Environment Health and Safety Publication. 54.

489

(9) Jager T, Albert C, Preuss TG, Ashauer R. 2011. General unified threshold model of survival

490

- A toxicokinetic-toxicodynamic framework for ecotoxicology. Environmental Science and

491

Technology. 45:2529–2540.

492 493

(10) Jonckheere, AR. 1954. A distribution-free k-sample test again ordered alternatives. Biometrika. 41:133–145.

494

(11) R Core Team. 2016. R: A language and environment for statistical computing. R

495

Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/

496

[cited 12 September 2017].

497

(12) Seshan V E. (2016). clinfun: Clinical Trial Design and Data Analysis Functions. R package

498

version 1.0.12. URL https://CRAN.R-project.org/package=clinfun [cited 12 September

499

2017].


24

Page 25 of 30


500

(13) Delignette-Muller ML, Ruiz P, Veber P. 2017. Robust fit of toxicokinetic-toxicodynamic

501

models using prior knowledge contained in the design of survival toxicity tests.

502

Environmental Science and Technology. 51(7):4038–4045.

503

(14) Forfait-Dubuc C, Charles S, Billoir E, Delignette-Muller ML. 2012. Survival data analyses

504

in ecotoxicology: Critical effect concentrations, methods and models. What should we use?

505

Ecotoxicology. 21:1072–1083.

506

(15) Baudrot V, Charles S, Delignette-Muller M L, Duchemin W, Kon-Kam-King G, Lopes C,

507

Ruiz P and Veber P. 2017. morse: MOdelling Tools for Reproduction and Survival Data in

508

Ecotoxicology. R package version 3.0.0. URL

509

https://github.com/pveber/morse/archive/v3.0.0_rc1.tar.gz [cited 02 January 2018].

510 511

(16) Plummer M. 2016. rjags: Bayesian Graphical Models using MCMC. R package version 46. URL https://CRAN.R-project.org/package=rjags [cited 12 September 2017].

512

(17) Albert C, Vogel S, Ashauer R. 2016. Computationally Efficient Implementation of a Novel

513

Algorithm for the General Unified Threshold Model of Survival (GUTS). PLoS

514

Computational Biology. 12:1–19.

515

(18) Billoir E, Delignette-Muller M, Péry A, Charles S. 2008. A Bayesian Approach to

516

Analyzing Ecotoxicological Data. Environmental Science and Technology. 42:879–884.

517

(19) Billoir E, Ene Delhaye HL, Clé Ment B, Delignette-Muller ML, Charles S. 2011. Bayesian

518

modelling of daphnid responses to time-varying cadmium exposure in laboratory aquatic

519

microcosms. Ecotoxicology and Environmental Safety. 74:693–702.

520

(20) Zhang J, Bailer AJ, Oris JT. 2012. Bayesian approach to potency estimation for aquatic

521

toxicology experiments when a toxicant affects both fecundity and survival. Environmental

522

Toxicology and Chemistry. 31:1920–1930.


25


523 524

Page 26 of 30

(21) Zhang J, Bailer AJ, Oris JT. 2012. Bayesian approach to estimating reproductive inhibition potency in aquatic toxicity testing. Environmental Toxicology and Chemistry. 31:916–927.

525

(22) Fox DR. 2010. A Bayesian approach for determining the no effect concentration and

526

hazardous concentration in ecotoxicology. Ecotoxicology and Environmental Safety.

527

73:123–131.

528

(23) Henry Y, Piscart C, Charles S, Colinet H. 2017. Combined effect of temperature and

529

ammonia on molecular response and survival of the freshwater crustacean Gammarus pulex.

530

Ecotoxicology and Environmental Safety. 137:42–48.

531 532

(24) Gelman A, Carlin JB, Stern HS, Dunson DB, Vehtari A, Rubin DB. 2014. Bayesian Data Analysis - Third Edition. CRC Press. 656pp.

533

(25) Raftery AE and Lewis SM. 1995. The number of iterations, convergence diagnostics and

534

generic Metropolis algorithms. In Practical Markov Chain Monte Carlo (W.R. Gilks, D.J.

535

Spiegelhalter and S. Richardson, eds.). London, U.K.: Chapman and Hall. 15pp.

536 537 538 539

(26) Gelman A, Shalizi CR. 2013. Philosophy and the practice of Bayesian statistics. British Journal of Mathematical and Statistical Psychology. 66:8–38. (27) Jager T, Heugens EHW, Kooijman S A LM. 2006. Making sense of ecotoxicological test results: towards application of process-based models. Ecotoxicology. 15:305–314.

540

(28) Ducrot V, Ashauer R, Bednarska AJ, Hinarejos S, Thorbek P, Weyman G. 2016. Using

541

toxicokinetic-toxicodynamic modeling as an acute risk assessment refinement approach in

542

vertebrate ecological risk assessment. Integrated Environmental Assessment and

543

Management. 12:32–45.


26

Page 27 of 30


544

(29) Ashauer R, Thorbek P, Warinton JS, Wheeler JR, Maund S. 2013. A method to predict and

545

understand fish survival under dynamic chemical stress using standard ecotoxicity data.

546

Environmental Toxicology and Chemistry. 32:954–965.

547

(30) Ashauer R, Albert C, Augustine S, Cedergreen N, Charles S, Ducrot V, Focks A, Gabsi F,

548

Gergs A, Goussen B, Jager T, Kramer NI, Nyman A-M, Poulsen V, Reichenberger S,

549

Schäfer RB, Van den Brink PJ, Veltman K, Vogel S, Zimmer EI, Preuss TG. 2016.

550

Modelling survival: exposure pattern, species sensitivity and uncertainty. Nature Scientific

551

Report. 6:29178.

552


27


Page 28 of 30

553 554 555

Figure 1. Observed numbers of survivors (dots) and simulated numbers of survivors associated

556

to the 95% credible band: (A) GUTS-SD model and (B) GUTS-IT model for Lab.14.

557

558 559

Figure 2. Average percentage of observed points in the 95% credible intervals of prediction.

560

Parameters are estimated from data sets of the x-axis, and 95% credible intervals of predictions


28

Page 29 of 30


561

are compared to data sets of the y-axis. Grey levels indicate the percentage of data lying within

562

the 95% credible intervals.

563 564 565

Figure 3. Time-course of the LC50 (black curves) and their 95% credible intervals (grey curves)

566

as estimated at target time (boxplots), with GUTS-SD (continuous lines) and GUTS-IT (dashed

567

lines) models for Lab.14. The observed number of death before day 14 was too low to perform

568

target-time analyses and estimate LC50 values.

569


29


Page 30 of 30

570

Table 1. Parameters and symbols used in the GUTS-SD and GUTS-IT models. Letter d stands

571

for the time unit in days, while (-) stands for dimensionless. Parameter

Symbol Unit

Dominant rate constant

Model

d-1

IT and SD

Background mortality

ℎ>

d-1

IT and SD

Threshold for effects (or no effect concentration)

,

µg.L-1

SD

Killing rate

L.µg-1.d-1 SD

Median of the threshold distribution

E

µg.L-1

IT

Shape of the threshold distribution

F

(-)

IT

572 573


30

New Insights to Compare and Choose TKTD Models for Survival

Recommend Documents