The application and limitations of exposure multiplication factors in sublethal effect modelling
The GUTS modelFor completeness, we will prove the result for the most general GUTS model, which will also prove the result for all reduced forms.Theorem 1
For any model version within the GUTS framework, let (S(t; alpha )) denote the survival probability at time t for a given non-zero exposure profile (C_w(t)) scaled by some EMF value (alpha). For any chosen (x > 0) percentage effect (exposure-induced mortality), model end time (t_{E}) and background mortality (h_b) low enough such that (S(t_E; 0) > 0) there exists a unique EMF (alpha _*) such that$$begin{aligned} S(t_{E}; alpha _*) = left( 1 – frac{x}{100}right) S(t_E; 0), end{aligned}$$
(2)
(alpha _*) is the (hbox {LP}_{x}) for the exposure profile.
Baudrot and Charles10 calculated (LC_{50}) values for GUTS-RED-SD and GUTS-RED-IT. Their results implied the result of Theorem 1 for the main regulatory models. Our work makes the result explicit and generalises it to the whole GUTS framework. Another result of Theorem 1 is that the (hbox {LP}_{x}) is monotonically increasing with respect to x. For example, the (hbox {LP}_{50}) will always larger than the (hbox {LP}_{10}) for the same exposure profile. This result comes directly from (S4) in the SI.DEB modelsDue to the additional complexity of the DEB model we split the result into multiple theorems and proofs, starting by showing continuity and monotonicity of the damage ODE.
Theorem 2
Let
(C_w)
be some external concentration over time. Assume an effects model where the effects of higher exposure on growth and/or reproduction are always adverse (or zero) at all points in time. Then, defining the scaled damage ODE as
$$begin{aligned} begin{aligned} frac{D(t; alpha )}{dt} =&k_d(x_u alpha C_w – x_eD) – (x_G + x_R)D. end{aligned} end{aligned}$$
(3)
Then, for any combination of feedbacks (varvec{X} = [varvec{X}_u, 0, varvec{X}_G, varvec{X}_R]), damage is monotonically increasing with respect to (alpha), and is continuous with respect to (alpha) as long as changes to L and R are continuous. Moreover, damage is strictly monotonically increasing with respect to (alpha) whenever (D(t;1) > 0).
The limitation that (varvec{X}_e = 0) will be discussed in greater detail later. However, depending on the pMoA of the stressor we can extend the result of Theorem 2 slightly.
Corollary 3.1
If the pMoA of a substance directly affects reproduction and does not affect growth, i.e. (varvec{S} = [0, 0, 0, s_{R}, s_{H}]) then the results of Theorem 2holds for any combination of feedbacks.
Finally, we can step from the results of Theorem 2 and Corollary 3.1 to show the existence and uniqueness of a critical multiplier ((hbox {EP}_{x}) or (hbox {LP}_{x})) for growth, reproduction and survival.
Theorem 3
Consider the DEB-TKTD model of Jager11 and a substance such that at least one of Theorem 2or Corollary 3.1hold. Further, let (C_w(t)) be a non-zero exposure profile where the time of first exposure is before (t_1) as defined in Table 1. Then, for any chosen percentage effect level (x > 0) there exists a unique EMF (alpha _* >0) such that$$begin{aligned} min left( frac{L(t_{E}; alpha _*)}{L_(t_{E}; 0)}, frac{R_c(t_{E}; alpha _*)}{R_c(t_{E};0)}, frac{S(t_{E}; alpha _*)}{S(t_{E}; 0)}right) = 1 – frac{x}{100} end{aligned}$$
(4)
this (alpha _*) is the (hbox {EP}_{x}) (or (hbox {LP}_{x})) for the exposure profile (C_w(t)).
Table 1 Table of the state variables and pMoAs (including combinations of pMoAs) in the DEB-TKTD model11.Full size table
The monotonicity of effects on all state variables in the DEB model means that, for the conditions described in Theorem 3, the (hbox {EP}_{x}) (or (hbox {LP}_{x})) is also monotonically increasing with respect to x.We should note here that one can either setup an algorithm to find the critical multiplier value for growth, reproduction and survival individually and then select the minimum or setup the algorithm to directly find the minimum critical multiplier as in (4). Both will produce the same result, but the second approach is likely to be faster.One could argue that ERA should consider the combined effects of lethal and sublethal stress on the individual’s fitness. This is possible using the continuous form of the Euler–Lotka equation24$$begin{aligned} B(t) = int _0^t B(t-a) l(a)b(a) da, end{aligned}$$
(5)
where B(t) is the number of births at time t, l(a) is the fraction of females which survive to age a and b(a) is the birth rate for mothers of age a. For the offspring of a test population which all have the same age (as is the standard in long-term toxicity experiments) this integral collapses to a single point, (B(t-a) = 1) when (t=a) and zero elsewhere. The DEB model provides exactly the values which we need to calculate B(t). Namely$$begin{aligned} l(a) = S(a), quad b(a) = frac{d}{dt}R_c(a). end{aligned}$$One can now find the births per individual per time predicted by the DEB model as$$begin{aligned} B(t) = S(t)frac{d}{dt}R_c(t). end{aligned}$$
(6)
Integrating (6) over the duration of the experiment gives the expected number of offspring produced per female alive at the start of the test.There are two clear options for how to proceed. Firstly, one could calculate (int _0^{t_E} B(t) dt) for each EMF and compare it to the control, similar to finding (hbox {EP}_{x}) values for individual endpoints. Alternatively, one can use B(t) as the basis to estimate the intrinsic population growth rate25. This quantity provides an estimation of population growth based on the survival and fecundity over time of individuals. Indeed, it is listed as a potential output value in the experimental guidelines for standard Daphnia magna reproduction tests26. For the first of these options we offer an extension to Theorem 3.
Corollary 3.2
Consider a DEB-TKTD model and exposure profile such that Theorem 3holds. The number of expected offspring per female, given by$$begin{aligned} mathrm {B}(t_E; alpha ) = int _0^{t_E} S(t; alpha )frac{d}{dt}R_c(t; alpha ) dt end{aligned}$$has a unique (hbox {EP}_{x}) (alpha _*) such that$$begin{aligned}frac{mathrm {B}(t_E; alpha _*)}{mathrm {B}(t_E; 0)} = 1 – frac{x}{100}end{aligned}$$
Our results provide a rigid boundary to the applicability domain of the EMF approach both in terms of existence and uniqueness. Existence relies on the initial time in the profile when external concentration is non-zero, as described in Table 1. While it is important to know about these conditions, they will rarely inhibit an ERA, since long initial periods with zero exposure are uncommon.Cases where uniqueness cannot be guaranteed require more caution and it is unwise to use root-finding algorithms. In the next subsection we explore what can happen outside of this domain and provide suggestions for how to still produce a single reliable (hbox {EP}_{x}) value.Surface:volume scaling of eliminationThere is a reason that in Theorem 2, (varvec{X}_e = 0) was specified. In some cases when (varvec{X}_e = 1) a higher multiplier is no guarantee of higher damage for all time. Consider a substance which acts on assimilation and has surface area:volume scaled elimination (i.e. (varvec{X} = [0, 1, 0, 0])). The damage ODE under some EMF (alpha) is then$$begin{aligned} frac{dD}{dt} = k_d left( alpha C_w – frac{L_m}{L} D right) , end{aligned}$$where (L_m) is the maximum length the organism can reach. The EMF has a positive direct effect on damage, but also an opposing indirect effect. Increasing damage decreases the size of the organism which, due to the surface area:volume elimination of damage, enables faster elimination of damage. As a result, not only does Theorem 2 no longer hold but in fact a larger multiplier value can cause lower damage at some points in an exposure profile. In other words, we observe a paradoxical result whereby more exposure translates to less effect some time after exposure.Figure 3 illustrates what we will refer to as the “more is less” scenario. The exposure consists of a single pulse early in the animal’s life, modelled for two multiplier values, (alpha _2 > alpha _1). During the exposure phase the direct effect of the higher exposure causes higher damage and greater effects on size. After the pulse, external exposure is zero, and therefore the external concentration and uptake remain zero regardless of (alpha). Regardless of the EMF, scaled damage can only decrease during this phase. However, the effects of the higher multiplier are still relevant. As Fig. 2 shows, the feedback processes still influence damage dynamics. The model organism exposed to (alpha _2C_w) is smaller and therefore able to eliminate damage more rapidly because (varvec{X}_e = 1). This eventually leads to lower damage for the model organism exposed to (alpha _2C_w) (i.e. (D(t; alpha _1) > D(t; alpha _2))). The more is less phenomenon can also impact growth and cumulative reproduction, as seen in Fig. 3b,c. Sometime after exposure (L(t;alpha _2) > L(t; alpha _1)) and (R(t;alpha _2) > R(t; alpha _1)). For survival, and any additional endpoints without recovery, this “crossover” is unlikely, mortality during the exposure phase (where (D(t; alpha _2) > D(t; alpha _1))) will almost certainly dominate any mortality during the recovery phase. Figure 3d shows that for certain (x%) effect levels (vertical axis) multiple (hbox {EP}_{x}) values exist.Figure 3An illustration of the issues which can occur using the EMF approach for substances with surface area:volume scaled elimination (i.e. (varvec{X} = [0, 1, 0, 0])). The (non-multiplied) exposure is a constant (1 mu g/L) for the first 14 days and zero thereafter and effects assimilation only ((varvec{S} = [1, 0,0, 0, 0])). (a) Scaled damage, (b) length over time, (c) cumulative reproduction. (d) Endpoint value as a proportion of control after 40 days. The shape of these curves show that certain effect levels can be caused by two distinct multiplier values. Parameter values are (L_0 = 0.1), (f = 1), (r_B = 0.1), (L_p = 0.6), (L_m = 1), (R_m = 15), (kappa = 0.8), (y_P = 0.64) (z_b = 0.1), (b_b = 1), (k_d = 0.05), (varvec{X} = [0, 1, 0, 0]). See the SI for the definitions of these parameter values.Full size imageIn practice, instances of non-uniqueness such as Fig. 3 will be rare since they rely on a sudden and significant decrease in external exposure. Moreover, EMF methods for DEB-TKTD models will include a moving time window method18 consisting of many exposures constructed sequentially and assessed. Each window will produce an (hbox {EP}_{x}) value, but only the lowest will be relevant for the ERA. A time window which starts slightly earlier in the broader exposure profile would feature the same pulse later in the model organism’s lifespan and thus not allow organism recovery. Depending on the exact endpoint used, one would expect those windows to have a lower (and unique) (hbox {EP}_{x}). However, the potential for multiple (hbox {EP}_{x}) values raises concerns across all areas which impose a multiplicative margin of safety. We cannot guarantee that a multiplier resulting in (x%) effects exists nor that any value found by the algorithm is unique.Although not pictured here, maintenance and growth pMoAs and combinations of feedbacks which include (varvec{X}_e = 1) can also produce the “crossover” in the damage values and the “more is less” phenomenon seen in Fig. 3. It can also arise for scenarios which do not feature a deviation from the standard rules for growth (e.g. a starvation phase) and for other DEB based models. The SI features a similar plot to Fig. 3 showing damage crossover for a standard DEB model.Knowing this, the obvious question is how to proceed? Certainly with caution when (varvec{X}_e = 1) is necessary in model calibration and validation. Under such circumstances algorithms must ensure that the (hbox {EP}_{x}) value found is the lowest multiplier which gives (x%) effects when there is a risk of non-uniqueness. The brute force approach, incrementing from zero until the desired effect level is met or exceeded, is one example. Whether it is realistic for higher EMF values to cause reduced effects in vivo then does not alter the conservatism of the approach for ERA.Table 2 summarises the domain where the margin of safety approach can be used in conjunction with a root-finding algorithm without concern in the DEB-TKTD model of Jager11. For model configurations where non-uniqueness could emerge using another method to find the (hbox {EP}_{x}) is advisable. For example, a brute-force approach starting from an EMF of 0 in small increments (e.g. by 0.1). Without good reason, calibration should first be attempted with no feedbacks. Under this guiding philosophy of pursuing model simplicity we expect that the problem cases will be rare.Table 2 A table to mark under which scenarios the EMF approach is and is not guaranteed to produce a unique (hbox {EP}_{x}).Full size tableOther issuesThe damage crossover illustrated in the previous subsection occurs more commonly, and to a greater extent, when the pMoA is assimilation effects. This is because, at least in this standard implementation, stress can cause (100%) effect and completely cease assimilation when (s_A ge 1) (see SI for details). When this is the case, higher exposure (even from an increased multiplier) does not translate to higher stress. This differs from other pMoAs, whose stress values are unbounded. Indeed, replacing (1 – s_A) with (1/(1 + s_A)) in the model ((S5) in the SI) reduces the occurrence and scale of “crossovers” such as Fig. 3. However, the formulation of the pMoA should not be based on how it might affect the algorithm or the EMF.Certain species require further deviations from the standard model. For instance, different life-stages, growth and/or reproduction rules might be introduced to explain observed phenomena. Before models featuring these deviations are used in an EMF approach one should consider the potential issues as we have done in this section. While a proof of existence and uniqueness of the (hbox {EP}_{x}) for each model variant is ideal it is also infeasible. However, modellers should ensure that their approach is robust enough to deal with issues around existence and uniqueness. Checking that the model endpoint is reduced by (x%) when the (hbox {EP}_{x}) is applied to the exposure profile is an easy way to check accuracy and existence. An argument (if not a full, formal proof) for uniqueness should also be considered. In cases where that is not possible, the algorithm must be set up to identify the lowest (hbox {EP}_{x}), or check that no lower values exist.One common addition is to DEB-TKTD models which feature starvation is to assume that there is some maximum amount of starvation/shrinking which an animal can survive. Once that point is met or exceeded death is instantaneous27. Such death mechanisms cause problems. They can introduce a discontinuity in the response versus multiplier value for a given time window (i.e. a “jump” in plots such as Fig. 3). For instance, if in the example given in Fig. 3 the animal was not allowed to shrink, and instead died, then the multiplier of 6 would result in (100%) effects on survival (and significant growth effects). In contrast, the exposure when the multiplier is 2 is survivable and the animal can recover. Presumably, for some critical (alpha _c in (2, 6)) the exact threshold for death is reached. This (alpha _c) is a discontinuity between partial and (100%) effects relative to control. Under some circumstances this will prohibit finding a multiplier which results in exactly (x%) effects, regardless of the method used.There are two readily apparent solutions to this at the individual level. One is to set (alpha _c) as the multiplier for the window, the second is to replace such discrete responses with graded responses. In this example for instance, shrinking could add to the lethal hazard h. It is not possible to universally recommend one approach over the other as it will depend on the species’ behaviour. Once that decision has been made these issues must be recognised and reported by the modellers. More
