The application and limitations of exposure multiplication factors in sublethal effect modelling

Thanks to growing interest and research in the field, toxicokinetic–toxicodynamic (TKTD) models are close to realising their potential in environmental risk assessment (ERA) of chemicals such as plant protection products. A fundamental application is to find a multiplicative scale factor which—when applied to an exposure profile—results in some specified effect relative to a control. The approach is similar to applying assessment factors to experimental results, common in regulatory frameworks. It also relies on the same core assumption: that increasing the scaling always produces more extreme effects. Unlike experimental approaches, TKTD models offer an opportunity to interrogate this assumption in a mathematically rigorous manner. For four well-known TKTD models we seek to prove that the approach guarantees a unique scale factor for any percentage effect. Somewhat surprisingly, certain model configurations may have multiple scale factors which result in the same percentage effect. These cases require a more cautious regulatory approach and generate open biological and mathematical questions. We provide examples of the violations and suggest how to deal with them. Mathematical proofs provide the strongest possible backing for TKTD modelling approaches in ERA, since the applicability of the models can be determined exactly.

Assessing the ecological consequences of introducing chemical compounds, such as plant protection products (PPPs), into an environment requires some form of environmental risk assessment (ERA). At the individual level, an assessment compares the effects that exposure to the substance has on the organism to the level of exposure that it is likely to experience in the field. Traditionally, laboratory studies on standard species under relatively constant exposure have been used to derive certain summary statistics 1 . For instance, a no observed effect concentration (NOEC) or an estimated concentration where x% effects relative to the control performance occur ( EC x ). These statistical results are then compared to a predicted peak or time-weighted exposure concentration to derive a toxicity exposure ratio (TER). The final step is to compare the TER to some assessment factor (AF). The AF is designed to cover uncertainties relating to extrapolation to different species, conditions and exposure. If the TER is greater than the AF then the use case of the PPP is considered low risk. Setting aside the known methodological flaws in NOECs and EC x values 2 , the approach cannot account for exposures which fluctuate over time. Variable exposure is typical in the field and can significantly alter the effects on the individual. Data collection and modelling methods exist which generate time series of concentrations present in the environment over time. For instance, in the aquatic environment, surface water models predict concentration levels over time in water bodies near treated fields. We refer to these time series concentrations as exposure profiles, since they determine the variable level of concentration organisms will be exposed to over time. ERA methods which accurately incorporate time-variable exposure are therefore valuable. However, it is impossible to test in the laboratory every exposure profile which may occur in the field.
In recent years, there has been growing interest in mechanistic models which can explain the effects of exposure to a harmful substance on an individual, population or community of organisms [3][4][5][6][7] . At the individual level, toxicokinetic-toxicodynamic (TKTD) models describe how given exposure conditions ultimately translate to effects on the organism over time 8 . Due to their mechanistic structure it is possible to extrapolate to predict effects under a different set of conditions, in particular to untested, time-variable exposure profiles. This ability to take into account the whole exposure profile is a great strength of TKTD models.
Within ERA, the pioneer of these models is the General Unified Threshold model of Survival (GUTS) 4 , a TKTD model for lethal effects which has been approved by the European Food Safety Authority (EFSA) for use

Materials and methods
The exposure multiplication factor approach. It is possible for TKTD models to incorporate many different aspects to increase realism, including modelling different temperature 14 and food availability, and mixtures of stressors 15 . However, there is currently no regulatory guidance on how to deal with the interaction of all of these sources of variability. Therefore, they are not currently included in TKTD models for standard ERA applications 7 . Instead, with the exception of the exposure concentration, all conditions in the model are typically fixed to the laboratory conditions maintained during the empirical bioassays. In essence, the model is simulating in silico versions of standard bioassays, but with many more exposure profiles than would be feasible in empirical, laboratory based experiments. We therefore refer to this strategy as a laboratory mimic approach. The approach limits complexity by focusing on a more thorough consideration of exposure compared to methods which produce single summary statistics.
Once a suitable TKTD model has been calibrated it is possible to run the model with any exposure profile and see the effects on the model organism relative to control conditions. For implementations within ERA, the model should also be properly validated according to the relevant guidelines 7 . An EMF scales the whole exposure profile according to that value. An EMF of 10 means that at any time the exposure in the model is 10 times greater than what is actually specified at that time in the unmodified exposure profile. The goal is to find an EMF value which imposes some specified percentage reduction in an aspect of organism performance. Typically, these will be the same as the measured outcomes in laboratory tests, for example survival probability or reductions in growth or fecundity.
The accepted use of GUTS in ERA finds the LP x value based on the survival probability at the end of the exposure profile. That is, LP 50 is the EMF which causes 50% toxicant induced mortality at the end of the exposure. We will extend the approach to the DEB-TKTD, algae and Lemna models, as is presumed by EFSA 7 . For sublethal effects, such as growth, reproduction or biomass, the EP 10 is the EMF which causes a 10% reduction in that quantity compared to control conditions in the model.
To find the LP x (or EP x ) an iterative algorithm tests different EMF values until exactly (within some numerical tolerance) x% effects occur at the end of the exposure profile. That EMF is the LP x . For efficiency, a root-finding algorithm is often used. Over many iterations the algorithm should converge to the LP x . However, if no LP x value exists, the algorithm could either run indefinitely or report the last tested EMF value before it stopped. In In order to construct the proofs, we introduce new notation to signify a state variable which may change in response to the EMF. This is best done with an example. Consider survival, denoted by the state variable S(t). This is obviously dependent on time t but also, perhaps implicitly, dependent on the EMF, which we will denote α . Survival when exposure has (potentially) been altered by the EMF α is denoted by S(t; α) . As a result, survival under control conditions (i.e. zero exposure) is equivalent to S(t; 0) for any exposure profile, whereas the response to the unmodified exposure profile, i.e. prior to any multiplicative changes, is equivalent to S(t; 1).
Proof methodology. For each model, the critical endpoint will be a change in one or more state variables in the model relative to the same variable under control conditions. Therefore we must show that a higher multiplier value applied to the exposure profile will always cause greater or equal effects on these state variables, a property known as monotonicity. Secondly, each variable affected by the EMF should be continuous with respect to the EMF, such that there are no values of the EMF which introduce discontinuities or undefined values in the state variables. Then, for some desired x% response, all that remains is to show that arbitrarily large EMFs cause greater than x% effects. When all these conditions are met, by the Intermediate Value Theorem (IVT) we know that some critical EMF leading to exactly x% exists. For uniqueness we also require strict monotonicity (without equality) in some neighbourhood of the critical EMF, otherwise the critical EMF will be some interval, rather than a single value.
Each of these proofs result in one or more theorems for the corresponding TKTD model. Only the theorems are presented in the main text. The details of each proof, and all model equations, are provided in the SI.

Models.
The GUTS framework. The GUTS framework is a set of nested models which quantify organism survival over time, typically under chemical stress 4 . In the full version of the model, uptake and elimination of the external concentration occurs according to a simple one compartment ODE which describes the internal concentration over time. Internal concentration then induces damage, another one-compartment model with accrual proportional to the internal concentration and repair proportional to the current damage. Each organism has some threshold for sensitivity drawn from a probability distribution, typically the log-logistic function. The individual's probability of death increases linearly according to some rate parameter (the killing rate) multiplied by the extent to which damage exceeds the threshold.
Typically, reduced forms of the model (referred to as GUTS-RED) are expected to be used in the ERA of PPPs 7,16 . These reduced models replace the internal concentration and damage ODEs with a single state variable, scaled damage, so called because it has the same dimensions as the external concentration 4 . As a further simplification, the death mechanism is simplified to one of two extreme cases. The first of these, stochastic death (GUTS-RED-SD), assumes that each organism has the same threshold value, once the scaled damage exceeds the threshold, death is a chance process. Secondly, the individual tolerance (GUTS-RED-IT) model takes the killing rate to infinity, such that death is instantaneous for the individual as soon as its threshold damage is exceeded. These simplifications enable GUTS-RED-SD and GUTS-RED-IT to be calibrated from standard mortality bioassays. A full schematic of the model, similar to that presented in 16 is presented in Fig. 1.
The output of the model is the survival probability over time. The goal is therefore to find the EMF value which causes an x% reduction in the survival probability due to toxicant induced mortality at the end of the exposure profile. www.nature.com/scientificreports/ DEB models. DEB also does not refer to a single model but rather a set of models which share many core assumptions. DEB theory uses a number of physical and biological laws to derive a mechanistic model of ODEs for the growth, reproduction and survival of an organism 17 . DEB-TKTD (sometimes referred to as DEBtox) models include additional equations to describe the uptake, transformation and elimination of a substance and the stress it causes. As such it can model the endpoints from standard aquatic animal bioassays and make predictions about untested exposures. The full suite of DEB-TKTD models has recently been reviewed for their suitability for standardised ERA applications in the near future 18 . For this analysis we will use the recent reserve-less DEBkiss model of Jager 11 . We chose this model as it is a simple, comprehensive implementation for an animal's post-birth lifecycle under some known exposure profile C w . This simplicity reduces the amount of data required for model calibration.
The derivation of the model is presented thoroughly in the original publication. A summary of the model and all state variables and model parameters are given in the SI. Here only a brief description is given.
All DEB models consider the assimilation, distribution and dissipation of energy within a single organism. In DEBkiss models, energy assimilated through food is immediately split between two branches. The somatic branch receives the fraction κ of all assimilated energy. From this energy, somatic maintenance (proportional to the volume of the organism) is paid, and the remainder goes toward growth of new structure (biomass). This structure is tracked with the state variable structural length, L. The remaining (1 − κ) fraction is spent on maturation processes. Energy in this branch is first used to pay maturity maintenance (again linked to volume in this model) with the remainder being invested in maturation (for juveniles) or reproduction (adults). Here, the model of Jager 11 makes two common simplifications to reduce model complexity. Firstly, we assume that size is a proxy for maturity, this means that maturation does not need to be modelled; the switch from juvenile to adult occurs once the animal has reached a given size. Secondly, the model assumes that reproduction is a continuous process rather than discrete. Although untrue, this assumption is useful for many species 19 . Cumulative reproduction R c is the second state variable.
When conditions worsen, for example due to increased exposure, it is possible that the organism will no longer have enough energy to pay somatic maintenance. This is starvation. Different species are likely to have different responses when this occurs. For our purposes, the model animal ceases growth and redirects sufficient energy from the reproductive branch to the somatic branch to pay somatic maintenance. If somatic maintenance costs exceed all energy assimilated from food then the animal must also shrink. It does this by burning structure to release energy with yield y P , between zero and one. Structural length, L, decreases during this phase.
While different core DEB model variants may be used for numerous reasons, e.g. to better suit the model organism, substance or model environment, the TKTD elements differ little between these 18 . The best approach to utilising DEB-TKTD in ERA is therefore also not species specific. Traditionally, stress in DEB-TKTD models was directly linked to the (scaled) internal concentration 17,20 . However, experience with GUTS has shown that often toxic effects are not well described by the internal concentration 21 , and thus a scaled damage state variable is now proposed for DEB as well. The uptake and elimination rate of damage is governed by a dominant rate constant, k d . Growth and reproduction processes can alter the uptake and elimination of damage, and thus may have an effect on damage dynamics. These influences on damage are known as "feedbacks" 11 . All feedback processes are represented as x i , where the subscript changes to describe each process. If uptake or elimination occur across a membrane (e.g. the skin) which grows as the organism grows, then the surface area to volume ratio will affect the rate of either process. This ratio changes over time. These feedbacks are denoted x u and x e respectively. Additionally, damage can also be diluted, either by growth ( x G ) or by reproduction ( x R ). Each process may or may not be relevant for any given stressor, this is represented by a binary vector Under this scheme, traditional DEB-TKTD models which used internal concentration, have X = [1, 1, 1, 0] . Suggestions on when each feedback is likely to be active ( X i = 1 ) or inactive ( X i = 0 ) are given by Jager 11 .
Above some tolerance threshold, damage creates stress, s, to the organism, which increases linearly, just as in the GUTS framework. This stress then perturbs one or more processes according to the physiological mode of action (pMoA) of the substance. Standard pMoAs are an inhibition to assimilation, s A , increased maintenance costs, s M , increased growth costs, s G , increased reproduction costs, s R and hazard to the embryo during oogenesis, s H . Stresses not active for a given pMoA are zero. We will sometimes refer to the pMoA as a binary vector S such that It is common to couple a DEB-TKTD model to a survival module following the GUTS-RED-SD framework 22,23 . We will also follow this approach to extend the utility of the results. The survival probability state variable is denoted by S. A full schematic of all elements of the model is given in Fig. 2.
Comparing the schematics in Figs. 1 and 2 shows the similarities in many elements of the model structures. For both models, the core elements of damage accumulation and repair are first-order processes. Indeed, if there are no feedback processes ( X = [0, 0, 0, 0] ) the survival probability over time predicted by the DEB-TKTD model will be identical to a GUTS-RED-SD model.
The extra complexity of the DEB model means that the model organism is no longer constant over time, it grows and reproduces. This introduces new state variables for these processes which can also be affected by stress, and investigated with the EMF approach. Additionally, the growth and reproduction processes can impact the damage uptake and elimination processes through the feedback processes. These processes are not included in a GUTS model, since there the model organism is not changing over time. For more information on GUTS and DEB models, see EFSA 2018 7 . (1)

Results
The GUTS model. For 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; α) denote the survival probability at time t for a given non-zero exposure profile C w (t) scaled by some EMF value α . 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 α * such that α * is the LP x for the exposure profile.
Baudrot and Charles 10 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 LP x is monotonically increasing with respect to x. For example, the LP 50 will always larger than the LP 10 for the same exposure profile. This result comes directly from (S4) in the SI.  Then, for any combination of feedbacks X = [X u , 0, X G , X R ] , damage is monotonically increasing with respect to α , and is continuous with respect to α as long as changes to L and R are continuous. Moreover, damage is strictly monotonically increasing with respect to α whenever D(t; 1) > 0.
The limitation that 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. Finally, we can step from the results of Theorem 2 and Corollary 3.1 to show the existence and uniqueness of a critical multiplier ( EP x or LP x ) for growth, reproduction and survival. 11 and a substance such that at least one of Theorem 2 or Corollary 3.1 hold. 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 α * > 0 such that

Theorem 3 Consider the DEB-TKTD model of Jager
The monotonicity of effects on all state variables in the DEB model means that, for the conditions described in Theorem 3, the EP x (or 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 equation 24 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 One can now find the births per individual per time predicted by the DEB model as 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 t E 0 B(t)dt for each EMF and compare it to the control, similar to finding EP x values for individual endpoints. Alternatively, one can use B(t) as the basis to estimate the intrinsic population growth rate 25 . 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 tests 26 . For the first of these options we offer an extension to Theorem 3.     11 . In each case, t 1 is first time point where the external concentration is non-zero. Each entry in the table shows the condition which must be met in order for x% effects to be possible at the end of the exposure profile ( t E ). N/A is used for pMoAs which do not affect the state variable. In cases where the condition is not met for a given variable, a multiplier may still be found for the others. Under assimilation stress (fourth column) there is a maximum rate at which length L can decrease. This rate is the von Bertalanffy growth rate, r B , divided by the yield of burning structure to provide energy, y P . *At least one of the two reproductive pMoAs must be active.

Endpoint of concern
S , pMoA  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 EP x value.
Surface:volume scaling of elimination. There is a reason that in Theorem 2, X e = 0 was specified. In some cases when 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. X = [0, 1, 0, 0] ). The damage ODE under some EMF α is then 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, α 2 > α 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 α . 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 α 2 C w is smaller and therefore able to eliminate damage more rapidly because X e = 1 . This eventually leads to lower damage for the model organism exposed to α 2 C w (i.e. D(t; α 1 ) > D(t; α 2 ) ). The more is less phenomenon can also impact growth and cumulative reproduction, as seen in Fig. 3b,c. Sometime after exposure L(t; α 2 ) > L(t; α 1 ) and R(t; α 2 ) > R(t; α 1 ) . For survival, and any additional endpoints without recovery, this "crossover" is unlikely, mortality during the exposure phase (where D(t; α 2 ) > D(t; α 1 ) ) will almost certainly dominate any mortality during the recovery phase. Figure 3d shows that for certain x% effect levels (vertical axis) multiple EP x values exist.
In 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 method 18 consisting of many exposures constructed sequentially and assessed. Each window will produce an 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) EP x . However, the potential for multiple 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 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 X e = 1 is necessary in model calibration and validation. Under such circumstances algorithms must ensure that the 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 Jager 11 . For model configurations where non-uniqueness could emerge using another method to find the 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.  www.nature.com/scientificreports/ Other issues. The 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 ≥ 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 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 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 EP x , or check that no lower values exist.

Scientific Reports
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 instantaneous 27 . 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 α c ∈ (2, 6) the exact threshold for death is reached. This α 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 α 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.

Discussion
Using mechanistic models to better understand and predict the effects of exposure to potentially harmful substances such as PPPs is clearly of great value. Especially for aquatic organisms, where we often have the luxury of a well defined predicted exposure, the path ahead may seem obvious. However, as we have shown, implementing models-especially for sublethal effects-requires care.
In most previous work, it has simply been assumed that the EMF approach will produce a unique value for any level of effects relative to the control. This is not the case. By pursuing rigorous mathematical proofs we have found certain restrictions and exceptions. For the GUTS model framework, any non-zero exposure profile has a unique LP x for any x. For the primary producer models, the major findings are the restrictions on the exposure profile. Since both the algae and Lemna models have some maximum rate of biomass or population decay, the first instance of non-zero external concentration in the exposure profile must occur with sufficient time before the end of the profile for the effect of the substance to exceed x% . For more information see the SI. Similar findings occur for some of the pMoAs in the DEB model, as shown in Table 1. The reason behind this is the same as for the primary producer models.
When using these models to make predictions it is important to realise that we must protect organisms (or populations in the case of algae) of different ages and sizes. For GUTS models this has no effect, but for DEB the time that exposure occurs within its lifespan is important. For instance, effects of exposure to a substance which increases growth costs will have more effect on a juvenile, growing organism than a fully grown adult. For this reason, moving time-windows have been suggested in ERA. The method constructs consecutive intervals of time from the exposure profile of concern and finds the EP x for each window. The EP x for the profile as a whole is the minimum of these values 18 . The approach also guarantees that some windows will have immediate exposure. This means that, although non-existence could occur for certain windows, a relevant EP x for the exposure profile will always exist. These results provide confidence in the EMF method and the use of the most efficient root-finding algorithms to calculate LP x and EP x values for these models.
However, for DEB models there are further limitations due to the possible feedbacks and model extensions. The feedback options in (1) are important to accurately describe damage dynamics in a one-compartment module. Finding the appropriate feedback mechanisms can significantly improve the correspondence between the model and the data 11 . Until recently the configuration X = [1, 1, 1, 0] was the default in DEB-TKTD applications 17,22 . This meant that every model application featured the problematic surface area:volume scaling of elimination. Although these feedbacks match our intuition and knowledge for dermal uptake and elimination the "more is less" multiplier paradox is a non-intuitive result which can arise in this scenario and indeed any other scenario with X e = 1 . However, as was shown in Theorem 3 and Corollary 3.1, uniqueness is guaranteed for all reproductive pMoAs, and substances with any pMoA as long as X e = 0.
Whether the behaviour illustrated in Fig. 3 can occur in real-world tests is unknown. Experimental data exploring this would be of value regardless of the results. The majority of ecotoxicology, especially ERA, is built on the presumed idea that scaling (i.e. multiplying) exposure will scale the effects in the same direction. Using the well-known and highly respected DEB theory we have shown that this might not always be the case. If the www.nature.com/scientificreports/ behaviour in the model were to be observed experimentally it could change our perspective on ranking different exposures in terms of effects, at least at the individual level. On the other hand, if it was clearly refuted it could signal a need to alter and improve the TKTD aspects of the model. The issues raised in this work are therefore important. It is imperative that modellers and anyone using the models are aware of the domain of applicability and demonstrate that they are within them. This applies to all TKTD models. The simplest solution to non-uniqueness in a regulatory context is to avoid root-finding algorithms in favour of (for example) a robust brute-force approach checking multiplier values from zero and increasing in small increments (e.g. 0.1) to approximate the lowest EP x . Accuracy and conservatism are thus maintained by sacrificing elegance and computational efficiency.
The EMF approach of replicating standard in vivo experiments in silico is appealing. Models are used to generate the same endpoints that are recorded at the end of real-world experiments specifically designed and performed for the ERA of PPPs. However, taking proportional effects at the end of the in silico experiment is only one option. We have shown in Corollary 3.2 one example of how a new endpoint can be derived from model state variables and used for ERA. Alternatively, one could define the endpoint for ERA as the maximum proportional effect at any point in the window. Although we cannot a priori determine at what time the maximum effects will be observed, the results of Theorems 1 and 3 as well as Theorem S1 and Theorem S2 in the SI hold for any value of t and so will also guarantee uniqueness and existence of α * as long as the time at which the difference is assessed remains constant for all EMFs. However, this approach is generally not recommended. To give one example, if reproduction is delayed at all compared to the control, during that delay there is a 100% reduction in reproduction compared to the control. This work has capitalised on a hitherto under-utilised aspect of TKTD models: the ability to provide rigorous proofs of certain properties. By doing so, we have conclusively shown under which circumstances LP x and EP x values exist and are unique. These model results also lend some verification to the use of multiplicative AFs applied to standard bioassay results. It may not always be possible to prove similar results for other TKTD models. Nonetheless, it is essential to consider all possible results and behaviours of a model before using it in any official context.