Predictive modelling of chemical mixture toxicity to amphibians in Discoglossus galganoi larvae and metamorphs

1 Introduction

Amphibians are currently not explicitly considered in the environmental risk assessment (ERA) of chemicals. It is rather assumed that assessments based on aquatic vertebrates (e.g. fish) and terrestrial vertebrates (birds and mammals) are also protective for amphibians. While sensitivity of amphibian larvae to chemicals is often well correlated to that of fish (Weltje et al., 2013), this does not apply to all chemicals (Glaberman et al., 2019), especially following long-term exposures. and the process of metamorphosis represents a vulnerable point in the life history of amphibians, which is not covered by using other vertebrates as proxies.

Furthermore, environmental exposure occurs in mixtures. The relative effect of mixtures within the metamorphosis phase is so far largely unexplored, and it cannot be expected that mixture toxicity data will be generated for a wide range of species and chemicals, necessitating the development of predictive models.

In this study, we demonstrate the possible use of toxicokinetic-toxicodynamic (TKTD) models to support a reliable ERA of chemicals for amphibians, including the extrapolation of chemical effects from larvae to metamorphs and the prediction of mixture effects. In general, this approach aims to maximize the amount of knowledge gained from experimental data, and to minimize the need for future animal testing.

1.1 Setting the scene: predictive modelling of chemical toxicity with Dynamic Energy Budgets

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).

1.2 Objective

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.

2 Methods

2.1 Model description

2.1.1 Amphibian DEB model

The DEB-TKTD model was derived from DEBkiss (Jager et al., 2013) and has been extended to account for amphibian metamorphosis (Hansul et al., 2026). In short, the model includes a metamorphic reserve compartment E^mt, which is defined as the amount of biomass that is used during metamorphosis to fuel maintenance costs and maturation. The parameter γ determines the allocation of assimilated resources into the metamorphic reserve compartment. Under ad libitum feeding conditions and at equilibrium, γ is equal to the fraction of biomass that is reserve.

During metamorphosis, feeding rates decline and are assumed to reach 0 with the completion of metamorphosis. The decline in the size-specific feeding rate is coupled to the depletion of metamorphic reserves:

(1) {I˙}maxmt={I˙}maxlrvEmtEmaxmt

Where {I˙}maxlrv is the maximum size-specific ingestion rate for larvae, {I˙}maxmt is the current maximum size-specific ingestion rate for metamorphs and Emaxmt is the internally tracked maximum reserve level, reached at metamorphic climax (Equation 1).

By definition, metamorphosis is completed when the metamorphic reserve is depleted. The full dynamics of the model (Hansul et al., 2026) are given by Equations 2–11:

(2) X˙emb={−I˙if embryo0otherwise

(3) X˙={0if embryoI˙otherwise

(4) I˙={0if embryo{I˙}max f([X])S2/3otherwise

(5) A˙=yAηIAI˙

(6) M˙=yMkMS

(7) J˙=yMkJH

(8) S˙={yGηAS(1−γ)(κA˙−kMS)if larvayGηASA˙if metamorphyGηAS(κA˙−M˙)otherwise

(9) E˙mt={yGηASγ(κA˙−kMS)if larva−(M˙+J˙+H˙)if metamorph0otherwise

(10) H˙={0if adult(1−κ)A˙−J˙otherwise

(11) R˙={yRηAR[(1−κ)A˙−J˙]if adult0otherwise

with the following rules for life stage transitions (Equations 12–16):

(12) embryo={trueif Xintemb>0falseotherwise

(13) larva={trueif Xintemb=0∧H<Hjfalseotherwise

(14) metamorph={trueif H≥Hj∧Emt>0falseotherwise

(15) juvenile={trueif H≥Hj∧Emt=0∧H<Hpfalseotherwise

(16) adult={trueif H≥Hpfalseotherwise

Symbols for the parameters and state variables of the baseline model are listed in Table 1. The values y_G, y_M, y_Aand y_κare intermediate quantities which account for the effects of stressors. Additionally,

Table 1. Parameters and state variables of the baseline DEB model.

the zoom factor Z, calculated as the ratio between the maximum structural masses of two organisms, is used to induce variability between individuals (Hansul et al., 2026). This leads to variability in growth rates, maximum size, time to reach metamorphosis, time to complete metamorphosis, and body mass at the beginning and end of metamorphosis, respectively.

DEB parameters, including σ_Z, have previously been inferred for Discoglossus galganoi. Baseline parameter values used in this study are taken from Hansul et al. (2026) and reported in Table 2.

Table 2. Parameter of the physiological baseline model estimated from the control data.

Best fit was the value associated with the smallest error in the calibration and used for simulation. Median refers to the median of the posterior distribution, P₀₅ and P₉₅ refer to the 5th to 95th percentiles (credible intervals). For further details, see Hansul et al. (2026).

σ_Z,M is the standard deviation of the zoom factor in relation to the maximum structural mass of the average individual.

2.1.2 Toxicokinetic-toxicodynamic model

To model the effects of chemicals, we first applied a minimal component for the damage dynamics:

(17) D˙z,j=kD(CW,z−Dz,j)

Equation 17 links the aqueous concentration C_W,zof stressor z to the scaled damage D_z,j, which is specific to the physiological mode of action (PMoA) j. This formulation ignores possible feedbacks between scaled damage and body size (Jager, 2020), in an attempt to first deploy the simplest possible model that can explain and predict observed effects.

For each possible PMoA j, the corresponding damage D_z,jis translated to a relative response y_z,j. We assume a log-logistic relationship between damage and response, corresponding to the standard assumption used in ecotoxicological analysis and resulting in a continuously differentiable TKTD model. For PMoAs with a monotonically decreasing relationship, the response can be directly calculated from the survival function of the log-logistic distribution:

(18) yz,j=11+(Dz,jez,j)βz,j

This function has two parameters, the sensitivity e_z,j(mg/L) and the slope β_z,j(dimensionless). In the present study, PMoAs with decreasing response are: Decrease in growth efficiency (G), decrease in assimilation efficiency (A) and decreased investment in growth/increased investment in maturation (κ). For these PMoAs, e is interpretable as a median effective damage, i.e. the damage level for which the corresponding flux is decreased by 50%.

For PMoAs with a monotonically increasing relationship, we use an affine transformation of the cumulative hazard function of the log-logistic distribution:

(19) yz,j=1−ln[11+(Dz,jez,j)βz,j]

In the present study, the only PMoA with an increasing relationship is an increase in maintenance costs (M). The parameter e_z,jnow corresponds to the damage level at which a 1.7-fold increase in maintenance costs occurs. Equations 18, 19 are defined so that the resulting relative response is a factor that can be inserted directly into the DEB equations. The intermediate variable stress (Jager, 2020) can be skipped without consequences for the model dynamics.

For mixtures of stressors with different PMoAs, their effects can be applied independently. For stressors with different PMoAs, we assume that an independent action (IA) model applies, where the combined relative response y_jis the product of the individual responses:

(20) yj=∏z=1nyz,j

Equation 20 yields the values y_G, y_M, y_Aand y_κused in Equation 11. The alternative to IA is a damage addition (DA) model. In the context of DEB-TKTD models, DA requires the estimation of an additional weight parameter (Jager et al., 2010). This means that the data is either fitted to the mixture data - which would contradict the goals of the present study - or that all TKTD parameters have to be estimated simultaneously from the combined single-stressor data, as previously done for lethal effect modelling (Bart et al., 2021). The assumption of IA is therefore a plausible and practical starting point. The model was implemented in the Julia programming language. All code and data is available on Github (see Supplementary Data).

2.1.3 Parameter inference

TKTD parameters were inferred from single-stressor effect data of 2,4-D and Flupyradifurone (González-López et al., 2026). For details on data generation and statistical analysis, see González-López et al. (2026). Chemical stressors were applied in their commercial formulations Primma ^®Dos and Sivanto ^®Prime, respectively.

Based on the available data, we cannot make statements about the possible contributions of other constituents than the active substance, but refer to the stressors as 2,4-D and Flupyradifurone for consistency with the original study. This dataset had been generated with considerations regarding mechanistic model development in mind, and therefore is likely the most suitable dataset currently available with regards to validating a mixture DEB-TKTD approach. González-López et al. (2026) investigated toxicity to Discoglossus galganoi and Pelophylax perezi, but since the observed sensitivity of P. perezi was much lower and effects were less consistent, the dataset for D. galganoi was deemed more relevant for the present study. For further information on experimental aspects, see González-López et al. (2026).

The calibration data consisted of time-resolved larval wet mass, the time to reach Gosner stage (GS, Gosner (1960)) 42 and the wet mass at GS 42. The time to reach GS 42 was used as the time-resolved fraction of tadpoles, calculated as the number of survivors which had not yet reached GS 42, relative to the total number of survivors at the given time point. A slower decrease in this value over time therefore indicates delayed metamorphosis, unrelated to survival. The highest tested 2,4-D treatment (100 mg/L) was omitted from the calibration, because it caused 100% mortality quickly after the onset of exposure. The presented analysis exclusively deals with sublethal effects.

For each stressor, we pre-selected plausible PMoAs and separately fitted a TKTD model based on each PMoA. The plausibility of PMoAs was evaluated by visual comparison. All parameter inferences were done using population Monte Carlo approximate Bayesian computation (PMC-ABC (Beaumont et al., 2009)), using the Euclidean distance as distance function, which is a default choice in ABC (Prangle, 2017). PMC-ABC is a practical choice because the model output was stochastic, and PMC-ABC is well-suited to deal with stochasticity in model outputs. The exact configurations of PMC-ABC runs are given in the calibration script files.

2.1.4 Model cross-validation and PMoA selection

To cross-validate the model, we attempted to predict the effects at GS 46 from effects on larvae up to GS 42. To do so, we used the point estimates of parameters inferred for PMoAs which plausibly reproduced the data during parameter inference. From the raw simulation output, i.e. the solution of the ODE system 11, we extracted the observed metamorphosis traits, mimicking the fixed-stage design.

We also used this part of the data to narrow down the selection of the PMoA, by comparing the predictions made by different PMoAs visually.

2.2 Model validation with mixture data

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:

• Predicting the effects of the binary mixture on larvae up to GS 42 from the corresponding single-stressor data.
• Predicting the effects of the binary mixture at GS 46 from the single-stressor data on larvae up to GS 42.

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:

(21) MAPE=100·[∑i=1n(y^i−yiyi)n−1]

y_iis 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.

3 Results

3.1 Model calibration

3.1.1 Effects of 2,4-D on larvae

Exposure to 2,4-D resulted in stagnation of growth at 30 mg/L, accompanied by a clear delay in metamorphosis. We fitted models based on PMoAs G (decrease in growth efficiency), M (increase in maintenance costs) and A (decrease in assimilation efficiency), and found that only PMoAs M and A can approximate the observations (Figure 1). PMoA G could not reproduce this pattern, because even a severe effect of growth implies only a slight delay of metamorphosis. PMoA A explained the timing of metamorphosis slightly better, whereas M explained effects on growth better. Based on this calibration, only PMoA G was ruled out. Quantitative evaluation supported this choice (Supporting Information Figure 1).

3.1.2 Effects of Flupyradifurone on larvae

Flupyradifurone caused effects of growth at 100 mg/L, but without a clear delay in metamorphosis. Based on the calibration to 2,4-D data, we could rule out PMoAs M and A, because both imply a delay in metamorphosis when growth is affected. Since the timing of metamorphosis was slightly earlier in the affected treatment, we also considered PMoA κ, which causes increased investment in maturation. PMoA G clearly explained effects on growth better, while implying a slight delay in metamorphosis (Figure 2). PMoA κ over-predicted the effects on growth, causing severe shrinking at high concentrations, while also over-predicting the effects on premature metamorphosis. Based on this calibration, PMoA G is more plausible. This was also reflected in quantitative metrics (Supporting Information Figure 2).

3.2 Cross-validation: Effects of 2,4-D at GS 46

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.

Table 3. TKTD parameter estimates for 2,4-D and PMoA M (increase in maintenance costs).

3.3 Cross-validation: Effects of Flupyradifurone at GS 46

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.

Table 4. TKTD parameter estimates for Flupyradifurone and PMoA G (decrease in growth efficiency).

3.4 Model validation: Predicting mixture effects

3.4.1 Predicting mixture effects on larvae

Based on the insights from calibration and cross-validation, we generated predictions of mixture toxicity by assuming PMoA M for 2,4-D and G for Flupyradifurone.

In the mixture treatments, the data indicated a dominance effect of 2,4-D with respect to growth and timing of metamorphosis. That is, additional exposure to Flupyradifurone did not alter the response to 2,4-D. This was predicted accurately by the model (Figures 5, 6). In the highest 2,4-D treatment, the model predicted severe shrinking and a complete inability to reach metamorphosis, which is generally in line with the observation that this treatment caused 100% mortality.

3.4.2 Predicting mixture effects on metamorphs

Similarly to the result of the cross-validation, the model predicted mixture effects on the timing of GS 46 with high accuracy, but low precision (Figure 7). This coincided with a high mortality at the corresponding treatments. The model predictions imply that individuals, if they would not have died from the exposure, may have taken a long time to complete metamorphosis, given that they reach GS 42 at all. Remarkably, the predicted variability in the timing of GS 46 was exclusively the result of the simulated individual variability, with the degree of variability estimated from control data (Hansul et al., 2026). This is different from the potential variability in simulation outputs that could result from the propagation of parameter uncertainty. The model under-predicted the effects on mass at GS 46 at the higher mixture treatment. When we compare this to the predictions of single-stressor effects, it becomes apparent that this mostly reflects a tendency of the model to under-predict effects on mass at GS 46, rather than under-predict the mixture effect.

4 Discussion

4.1 Calibration and cross-validation: Inference of PMoAs from life history data

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 η_ASdoes 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_Dcan 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)).

4.2 Validation: Predicting mixture effects

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.

4.3 Potential implications for amphibian population dynamics

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.

4.4 Conclusions

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.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/famrs.2026.1909160/full#supplementary-material

Quantitative metrics for physiological mode of action selection. Supplementary material at the source

Footnotes

Edited by: Giulia Simbula, Universidade do Porto, Portugal

Reviewed by: Miguel A: Carretero, Centro de Investigacao em Biodiversidade e Recursos Geneticos (CIBIO-InBIO), Portugal

Muammer Kurnaz, Gumushane University, Türkiye