This document describes the statistical models used in morse to analyze survival and reproduction data, and as such serves as a mathematical specification of the package. For a more practical introduction, please consult the Tutorial vignette ; for information on the structure and contents of the library, please consult the reference manual.
Model parameters are estimated using Bayesian inference, where posterior distributions are computed from the likelihood of observed data combined with prior distributions on the parameters. These priors are specified after each model description.
In a survival toxicity test, subjects are exposed to a measured concentration of a contaminant over a given period of time and the number of surviving organisms is measured at certain time points during exposure. In most standard toxicity tests, the concentration is held constant throughout the whole experiment, which we will assume for Analysis of target time survival toxicity tests, but not for Toxicokinetic-Toxicodynamic modeling which can handled time variable exposure. In the case of constant exposure, an experiment is generally replicated several times and also repeated for various levels of the contaminant. For time-variable exposure, a profile of exposure is usually unique, and the experiment is repeated with several profiles of exposures.
In so-called target time toxicity tests, the mortality is usually analyzed at the end of the experiment. The chosen time point for this analysis is called target time. Let us see how this particular case is handled in ‘morse’.
A dataset from a target time survival toxicity test is a collection D = {(ci, niinit, ni)}i of experiments, where ci is the tested concentration, niinit the initial number of organisms and ni the number of organisms at the chosen target time. Triplets such that ci = 0 correspond to control experiments.
In the particular case of target time analysis, the model used in ‘morse’ is defined as follows. Let t be the target time in days. We suppose the mean survival rate after t days is given by a function f of the contaminant level c. We also suppose that the death of two organisms are two independent events. Hence, given an initial number niinit of organisms in the toxicity test at concentration ci, we obtain that the number Ni of surviving organisms at time t follows a binomial distribution: Ni ∼ ℬ(niinit, f(ci)) Note that this model neglects inter-replicate variations, as a given concentration of contaminant implies a fixed value of the survival rate. There may be various possibilities for f. In ‘morse’ we assume a three parameters log-logistic function: $$ f(c) = \frac{d}{1 + (\frac{c}{e})^b} $$ where b, e and d are (positive) parameters. In particular d corresponds to the survival rate in absence of contaminant and e corresponds to the LC50. Parameter b is related to the effect intensity of the contaminant.
Posterior distributions for parameters b, d and e are estimated using the JAGS software [@rjags2016] with the following priors:
The prior on d is chosen as follows: if we observe no mortality in control experiments then we set d = 1, otherwise we assume a uniform prior for d between 0 and 1.
For datasets featuring time series measurements, more complete models can be used to estimate the effect of a contaminant on survival. We assume the toxicity test consists in exposing an initial number ni0 of organisms to a concentration ci(t) of contaminant (constant or time-variable), and following the number nik of survivors at time tk (with t0 < t1 < … < tm and t0 = 0), thus providing a collection D = (ci, tk, nik)i, k of experiments. In ‘morse’, we propose two Toxicokinetic-Toxicodynamic (TKTD) models belonging to the General Unified Threshold model for Survival (GUTS) [@jager2011; @Jager2018GUTSbook]. One is known as the reduced stochastic death model [@nyman2012] or GUTS-SD and the other is the reduced organism tolerance model or GUTS-IT, which we describe now.
| Parameters | Symbols | Alternative symbols | Units | Models |
|---|---|---|---|---|
| Background hazard rate | hb | m0 | time−1 | SD and IT |
| Dominant toxicokinetic rate constant | kd | NEC | time−1 | SD and IT |
| Threshold for effects | zw | m0 | [D] | SD |
| Killing rate constant | bw | kk | [D]−1 | SD |
| Median of the threshold effect distribution | mw | α | [D] | IT |
| Shape of the threshold effect distribution | β | − | n.d. | IT |
The number of survivors at time tk given the number of survivors at time tk − 1 is assumed to follow a binomial distribution: Nik ∼ ℬ(nik − 1, fi(tk − 1, tk)) where fi is the conditional probability of survival at time tk given survival at time tk − 1 under concentration ci(t). Denoting Si(t) the probability of survival at time t, we have: $$ f_i(t_{k-1}, t_k) = \frac{S_i(t_k)}{S_i(t_{k-1})} $$
The formulation of the survival probability Si(t) in GUTS [@jager2011] is given by integrating the instantaneous mortality rate hi: Si(t) = exp (∫0t − hi(u)du)
In the model, function hi is expressed using the internal concentration of contaminant (that is, the concentration inside an organism) $C^{\mbox{${INT}$}}_i(t)$. More precisely: $$ h_i(t) = b_w \max(C^{\mbox{${\tiny INT}$}}_i(t) - z_w, 0) + h_b $$ where (see Table of parameters):
The internal concentration is assumed to be driven by the external concentration, following:
$$ \frac{\mathop{\mathrm{d}\!}C^{\mbox{${\tiny INT}$}}_i}{\mathop{\mathrm{d}\!}t}(t) = k_d (c_i(t) - C^{\mbox{${\tiny INT}$}}_i(t)) \tag{1} $$
We call parameter kd of Eq.(1) the dominant rate constant (expressed in time−1). It represents the speed at which the internal concentration in contaminant converges to the external concentration. The model could be equivalently written using an internal damage instead of an internal concentration as a dose metric [@jager2011].
If we denote fz(zw) the probability distribution of the no effect concentration threshold, zw, then the survival function is given by:
S(t) = ∫0tSi(t)fz(zw)dzw = ∫exp (∫0t − hi(u)du)fz(zw)dzw
Then, the calculation of S(t) depends on the model of survival, GUTS-SD or GUTS-IT [@jager2011].
In GUTS-SD, all organisms are assumed to have the same internal concentration threshold (denoted zw), and, once exceeded, the instantaneous probability to die increases linearly with the internal concentration. In this situation, fz(zw) is a Dirac delta distribution, and the survival rate is given by Eq.(2).
In GUTS-IT, the threshold concentration is distributed among all the organisms, and once exceeded for one organism, this organism dies immediately. In other words, the killing rate is infinitely high (e.g. kk = +∞), and the survival rate is given by: $$ S_i(t) = e^{-h_b t} \int_{\max\limits_{0<\tau <t}(C^{\mbox{${\tiny INT}$}}_i(\tau))}^{+\infty} f_z(z_w) \mbox{d} z_w= e^{-h_b t}(1- F_z(\max\limits_{0<\tau<t} C^{\mbox{${\tiny INT}$}}_i(\tau))) $$ where Fz denotes the cumulative distribution function of fz.
Here, the exposure concentration ci(t) is not supposed constant. In the case of time variable exposure concentration, we use an midpoint ODE integrator (also known as modified Euler, or Runge-Kutta 2) to solve models GUTS-SD and GUTS-IT. When the exposure concentration is constant, then, explicit formulation of integrated equations are used. We present them in the next subsection.
If ci(t) is constant, and assuming $C^{\mbox{${INT}$}}_i(0) = 0$, then we can integrate the previous equation (1) to obtain:
$$ C^{\mbox{${\tiny INT}$}}_i(t) = c_i(1 - e^{-k_d t}) \tag{4} $$
In the case ci < zw, the organisms are never affected by the contaminant:
Si(t) = exp (−hbt)
When ci > zw, it takes time tiz before the internal concentration reaches zw, where: $$ t^z_i = - \frac{1}{k_d} \log \left(1 - \frac{z_w}{c_i} \right). $$ Before that happens, Eq.(3) applies, while for t > tiz, integrating Eq.(2) results in: $$ S_i(t) = \exp \left(- h_b t - b_w(c_i - z_w) (t - t^z_i) - \frac{b_w c_i}{k_d} \left(e^{- k_d t} - e^{-k_d t^z_i} \right) \right) $$
In brief, given values for the four parameters hb, bw, kd and zw, we can simulate trajectories by using Si(t) to compute conditional survival probabilities. In ‘morse’, those parameters are estimated using Bayesian inference. The choice of priors is defined hereafter.
With constant concentration, Eq.(4) provides that $C^{\mbox{${INT}$}}_i(t)$ is an increasing function, meaning that:
$$ \max\limits_{0 < \tau < t} (C^{\mbox{${\tiny INT}$}}_i(\tau)) = c_i(1 - e^{-k_d t}) $$
Therefore, assuming a log-logistic distribution for fz yields:
$$ S_i(t) = \exp(- h_b t) \left( 1 - \frac{1}{1+ \left( \frac{c_i(1-\exp(-k_d t ))}{m_w} \right)^{- \beta}} \right) $$
where mw > 0 is the scale parameter (and also the median) and β > 0 is the shape parameter of the log-logistic distribution.
Posterior distributions for all parameters hb, bw, kd, zw, mw and β are computed with JAGS [@rjags2016]. We set prior distributions on those parameters based on the actual experimental design used in a toxicity test. For instance, we assume zw has a high probability to lie within the range of tested concentrations. For each parameter θ, we derive in a similar manner a minimum (θmin) and a maximum (θmax) value and state that the prior on θ is a log-normal distribution [@delignette2017]. More precisely: $$ \log_{10} \theta \sim \mathcal{N}\left(\frac{\log_{10} \theta^{\min} + \log_{10} \theta^{\max}}{2} \, , \, \frac{\log_{10} \theta^{\max} - \log_{10} \theta^{\min}}{4} \right) $$ With this choice, θmin and θmax correspond to the 2.5 and 97.5 percentiles of the prior distribution on θ. For each parameter, this gives:
In a reproduction toxicity test, we observe the number of offspring produced by a sample of adult organisms exposed to a certain concentration of a contaminant over a given period of time. The offspring (young organisms, clutches or eggs) are regularly counted and removed from the medium at each time point, so that the reproducing population cannot increase. It can decrease however, if some organisms die during the experiment. The same procedure is usually repeated at various concentrations of contaminant, in order to establish a quantitative relationship between the reproduction rate and the concentration of contaminant in the medium.
As already mentionned, it is often the case that part of the organisms die during the observation period. In previous approaches, it was proposed to consider the cumulated number of reproduction outputs without accounting for mortality [@OECD2004; @OECD2008], or to exclude replicates where mortality occurred [@OECD2012]. However, organisms may have reproduced before dying and thus contributed to the observed response. In addition, organisms dying the first are probably the most sensitive, so the information on reproduction of these prematurely dead organisms is valuable ; ignoring it is likely to bias the results in a non-conservative way. This is particularly critical at high concentrations, when mortality may be very high.
In a toxicity test, mortality is usually regularly recorded, i.e. at each time point when reproduction outputs are counted. Using these data, we can approximately estimate for each organism the period it has stayed alive (which we assume coincides with the period it may reproduce). As commonly done in epidemiology for incidence rate calculations, we can then calculate, for one replicate, the total sum of the periods of observation of each organism before its death (see next paragraph). This sum can be expressed as a number of organism-days. Hence, reproduction can be evaluated through the number of outputs per organism-day.
In the following, we denote Mijk the observed number of surviving organisms at concentration ci, replicate j and time tk.
We define the effective observation period as the sum for all organisms of the time they spent alive in the experiment. It is counted in organism-days and will be denoted NIDij at concentration ci and replicate j. As mentionned earlier, mortality is observed at particular time points only, so the real life time of an organism is unknown and in practice we use the following simple estimation: if an organism is alive at tk but dead at tk + 1, its real life time is approximated as $\frac{t_{k+1}+t_k}{2}$.
With this assumption, the effective observation period at concentration ci and replicate j is then given by: $$ NID_{ij} = \sum_k M_{ij(k+1)} (t_{k+1} - t_k) + (M_{ijk} - M_{ij(k+1)})\left( \frac{t_{k+1}+t_k}{2} - t_k \right) $$
In this paragraph, we describe our so-called target time analysis, where we model the cumulated number of offspring up to a target time as a function of contaminant concentration and effective observation time in this period (cumulated life times of all organisms in the experiment, as described above). A more detailed presentation can be found in [@delignette2014].
We keep the convention that index i is used for concentration levels and j for replicates. The data will therefore correspond to a set {(NIDij, Nij)}i of pairs, where NIDij denotes the effective observation period and Nij the number of reproduction output. These observations are supposed to be drawn independently from a distribution that is a function of the level of contaminant ci.
We assume here that the effect of the considered contaminant on the reproduction rate 1 does not depend on the exposure period, but only on the concentration of the contaminant. More precisely, the reproduction rate in an experiment at concentration ci of contaminant is modelled by a three-parameters log-logistic model, that writes as follows:
$$ f(c;\theta)=\frac{d}{1+(\frac{c}{e})^b} \quad \textrm{with} \quad \theta=(e,b,d) $$ Here d corresponds to the reproduction rate in absence of contaminant (control condition) and e to the value of the EC50, that is the concentration dividing the average number of offspring by two with respect to the control condition. Then the number of reproduction outputs Nij at concentration ci in replicate j can be modelled using a Poisson distribution: Nij ∼ Poisson(f(ci; θ) × NIDij) This model is later referred to as Poisson model. If there happens to be a non-negligible variability of the reproduction rate between replicates at some fixed concentrations, we propose a second model, named gamma-Poisson model, stating that: Nij ∼ Poisson(Fij × NIDij) where the reproduction rate Fij at ci in replicate j is a random variable following a gamma distribution. Introducing a dispersion parameter ω, we assume that: $$ F_{ij} \sim gamma\left( \frac{f(c_i;\theta)}{\omega}, \frac{1}{\omega} \right) $$ Note that a gamma distribution of parameters α and β has mean $\frac{\alpha}{\beta}$ and variance $\frac{\alpha}{\beta^2}$, that is here f(ci; θ) and ωf(ci; θ) respectively. Hence ω can be considered as an overdispersion parameter (the greater its value, the greater the inter-replicate variability)
Posterior distributions for parameters b, d and e are estimated using JAGS [@rjags2016] with the following priors:
we assume the range of tested concentrations in an experiment is chosen to contain the EC50 with high probability. More formally, we choose: $$\log_{10} e \sim \mathcal{N} \left(\frac{\log_{10} (\min_i c_i) + \log_{10} (\max_i c_i)}{2}, \frac{\log_{10} (\max_i c_i) - \log_{10} (\min_i c_i)}{4} \right)$$ which implies e has a probability slightly higher than 0.95 to lie between the minimum and the maximum tested concentrations.
we choose a quasi non-informative prior distribution for the shape parameter b: log10b ∼ 𝒰(−2, 2)
as d corresponds to the reproduction rate without contaminant, we set a normal prior 𝒩(μd, σd) using the control: $$ \begin{align*} \mu_d & = \frac{1}{r_0} \sum_j \frac{N_{0j}}{NID_{0j}}\\ \sigma_d & = \sqrt{\frac{\sum_j \left( \frac{N_{0j}}{NID_{0j}} - \mu_d\right)^2}{r_0(r_0 - 1)}}\\ \end{align*} $$ where r0 is the number of replicates in the control condition. Note that since they are used to estimate the prior distribution, the data from the control condition are not used in the fitting phase.
we choose a quasi non-informative prior distribution for the ω parameter of the gamma-Poisson model: log10(ω) ∼ 𝒰(−4, 4)
For a given dataset, the procedure implemented in ‘morse’ will fit both models (Poisson and gamma-Poisson) and use an information criterion known as Deviance Information Criterion (DIC) to choose the most appropriate. In situations where overdispersion (that is inter-replicate variability) is negligible, using the Poisson model will provide more reliable estimates. That is why a Poisson model is preferred unless the gamma-Poisson model has a sufficiently lower DIC (in practice we require a difference of 10).
that is, the number of reproduction outputs during the experiment per organism-day↩︎