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 is assumed for
Analysis of target time survival toxicity tests in
morseDR package, 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.
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 \(n_i^0\) of organisms to a concentration \(c_i(t)\) of contaminant (constant or time-variable), and following the number \(n_i^k\) of survivors at time \(t_k\) (with \(t_0 < t_1 < \dots < t_m\) and \(t_0 = 0\)), thus providing a collection \(D = {(c_i, t_k, n_i^k)}_{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 | \(h_b\) | \(m_0\) | \(\text{time}^{-1}\) | SD and IT |
| Dominant toxicokinetic rate constant | \(k_d\) | \(\mbox{NEC}\) | \(\text{time}^{-1}\) | SD and IT |
| Threshold for effects | \(z_w\) | \(m_0\) | \([D]\) | SD |
| Killing rate constant | \(b_w\) | \(k_k\) | \([D]^{-1}\) | SD |
| Median of the threshold effect distribution | \(m_w\) | \(\alpha\) | \([D]\) | IT |
| Shape of the threshold effect distribution | \(\beta\) | \(-\) | \(n.d.\) | IT |
The number of survivors at time \(t_k\) given the number of survivors at time \(t_{k-1}\) is assumed to follow a binomial distribution: \[ N_i^k \sim \mathcal{B}(n_i^{k-1}, f_i(t_{k-1}, t_k)) \] where \(f_i\) is the conditional probability of survival at time \(t_k\) given survival at time \(t_{k-1}\) under concentration \(c_i(t)\). Denoting \(S_i(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 \(S_i(t)\) in GUTS [@jager2011] is given by integrating the instantaneous mortality rate \(h_i\): \[ S_i(t) = \exp \left( \int_0^t - h_i(u)\mbox{d}u \right) \tag{2} \]
In the model, function \(h_i\) is expressed using the internal concentration of contaminant (that is, the concentration inside an organism) \(C^{{\scriptsize 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 \(k_d\) 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 \(f_z(z_w)\) the probability distribution of the no effect concentration threshold, \(z_w\), then the survival function is given by:
\[ S(t) = \int_0^t S_i(t) f_z(z_w) \mbox{d} z_w= \int \exp \left( \int_0^t - h_i(u)\mbox{d} u \right) f_z(z_w) \mbox{d} z_w \]
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 \(z_w\)), and, once exceeded, the instantaneous probability to die increases linearly with the internal concentration. In this situation, \(f_z(z_w)\) 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. \(k_k = + \infty\)), 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 \(F_z\) denotes the cumulative distribution function of \(f_z\).
Here, the exposure concentration \(c_i(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 \(c_i(t)\) is constant, and assuming \(C^{{\scriptsize 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 \(c_i < z_w\), the organisms are never affected by the contaminant:
\[ S_i(t) = \exp( - h_b t ) \tag{3} \]
When \(c_i > z_w\), it takes time \(t^z_i\) before the internal concentration reaches \(z_w\), 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 > t^z_i\), 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 \(h_b\), \(b_w\), \(k_d\) and \(z_w\), we can simulate trajectories by using \(S_i(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 \(f_z\) 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 \(m_w>0\) is the scale parameter (and also the median) and \(\beta>0\) is the shape parameter of the log-logistic distribution.
Posterior distributions for all parameters \(h_b\), \(b_w\), \(k_d\), \(z_w\), \(m_w\) and \(\beta\) 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 \(z_w\) has a high probability to lie within the range of tested concentrations. For each parameter \(\theta\), we derive in a similar manner a minimum (\(\theta^{\min}\)) and a maximum (\(\theta^{\max}\)) value and state that the prior on \(\theta\) 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, \(\theta^{\min}\) and \(\theta^{\max}\) correspond to the 2.5 and 97.5 percentiles of the prior distribution on \(\theta\). For each parameter, this gives: