Toxicokinetic–Toxicodynamic Model to Assess Thermal Stress

Temperature is a crucial environmental factor affecting the distribution and performance of ectothermic organisms. This study introduces a new temperature damage model to interpret their thermal stress. Inspired by the ecotoxicological damage model in the General Unified Threshold model for Survival (GUTS) framework, the temperature damage model assumes that damage depends on the balance between temperature-dependent accumulation and constant repair. Mortality due to temperature stress is driven by the damage level exceeding a threshold. Model calibration showed a good agreement with the measured survival of Gammarus pulex exposed to different constant temperatures. Further, model simulations, including constant temperatures, daily temperature fluctuations, and heatwaves, demonstrated the model’s ability to predict temperature effects for various environmental scenarios. With this, the present study contributes to the mechanistic understanding of temperature as a single stressor while facilitating the incorporation of temperature as an additional stressor alongside chemicals in mechanistic multistressor effect models.


S01 Derivation of the model equations
Here we describe the development of the temperature damage module presented in this study, starting from the state-of-the-art damage module for chemicals as part of the widely used General Unified Threshold model of Survival (GUTS).The previously described GUTS approach (Ashauer et al. 2011;Jager and Ashauer 2018) simulates the probability of death of individuals over time-based on measured survival data.As proposed by Jager and Ashauer, the mechanistic approach of GUTS, though developed and increasingly used for chemical risk assessment, can be applied to other stressors that cause effects on survival (Jager and Ashauer 2018).
We perform a systematic analysis of the dimensions of the damage equation of the GUTS model (equation S1) and damage in our temperature approach (equation S7).The purpose of such analysis is to eliminate parameters that are not identifiable because they have the same effect on the model as another parameter.For instance if we want to do a regression of the following equation we cannot identify both the parameters a and b using data of the variables Y and X, because the model is overparameterized: Similarly, if variable Z is a latent variable that cannot be measured, we cannot determine parameter b in the following model because this parameter is dependent on the scale of Z: Y=a*X+b*Z+c In our approach, we cannot measure damage, so our purpose is to eliminate parameters that can only be determined if we would be able to quantify damage.

Derivation example based on the damage equation of the GUTS model
For the damage module in GUTS, a one-compartment model is assumed with first-order kinetics (eq.S1).
Here, the damage accrual is proportional to the external water concentration of the chemical (Cw), and the damage repair is proportional to the amount of damage (D).With those state variables and the damage accrual rate (ka) and the damage repair rate (kr), we can describe the damage dynamics for chemicals over time (t) with eq.S1. eq.S5 We then arrive at the following model (in which we dropped the *).
eq. S6 Because we chose   ̂= 1, this means that () should have the same dimensions as   ().Therefore, the damage as presented in the GUTs model is scaled to the chemical concentration and consequently, the unit of   ̂.

Derivation of the damage equation of the temperature damage model
The temperature damage model presented in this study was inspired the temperature injury model as presented by Jørgensen et al. (2021).With the aim to develop "a unifying model to estimate thermal tolerance limits in ectotherms across static, dynamic, and fluctuating exposures to thermal stress", their model enables the cross-study comparison of thermal tolerance measures.Their model is constructed of similar parts to the GUTS model (Table S1).However, it comes with its own limitations.The model of Jørgensen et al. (2021) only considers fixed temperatures or a linearly increasing temperature, and ignores damage repair.This considerably limits the use of this model to be applied on realistic temperature scenarios, varying from temperatures in-and outside of the temperature niche of the organism.We addressed these limitations in the temperature damage model by using a structure closely related to the damage equation of GUTS.In the temperature damage model, the damage accrual depends on the external temperature condition over time, expressed as some function g(T(t)).Thus our starting point is, to replace the chemical concentration as the stressor in eq.S1 with this temperature stress function (g(T(t)).Jørgensen et al. (2021) argue that temperatures' effect on survival is best described by an exponential relationship.There is namely much evidence that time to death (or knock-down time) increases exponentially with increasing temperature.This is first shown by (Bigelow 1921)

eq. S14
Note that this is completely analogous to the GUTS approach, but in this case it implies that   is dimensionless, because as an e-power is dimensionless.Which gives the used equation (in which we dropped the * and renamed kr to kT): ()  =   ( (()−  ) −   ()) eq.S15 This new temperature-related damage model describes the effect on survival induced by stressful (in contrast to modulating) temperature conditions.To estimate the model parameters based on the measured survival, the rest of the model is defined as in the main script (eq.1-3).

Symbol pairs
(GUTS ~ Jørgensen) show the parameter sets within the critical value (i.e., under the horizontal black line).The parameter symbols and units can be obtained from Table 1 in the main script.Solid lines show the model for the respective exposure scenarios (i.e.,10, 15, 20, 25 °C)  show the parameter sets within the critical value (i.e., under the horizontal black line).The parameter symbols and units can be obtained from Table 1 in the main script.

Function
Figure S1: Parameter space plot, calibration (for Tc = 11 °C).The plots on the diagonal show the profile likelihoods for the individual parameters.The scatter plots underneath are the 95% joint confidence regions.The yellow dots mark the best-fit values for each parameter and green dots

Figure S2 :
Figure S2: Predicted and observed survival probability for model calibration (for Tc = 11 °C).Vertical error bars represent the uncertainty of the model prediction (their 95% confidence intervals).Horizontal error bars are the Wilson scores intervals of the observed survival data.

Figure S4 :
Figure S4: Parameter space plot, calibration (for Tc = 14 °C).The plots on the diagonal show the profile likelihoods for the individual parameters.The scatter plots underneath are the 95% joint confidence regions.The yellow dots mark the best-fit values for each parameter and green dots

Figure S5 :
Figure S5: Predicted and observed survival probability for model calibration (for Tc = 14 °C).Vertical error bars represent the uncertainty of the model prediction (their 95% confidence intervals).Horizontal error bars are the Wilson scores intervals of the observed survival data.

Figure S6 :
Figure S6: Model simulations of damage and survival probability for different temperature scenario types.The different temperature scenario types are constant temperature scenarios (top row), daily temperature fluctuation (middle row), and daily temperature fluctuations with heatwaves (bottom row).Tc with 14 °C is marked with a dotted horizontal line.NOTE: While the survival probability is plotted for the whole simulation time (i.e., 200 days), the temperature and damage are plotted only for a representative period of the simulation (i.e., 2 and 30 days) for the daily temperature fluctuation and heatwave scenarios.Simulations were done with hb = 0.

Figure S7 :
Figure S7: Heatmap for the survival probability at the end of the simulation period (t=200 days) depending on different heatwave intensities and durations.During the heatwave the base water temperature (daily temperature fluctuations of 4 °C around the average of 12 °C) was increased by the intensity.For each scenario, the heatwave start was at day 10.Simulations were done with hb = 0 and Tc = 14°C.
(Schulte, Healy, and Fangue 2011;frequently confirmed for other species in the scientific literature literature(Schulte, Healy, and Fangue 2011; Sutcliffe, Carrick, and Willoughby 1981; Mundim et al. 2020; Jørgensen, Malte, and    Overgaard 2019; Aquilanti et al. 2010).Moreover, we know that within the temperature niche of the organism, there is no damage caused by temperature to be expected.Thus, we should define a threshold decided to use the exponential expression of Jørgensen et al. as the function to express the influence of temperature-related damage accrual: (Jager and Ashauer, 2018)d-Döring et al. 2022;Schulte, Healy, and Fangue 2011) temperature probably still modulates, as discussed elsewhere(Huang et al. 2023;Mangold-Döring et al. 2022;Schulte, Healy, and Fangue 2011).*̂−* ) −   ) eq. S13 Similar as the GUTS approach(Jager and Ashauer, 2018), we wish to keep normal units for time ( = 1) and temperature ( ̂= 1).With the following choice for   ̂, we remove one parameter (  ):

Table S2 :
Temperature damage model variable and parameter symbols and explanations.For model variables the respective model equation is provided, and for the parameters their best fit value alongwith their 95% confidence interval is provided.For this calibration, Tc was set to 14 °C.
*Boundary of the parameter space explorer