Dynamic Energy Budget (DEB) models are by now the most established models for predictive and interpretable modelling of sublethal chemical toxicity (Jager and Zimmer, 2012; Martin et al., 2013; Sherborne and Galic, 2020; Romoli et al., 2024a). These models describe the individual as a mass balance, starting with the amount of ingested food and ending with the central metabolic processes determining an individual’s life history: maintenance, growth, maturation and reproduction (Jager et al., 2013). DEB models can be extended to account for the idiosyncrasies of the amphibian life cycle (Romoli et al., 2025; Hansul et al., 2026) and effects of chemicals can be incorporated into DEB models by coupling them to a toxicokinetic-toxicodynamic (TKTD) submodel via physiological modes of action (PMoA). These translate the external chemical exposure to an effect on the energy budgets, e.g. a decrease in a conversion efficiency, increase in maintenance costs or change in resource allocation. A priori, there is no good reason to deviate from the standard approach for the application to amphibians. In general, it can be necessary to include other PMoAs to explain observed effects, or to consider the simultaneous action of multiple PMoAs (Koch and De Schamphelaere, 2021). While the toxicokinetic (TK) component can be expressed in different levels of detail, the most common approach in the context of DEB-TKTD models is a reduced TK submodule, where the external concentration is directly linked to an abstract quantity called scaled damage (Jager, 2020). The scaled damage is assumed to be proportional to a toxicologically relevant internal concentration, but does not require measured internal concentrations to estimate the corresponding rate constants, since the rate at which scaled damage is accumulated or repaired is primarily linked to the dynamics of effects after exposure or non-exposure. The coupling of such a reduced TK module to a DEB model further allows to include feedbacks with growth and reproduction. While these are often included by default, it is also reasonable to start with the simplest possible model, i.e. exclude all feedbacks and evaluate a posteriori whether a more complex model structure may be needed (Jager, 2020; Ockleford et al. 2018) Panel on Plant Protection Products and their Residues (PPR) et al., 2018). The application of DEB-TKTD models to the analysis of toxicity data has multiple advantages over a purely statistical approach (Jager, 2011, 2025). For example, it allows to test hypotheses regarding the mechanisms through which a substance may affect the organism (Koch and De Schamphelaere, 2021; Weighman et al., 2023). Based on such insights, the coupling to population models is possible (Jager and Klok, 2010; Martin et al., 2013) and such approaches have been used to predict differences between individual-level and population-level sensitivity (Pereira et al., 2019), as well as mixture effects on the population level (Vlaeminck et al., 2022). Further relevant applications of DEB and DEB-TKTD models include the modelling of interspecific sensitivity with regards to apical endpoints (Zubrod et al., 2024) and prediction of effects under time-varying exposure (Romoli et al., 2024b).

The objective of the present study was to test the applicability of DEB-TKTD models to analyze toxicity data in amphibians, specifically the anuran Discoglossus galganoi, in order to draw conclusions about mechanisms of toxicity and synthesize effects on different endpoints and on different life stages. We briefly describe and motivate the model structure, and report results on three case studies in D. galganoi, concerning the estimation of DEB parameters from control data with account for biological variability, the comparative effects of the herbicide 2,4-dichlorphenoxyacetic acid (2,4-D) and the insecticide Flupyradifurone.

The model was validated by predicting the effect of the binary mixture of 2,4-D and Flupyradifurone from the single-stressor effects. This involved two sets of predictions:

The latter prediction combines the prediction of mixture effects with the extrapolation of effects across life stages. To assess the model’s predictive capacity quantitatively, we calculated the mean absolute percent error (MAPE) for the selected PMoA:

y_i is the ith observed value and y^i is the ith simulated value.

For PMoA selection, additional quantitative metrics (MAPE and Nash-Sutcliffe Efficiency, Nash and Sutcliffe (1970)) are presented as Supporting Information, whereas the discussion in the main text is largerly based on visual comparison.

2,4-D caused delay of GS 46 and a decrease in the mass at GS 46. The effect of timing was predicted imprecisely, but accurately. The prediction accuracy for the timing of GS 46 was similar for PMoAs M and A (Figure 3). Effects on mass at GS 46 were under-predicted by both PMoAs, but less so by M. We therefore concluded that M is the more plausible PMoA to explain effects of 2,4-D. The associated MAPE (Equation 21) was 22% for the timing of GS 46 and 17% for mass at GS 46. The corresponding parameter values are reported in Table 3.

Effects of Flupyradifurone on the timing of GS 46 were predicted well by PMoA G (decrease in growth efficiency), but not by κ (shift in resource allocation from growth to maturation). While both PMoAs correctly predicted the effect on wet mass at GS 46, the difference between PMoAs was clearly visible in the predicted effects on the timing of GS 46. While we observed no clear effect of Flupyradifurone on the timing of GS 46, PMoA κ implied that GS 46 would be reached sooner at high exposure concentrations (a shortening of metamorphosis), which was not observed. PMoA G correctly predicted that the timing of GS 46 is unaffected by Flupyradifurone exposure (Figure 4). Therefore, we concluded that the effects of Flupyradifurone were more likely caused by a decrease in growth efficiency, rather than a change in the resource allocation from growth to maturation. The corresponding parameter values are reported in Table 4.

By fitting alternative PMoAs to effect data for larvae, we could draw some conclusions regarding the plausibility of PMoAs. For a given effect on growth, PMoAs M and A are associated with a clear delay of metamorphosis, whereas PMoA G only implies a subtle delay of metamorphosis. This can be explained by the interaction of these PMoAs with mass fluxes in the DEB model.

PMoA G, i.e. a decrease in growth efficiency, only causes an indirect effect on the maturation rate. Furthermore, any reduction in growth also reduces maintenance costs, and the effect of G approaches zero as an individual reaches its maximum length (the growth efficiency η_AS does not affect the equilibrium conditions of the model). PMoAs M and A on the other hand have a direct effect on growth and maturation. PMoA κ represents a special case, as it does not imply a loss of energy, but only a change in the allocation of energy. In the present study, κ did not explain observed effects. Generally, a decrease in κ might be a plausible candidate if premature metamorphosis is observed in combination with adverse effects on growth. The prediction of effects at GS 46 depended heavily on the choice of PMoA. The message here is two-fold: At the one hand, extrapolations from larval effect data to GS 46 should be done with caution. A model averaging approach may be useful if the PMoA identification from larval data is ambiguous and no further data can be obtained (Hansul et al., 2024). At the other hand, effect data at GS 46 apparently carries a lot of information regarding the PMoA. Generating such effect data is resource-intensive, but for cases where an unequivocal PMoA identification is crucial and additional data generation can be ethically justified, it could be worth the effort.

Selection of PMoAs here relied heavily on visual evaluation. This choice of a primarily visual evaluation was also attributable to the fact that easily applicable quantitative metrics like MAPE compare the average goodness-of-fit, but in no way capture how well a model describes the temporal dynamics. If the difference in PMoAs clearly manifests in the presence/absence of effects on certain endpoints, this may be negligible. If the difference in PMoAs is more nuanced and related to the biological plausability of the produced temporal dynamics, assuming independence for time-series data becomes more problematic. While more objective measures for PMoA selection are desirable, reliance on such metrics may require that the sequential ordering of time-series data is taken into account. Interesting developments in statistics that may help approach this issue in the future include, for example, the use of path signatures as summary statistics Dyer et al. (2023). Uncertainties of the inferred parameters were notably higher for k_D (dominant rate constant) and b (Hill’s slope) than for e (median effective damage/sensitivity). It is known that the identifiability of k_D can be an issue (Billoir et al., 2008). This issue is to be expected when the data implies fast damage dynamics, as observed here for 2,4-D due to the immediate onset of effects on growth. Limited identifiability of b can be explained by the fact that in both datasets, the relevant effects occured in a single test concentration, which differed by one order of magnitude from the next-lowest test concentration. Consequently, the amount of information in the data that can be used to derive a slope was limited.

The model’s tendency to under-predict mass at GS 46 cannot be attributed to a single explanation, based on the data analyzed here. Based on the structure of the baseline model, possible modifications that could explain a lower relative mass change during metamorphic climax are differences in size-specific maintenance costs during larval development and climax, as well as a different energy density of the biomass consumed during climax, compared to the remaining biomass. While the baseline DEB model could overall capture the effects of chemical stressors across life stages, we suggest that future refinements should focus on mass fluxes during climax (further discussed in Hansul et al. (2026)).

The dominance effect of 2,4-D over Flupyradifurone could be predicted with high accuracy. Given that the substances acted via different PMoAs, a differentiation between IA and DA models was obsolete. The presence of a dominance effect can be explained through the combination of PMoAs, based on the same rationale that was used to select the PMoAs in the first place. We identified G as most likely PMoA for Flupyradifurone, implying that effects weaken over time and only indirect effects on maturation. We identified M as most likely PMoA for 2,4-D, implying effects across the entire growth curve, as well as direct effects on the maturation rate.

While the present study focused on the modelling of sublethal effects on the individual level, one of the goals of the amphibian DEB-TKTD model development is to lay out the basis to extrapolate effects to the population level. Achieving this would at least require that the present analysis is extended to also account for lethal effects. In fact, the severity of lethal effects is highly life stage-specific (González-López et al., 2026), with the highest mortality rate occurring during metamorphosis, manifesting metamorphosis as a bottleneck in the amphibian life history.

Mortality during climax at the one hand reduced sample size for endpoints at GS 46, somewhat weakening the conclusions that can be drawn regarding predictive performance for GS 46 endpoints. At the other hand, this highlights the necessity to include lethal effects for a complete description of effects on individual life history and possible integration into population modelling approaches.

Apart from direct lethal effects, the ecology of amphibians implies that sublethal effects in laboratory tests may translate to lethal effects in the field.

In ephemeral aquatic habitats such as those inhabited by D. galganoi, a delay in metamorphosis could lead to death if individuals fail to reach metamorphosis before the aquatic habitat has dried out. This means that the results presented in the present study could lead to dramatically different effects on population dynamics if the effect of drying is considered. Simultaneously, it is known that amphibians may accelerate development under drying conditions, counter-acting the effects described here. Such acceleration of development may in turn lead to adverse effects in later life stages. The matter may be further complicated by environmental temperatures, which simultaneously affect aquatic habitat quality as well as metabolic rates. These are relationships between chemical exposure and amphibian life history that have not been considered in the present study, but could be integrated into a DEB-TKTD modelling approach.

We have demonstrated that DEB-TKTD models can be used to perform predictive modelling of pesticide effects in amphibian larvae and metamorphs. While this study was limited to a single binary mixture combination with substances that showed different PMoAs, validation studies with substances that share PMoAs would be extremely helpful to further test the robustness of such predictions. Furthermore, we tested the approach for a single species. While the model structure did not contain species-specific assumptions, validation of the approach with other species, ideally those exhibiting contrasting life-histories, is needed for more generalized claims regarding the model’s predictive power. The extrapolation of effects across life stages depends heavily on the selection of the PMoA, which currently relies heavily on visual evaluation. Quantitative metrics may support the selection, but the fact that these typically assume independence of observations in time-series data deserves further attention. DEB-TKTD models can be fitted to larval data up to GS 42 and some inference about the PMoA can be made, but data on effects at GS 46 is highly informative regarding the PMoA.

Generally, the model has a tendency to under-predict effects on the mass at metamorphosis. The effects of a binary mixture of pesticides with different PMoAs could be predicted with high accuracy.