A Bayesian Neural Ordinary Differential Equations Framework to Study the Effects of Chemical Mixtures on Survival

A Bayesian Neural Ordinary Differential Equations Framework to Study the Effects of Chemical Mixtures on Survival | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results A Bayesian Neural Ordinary Differential Equations Framework to Study the Effects of Chemical Mixtures on Survival View ORCID Profile Virgile Baudrot , View ORCID Profile Nina Cedergreen , Thomas Kleiber , View ORCID Profile André Gergs , View ORCID Profile Sandrine Charles doi: https://doi.org/10.1101/2025.03.14.643272 Virgile Baudrot 1 Qonfluens , Montpellier, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Virgile Baudrot For correspondence: virgile.baudrot{at}qonfluens.com Nina Cedergreen 2 University of Copenhagen , Copenhagen, Denmark Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Nina Cedergreen Thomas Kleiber 1 Qonfluens , Montpellier, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site André Gergs 3 Bayer Crop Sciences , Monheim, Germany Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for André Gergs Sandrine Charles 4 Université Claude Bernard Lyon 1 , Villeurbanne, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Sandrine Charles Abstract Full Text Info/History Metrics Supplementary material Preview PDF Abstract In the field of Plant Protection Products (PPP), the combination of multiple active ingredients is a common strategy to improve the efficacy against target species. However, these mixtures may expose non-target species to untested chemical combinations due to varied application and degradation pathways. Developing a method to predict the effects of such mixtures without extensive testing is crucial for enhancing risk assessment for non-target organisms. Toxicokinetic (TK) and toxicodynamic (TD) models currently provide predictive power for the toxicity of mixture over time, utilizing the concentration addition and independent action principles. However, their capacity to capture complex, non-linear interactions — such as synergistic and antagonistic effects — remains limited. To overcome these limitations, we introduce a new methodology that integrates TKTD models using Ordinary Differential Equations (ODE) with neural networks (NN), offering a robust framework for modeling complex substance interactions. The ODE component encodes fundamental biological principles that constrain the neural network to biologically plausible solutions. Bayesian inference further refines the model by addressing uncertainties in data, model parameters, and biological processes, while also allowing quantification of prediction uncertainty, highlighting gaps in experimental data that warrant further investigation. This hybrid model was evaluated across 99 acute toxicity studies that included various PPP mixtures, testing its ability to identify and forecast deviations from expected mixture behaviors. Our findings highlight the potential of combining mechanistic models with machine learning techniques to advance predictive accuracy in environmental toxicology, providing a framework to anticipate the combined effects of untested chemical mixtures and inform risk assessment practices. 1 Introduction The use of synthetic pesticides has expanded significantly since its inception in the late 1930s, both in total volume and in the diversity of available active ingredients [ Bernhardt et al., 2017 ]. Today, more than 500 active ingredients are approved in Europe, distributed in thousands of commercial products [ European Commission, 2024 ], resulting in a vast array of chemical mixtures present in ecosystems [ Schreiner et al., 2016 ]. These substances may also spread through trophic networks, affecting a wide range of non-target organisms [ Fritsch et al., 2022 ]. Active ingredient mixtures arise through multiple pathways: the use of different products in various crops, the combination of products within a single crop, or the incorporation of multiple active ingredients within a single product [ Schreiner et al., 2016 ]. While formulated mixtures are generally evaluated for efficacy on target species and regulated for selected non-target species, the effects of unformulated, incidental mixtures in the environment are less well understood. Given the extensive range of substances, concentrations, and combinations involved, it is impossible to empirically test every possible mixture and concentration ratio [ Belden et al., 2007 ]. For mixtures, the joint effects on exposed organisms are frequently described by two main toxicity models: Concentration Addition (CA) and Independent Action (IA). These models provide reasonable estimates of the effects of the mixture when single-component toxicity data is available for the organism of interest. Reviews indicate that these models accurately predict the effects of joint pesticides in approximately 90% of cases [ Belden et al., 2007 , Cedergreen, 2014 , Martin et al., 2021 ]. However, the remaining cases exhibit interactions that lead to greater or lesser toxicity than predicted, known as synergism and antagonism, respectively. Such interactions often occur when one chemical (the active substance or the co-formulant) affects the bioavailability, uptake, distribution, or biotransformation of another [ Cedergreen, 2014 ]. Predicting the joint effects of active ingredients is critical for environmental risk assessments, as most plant protection products consist of multiple active ingredients designed to broaden their efficacy and delay resistance. The course of toxic effects over time provides valuable insights for both single substances and chemical mixtures. Toxicokinetic (TK) and toxicodynamic (TD) models are particularly useful here, as they account for the dynamic nature of toxicity. A widely used framework, the General Unified Threshold model of Survival (GUTS), offers a robust method to assess time-dependent mortality trends [ Jager et al., 2011 , Jager and Ashauer, 2023 ], and recent adaptations have extended GUTS to analyze mixture toxicity [ Ashauer et al., 2017 , Cedergreen et al., 2017 , Bart et al., 2021 ]. These adaptations aggregate the internal damage caused by each active ingredient, producing a single latent variable to describe the toxicodynamic state of the organism, a concept known as “Damage Addition” [ Ashauer et al., 2017 , Cedergreen et al., 2017 , Bart et al., 2021 , Cedergreen et al., 2021 ]. In this approach, the internal damage of each substance is combined additively, providing an input to the TD component of the framework. This additive approach allows for predictions across varied exposure profiles, although it does not capture non-additive interactions like synergy and antagonism. To refine this approach, some models have adjusted key parameters to represent deviations from additivity. For example, Ashauer et al. [2017] generalizes “Damage Addition” to n active ingredients using a weighted sum of TK outputs, while Cedergreen et al. [2017] introduce a deviation parameter that affects rate constants to reflect how one substance may alter the elimination or damage repair rates of another. However, these approaches rely on assumptions that (1) the CA concept implemented through damage addition is suitable for integrating mixtures in time-resolved models, and (2) deviations are likely to stem from interactions with rate-limiting constants, assumptions which remain inadequately tested. To address these limitations, we propose an alternative approach that minimizes predefined assumptions, utilizing a free-constraint model for damage aggregation between the TK and TD stages by means of Neural Networks (NN). This hybrid model couples the mechanistic TKTD parts of GUTS with a flexible, data-driven NN component to capture complex, non-linear interactions without imposing rigid mechanistic constraints. This “knowledge-embedded” approach, also known as physics-informed neural networks [ Raissi et al., 2019 ] or Neural ODE [ Chen et al., 2018 ], allows us to retain mechanistic insights for well-understood subcomponents ( e . g ., absorption processes in TK, mortality mechanisms in TD) while treating more ambiguous interactions as flexible, data-driven functions [ Glassey and Von Stosch, 2018 , Gruber et al., 2024 ]. The NN, constrained by ODE-driven biological laws, is effectively calibrated with large data sets to approximate the complex functional relationships between TK and TD phases. We designed this hybrid model within a Bayesian framework to incorporate and quantify uncertainties in data, parameters, and biological processes. Bayesian inference, by framing model dependencies as conditional probabilities, allows robust calibration with standard experimental toxicity data (OECD protocols), while also identifying areas where further data collection is needed to refine predictions [ Kaikkonen et al., 2021 , Lee, 2024 ]. In this regard, our study has three main objectives: (1) implement a Bayesian-calibrated hybrid model compatible with routine toxicity data; (2) evaluate the predictive capacity of the model in mixtures, including variations in component ratios or entirely new combinations of active ingredients; and (3) assess the model’s ability to quantify deviations from CA and IA, characterizing potential synergistic or antagonistic tendencies within these mixtures. To this end, we tested the model using 99 studies on fish survival conducted according to the OECD 203 guidelines [ OECD, 2019 ], covering herbicides, fungicides, and insecticides. Although the parameters of the TKTD model were based on previous work, the NN required extensive testing across multiple configurations, exploring the layer structure, parameter size, and activation functions. This calibration-validation process, with subsets for model calibration and prediction validation, reflects the challenges inherent to data-driven models constrained by smaller effective training sets. Overall, we expect that this hybrid framework advances mixture toxicity modelling by integrating mechanistic insights with the adaptability of neural networks, enhancing the predictive accuracy of untested chemical combinations and supporting more reliable environmental risk assessment. 2 Methods 2.1 Data sets The inference process was tested on 99 acute toxicity studies (data in the Supplementary Material) with rainbow trout (according to the OECD 203 standards) [ OECD, 2019 ], divided into four data sets comprising 38, 19, 21, and 21 studies, respectively. The experiments consisted of 96 hours of exposure, in which survival was monitored daily over a range of 4 to 8 concentrations tested, constant over time. For most of the data sets, the concentration range allowed one to cover low to high mortality. Each data set facilitated the analysis of interactions among 5 to 6 different active ingredients. The data set 1 consists of 244 time-series covering 21 mixtures, resulting in 1482 survival data points. The active ingredients are: bixafen, fluopyram, prothioconazole, spiroxamine, tebuconazole and trifloxystrobin. For instance, in this large data set, there are 13 times series with only bixafen as active ingredients, 6 time-series of a mixture of bixafen and fluopyram, 6 other times-series with the three active ingredients, bixafen, fluopyram, prothioconazole, another mixture with prothioconazole and spiroxamine, and others as detailed in the Supplementary Material. In the same way, the data set 2 consists of 120 time-series covering 9 mixtures, resulting in 662 survival data points. The active ingredients are deltamethrin, flupyradifurone, spiromesifen, tralomethrin and triazophos. The data set 3 consists of 130 time-series covering 11 mixtures, given 727 survival data points. The active ingredients are beta-cyfluthrin, clothianidin, cyfluthrin, imidacloprid, thiacloprid and thiodicarb. Finally, the data set 4 consists of 126 time-series covering 12 mixtures, given 758 survival data points. The active ingredients are aclonifen, diflufenican, flufenacet, flurtamone and metribuzin. Following guidelines, the survival endpoint refers to no visible movements such as gill movements and no reaction when touching of the caudal peduncle, while observation of immobility was considered as abnormal swimming behaviour. To further evaluate how model calibrations under the Concentration Addition hypothesis would behave in scenarios of synergism or antagonism compared to purely additive interactions, we generated three artificial data sets, as illustrated in Figure 2 . Download figure Open in new tab Fig 1. Graphical Abstract. Download figure Open in new tab Fig 2. Artificial data sets were generated for a mixture containing two active ingredients, A and B, with interaction patterns represented by blue effect isoclines, indicating equal survival probabilities over time (0, 4, 24, 48, 72, and 96 hours, displayed across columns). The concentration and constituent proportions of the mixture are illustrated by dots. When compared to Concentration Addition hypothesis, the top row displays the additive interaction pattern, the middle row shows antagonism, and the bottom row depicts synergism. 2.1.1 Potential effect of co-formulant To ensure that co-formulants do not affect survival, we compared the LC50 predicted with the formulation to the LC50 predicted with the active ingredients alone. This comparison was conducted using single-compound formulations that share the same co-formulants as those used in formulations containing mixtures of several active ingredients. No additional effects of the co-formulants were detected (see Supplementary Material). 2.2 Model The model is organized into three submodels, as illustrated in Figure 3 . The first and third submodels correspond to the toxicokinetic (TK) and toxicodynamic (TD) components, respectively. The TK sub-model describes the dynamics of the internal chemical concentration within the organisms, while the TD sub-model translates this internal concentration into biological damage and calculates the survival probability of the organisms. Both sub-models utilize mixtures from the reduced General Unified Threshold model of Survival (GUTS) [ Jager et al., 2011 ]. The second intermediary submodel bridges TK and TD by translating the TK output (internal chemical concentrations) into a single input for the TD model (biological damage) using several aggregation mixtures. Throughout this paper, we use the term “aggregation” of substances to describe this intermediary process, distinguishing it from “interaction” which refers to the biochemical interactions between substances. Download figure Open in new tab Fig 3. Diagram of the generic models illustrating the sequential steps, beginning with (1) the toxicokinetic modeling, which describes the uptake and elimination of n compounds; followed by (2) the bridging component between the toxicokinetic and toxicodynamic sections, which aggregates the compounds into a one-dimensional, time-dependent damage variable, D ( t ); and concluding with (3) the toxicodynamic modeling that translates D ( t ) into survival probability. 2.2.1 Toxicokinetics The toxicokinetic (TK) component of the model is mechanistic, defined by Ordinary Differential Equations (ODE) to capture the time-dependent dynamics of the exposure profiles. For each substance i with an exposure profile X i ( t ) in water, the TK model describes the kinetics of the substance as a proxy for the scaled internal concentration, denoted C i ( t ), since direct measurements of internal concentrations are not available. This model includes the dominant toxicokinetic rate constant, k d,i , which is specific to each substance i . The TK model is represented by the following equation: 2.2.2 TKTD bridge The second part of the model serves as the bridge between the TK and TD components, with multiple variants detailed in Table 1 . View this table: View inline View popup Download powerpoint Table 1. Labels and mathematical equations for the TKTD bridge models used in calibration and depicted in Figure 4 . The bridging function operates on the scaled internal concentration of each substance, C i ( t ). To calculate the damage, there are two possible approaches: using either the derivative of C i ( t ), or the variable itself. For n exposure compounds, using the derivative approach yields the following: where f d () is a function of the derivatives for each substance i . Alternatively, working directly with the concentrations gives: where f () is a function of all concentrations of substances. We chose the second option as a straightforward rule that ensures that the variable D ( t ) remains positive, which we achieved by applying an exponential transformation directly to D ( t ). Linear model If the model is linear, denoted n in Table 1 , the equations are similar to those used in previous studies [ Ashauer et al., 2017 , Cedergreen et al., 2017 , Bart et al., 2021 ], with the addition of a bias b to obtain a true linear model. Neural Network To address non-linear aggregations, we implemented a Neural Network (NN) based on the Universal Approximation Theorem, which states that a two-layer NN with a non-linear activation function can approximate any function [ Cybenko, 1989 ]. The input for the TKTD bridge is the scaled internal concentration vector of each substance i , represented as C ( t ) = ( C 1 ( t ), …, C i ( t ), … C n ( t )), while the output is the calculated damage, D ( t ). The first layer performs a linear aggregation of input features, adding a bias term: where W 1 = { w i,j } i,j =1, n is an n × n weight matrix, and b 1 = { b 1, i } i =1, n the bias vector of size n . This yields z 1 ( t ), a vector of size n with each coordinate given by: . At this stage, the model remains linear. To introduce non-linearity, we apply an activation function, σ 1 , to the resulting z 1 ( t ): In the final model configuration, we used a rectified linear unit function (ReLU) as the activation function, defined as σ 1 ( z 1 ( t )) = max(0, z 1 ( t )). The combination of the linear operation in (4) with the activation function, (5), creates a perceptron defining the first layer of the neural network as L 1 ( C ( t )) = σ 1 ( W 1 C ( t ) + b 1 ). For the second layer, a similar approach is applied: All models used for calibration are presented in Table 1 . 2.2.3 Toxicodynamics The component of Toxicodynamics (TD) is the final stage of the model, receiving as input the latent damage variable D ( t ), derived from the combined output variables of the component of Toxicokinetics (TK) through the TKTD bridge (see Section 2.2.2 ). This aggregate damage, D ( t ), represents the cumulative marginal damage D i ( t ) for each active ingredient i . The Individual Tolerance version of the General Unified Threshold model of Survival (GUTS-IT) framework is then employed, incorporating both the baseline survival probability S b ( t ) and the additional survival probability due to exposure to each substance i . The role of the TD component is only to convert cumulative damage D ( t ) into a global survival probability S ( t ), expressed as follows: This equation satisfies the initial conditions for damage, specifically that the damage is zero at the beginning of the experiment, D ( t = 0) = 0, and remains positive for all time points t, D ( t ) ≥ 0 ∀ t ). In its simplest form, the TD component relies on three parameters: background mortality, h b , median α and shape β of the threshold distribution under the GUTS-IT framework. Splitting α In a previous study, we demonstrated that the parameter α is correlated with the parameter that represents the no-effect concentration [ Baudrot and Charles, 2019 ]. To account for such a correlation within the context of mixtures, we propose an alternative model, termed “splitting”, in which α is expressed as a function of the marginal α i for each individual compound i . This approach utilizes the same function used to model damage from scaled internal concentrations, as follows: where f () is the same function as the one defined in Equation (3) . 2.3 Interaction hypotheses 2.3.1 Concentration Addition (CA) Let us consider two compounds A and B , with concentrations C A and C B . The concept of Concentration Addition (CA) states that the combined effect of the two compounds can be predicted by adding their concentrations in proportion to their isolated effect concentrations. Mathematically writing, denoting s () the survival probability, then the lethal concentrations that lead to X % of effects, LC X,A and LC X,B , are defined as s ( LC X,A ) = s ( LC X,B ) = 1 − X/ 100. The CA hypothesis is satisfied if for any proportion p : s ( p LC X,A + (1 − p ) LC X,B ) = s ( LC X,A ) = s ( LC X,B ). Therefore, assuming s ( LC X,A ) = s ( LC X,B ) = s ( C A , C B ), then CA explains the mixture if 1. This approach has been generalized by Berenbaum [1985] and Faust et al. [2000] to provide a generalized CA hypothesis with n compounds and to measure the deviation from it: In the situation where survival under exposure to the mixture is lower than predicted for compound A alone, s ( p LC X,A + (1 − p ) LC X,B ) < s ( A ). This effect is termed “CA synergism”, where the observed effect exceeds what would be expected by CA. If survival is higher than predicted, this is referred to as “CA antagonism”, and is expressed as s ( p, LC X,A + (1 − p ), LC X,B ) > s ( A ). The “Damage Addition” approach, as applied in previous studies [ Ashauer et al., 2017 , Cedergreen et al., 2017 , Bart et al., 2021 ] (model n in Table 1 ), is based on the sum of the internal concentrations estimated by the TK component, formulated as .This aggregate value is then scaled according to the effect of each chemical using the values LC X , as described above. Using the simplified TK model of the reduced GUTS, where there is only one parameter per substance, k d,i , relating internal to external concentration, a single k d value can be found, allowing multiple substances to share the same TK characteristics. Thus, in our model, “Damage Addition” aligns directly with CA, along with all of its underlying assumptions. 2.3.2 Independent Action (IA) The IA hypothesis is expressed by the equation s ( C A + C B ) = s ( C A )× s ( C B ), which is often represented in terms of effect e as: e ( C A + C B ) = 1 − (1 − e ( C A )) * (1 − e ( C B )) [ Faust et al., 2000 ]. Accordingly, deviation from the IA hypothesis is defined differently than for CA. Specifically, there is “IA synergism” if s ( C A + C B ) s ( C A ) × s ( C B ). Faust et al. [2000] provided a general definition for the deviation in effect when considering n compounds as follows: 2.3.3 Deviation To quantify deviations from the CA and IA hypotheses, we adopted a deviation equation used in previous studies [ Faust et al., 2000 , Belden et al., 2007 , Cedergreen et al., 2017 ], where a ratio is defined between the expected outcomes under each hypothesis and the observed outcomes. To facilitate comparisons across multiple mixtures at varying time points and effect levels, we defined a symmetric and normalized deviation ratio, expressed as follows: The term observed i refers to the value calculated by our model, while predicted i corresponds to the expected value under the CA or IA hypotheses. Given that both observed i and predicted i are positive, the resulting deviation falls within the interval [−1, 1]. A deviation of 0 indicates the match to the CA or IA hypothesis; a negative deviation indicates a tendency toward antagonism and a positive deviation suggests synergism. The term uncertainty i represents the uncertainty of predicted i , estimated by our model within a Bayesian framework and captured by the range between the 97.5 and 2.5% quantiles of the posterior distribution. 2.4 Implementation Simulation and inference were performed using two programming approaches: one in Julia [ Bezanson et al., 2017 ] and the other in R [ R Core Team, 2024 ]. For the Julia implementation, we utilized the Scientific Machine Learning suite (SciML), specifically the DifferenceEquations.jl ODE component [ Rackauckas and Nie, 2017 ], the Lux.jl library for the neural network component [ Pal, 2023 ], and the Turing.jl library for Bayesian inference across the entire model [ Ge et al., 2018 ]. For the R-based approach, we employed the deSolve R-package for simulation [ Soetaert et al., 2010 ] and the JAGS software [ Plummer, 2004 ] via the rjags R-package for Bayesian inference [ Plummer, 2023 ]. All scripts and data are available in the Supplementary Material. 2.5 Goodness-of-fit To assess the quality of the calibration output, we used the model performance criteria recommended by EFSA in its Scientific Opinion on TKTD modeling [ EFSA et al., 2018 ], including the Posterior Prediction Check (PPC), the Normalized Root Mean Square Error (NRMSE), and the Survival Probability Prediction Error (SPPE). The PPC enables a comparison between the predicted and observed survival counts, summarizing the predictions by their mean and the 95% uncertainty range. The percentage of predictions that capture the observed data serves as a quantitative goodness-of-fit criterion. The NRMSE is defined by the equation: where y obs,i and y pred,i , i = 1, N , are the observed and predicted values, respectively, across N data points. The SPPE focuses on the survival probabilities at the end of the experiment, quantifying the accuracy of the model at the end of the exposure period. It is defined as: where and represent the observed and predicted survival counts at the end of the exposure period. Since Bayesian inference was applied, we included the Widely Applicable Information Criterion (WAIC) to estimate the effective number of parameters, accounting for overfitting [ Watanabe, 2013 ]. The WAIC is defined as: where lpd is the log-point-wise predictive density and pWAIC represents the effective number of parameters. 2.6 Calibration During the calibration phase, we fitted the four data sets using eight different models, each sharing the same TK and TD components but differing in their specific TKTD bridge expressions, as described in Table 1 . The simplest mode model, denoted n in Table 1 , is linear. From this baseline, the complexity of the model was incrementally increased by adding one or more layers, incorporating weight matrices, biases, and activation functions. 2.7 Cross-Validation The cross-validation phase involved dividing each data set into two subsets: one for calibration, known as the training data set, used for parameter inference, and the other for validation, where observed data were compared with model predictions post-training. We implemented two approaches for splitting the data sets into training and validation subsets. In the random time-point approach, 90% of the data points were chosen randomly for calibration, while the remaining 10% were reserved for validation. In the random time-series approach, selection occurred at the level of entire time series, maintaining the same 90/10% ratio for training and validation (see Supplementary Material). 2.8 Validation on independent mixtures We conducted an additional validation phase using independent mixtures, specifically data sets that were not included in the calibration phase. To achieve this, we removed all time series corresponding to mixtures containing pairs of active ingredients prior to the calibration phase. The model was then fitted to the remaining training data set, and the predictions were subsequently compared with the removed time series during the validation phase (Figure in Supplementary Material). For illustration, consider a data set with three active ingredients: A, B, and C. The complete data set would likely include survival time series for various mixtures with one, two, or three active ingredients: A, B, C, AB, AC, BC, ABC . If we designate the pair AB for removal from the training subset, leaving A, B, C, AC, BC, ABC , then the validation phase would involve comparing the model predictions with the observed survival under the AB exposure profiles. 2.9 Prediction Unlike the validation phase, the prediction phase lacks observational data to assess the prediction quality. There are three types of prediction (see Supplementary Material): (1) predictions for the same mixture as in calibration but at new time points; (2) predictions for the same mixture as in calibration but with different exposure ratios; and (3) predictions for newly composed mixtures of tested active ingredients not yet tested. 3 Results 3.1 Calibration Model calibration was performed separately for the three artificial data sets and the four empirical data sets. Figure 4 provides an overview of the quality of model fit based on two goodness-of-fit measures: the Widely Applicable Information Criterion (WAIC) and the Normalized Root Mean Square Error (NRMSE). Models with dots located in the bottom left corner of each graph indicate a good fit. Download figure Open in new tab Fig 4. Scatter plots illustrate the quality of calibration results for three artificial data sets fitted to five generic models (upper panel) and for four experimental data sets fitted to eight different models (lower panel). The WAIC ( x -axis) and NRMSE ( y -axis) serve as goodness-of-fit metrics, with lower values indicating higher fit quality . In the “Models of aggregations” legend, n stands for the simplest linear model with n compounds ; n exp for the linear model with an exponential transformation ; n exp split is the n exp model with the splitting of α ; nn n exp is the application of a perceptron matrix of size n × n without activation function before the n exp model; and nn ReLU n exp is the application of an n × n layer with a ReLU function before the n exp model. The remaining models with the nn ReLU prefix involve a sequence of perceptron layers before the n exp model. See Table 1 for details. The data sets called “add”, “ant” and “syn” were simulated under additive, antagonism and synergism interaction hypotheses, respectively. See Figure 2 for details. For artificial data sets (upper panel in Figure 4 , the linear models n and n_exp best fit the data simulated with the additive model. In contrast, these two models perform poorly when calibrated on antagonistic and synergistic data, displaying high NRMSE and WAIC values. The most complex model, with three NN layers (nn_ReLU_nn_ReLU_nn_ReLU_n_exp) , is the least effective for the additive data but performs best for the antagonistic and synergistic data ( Figure 4 ). Although NRMSE quantifies the fit between observations and predictions, WAIC accounts for the cost associated with the number of model parameters. As shown in Figure 4 , calibration results for artificial data reveal that the high parameter demand and computation time of the NN models is effectively offset when non-linearities are involved in compound aggregation. In simple terms, although NN models are complex and time-intensive, their application proves highly advantageous and thus warranted when handling mixture data with non-linear aggregations (e.g., involving antagonism and synergism). In contrast, when the concentration addition (CA) hypothesis is met, using an NN model for damage aggregation becomes excessive, and a simple linear model suffices. Calibration on empirical data also yields NRMSE values, indicating that NN models consistently outperform linear models. When including a cost adjustment for model parameters (as in WAIC), linear models (blue dots in Figure 4 ) show slightly better performance for data sets 2 and 4, while performance is comparable for the others. To determine whether the empirical data satisfy the CA hypothesis, we compared the calibration performance in Figure 5 with the results of the artificial data ( Figure 4 ). This comparison suggests the presence of non-linearities in each of the four empirical data, as NN models consistently perform better. If the CA hypothesis were largely valid, models with linear aggregation (without NN) would have demonstrated the best performance. Download figure Open in new tab Fig 5. Dots represent NRMSE values obtained by fitting the four models (distinguished by different colors) across the four experimental data sets: lower NRMSE values indicate higher-quality cross-validation results . The x-axis is not quantitative (points are spread). Two sampling approaches are shown: random selection of time-points (circle dots) and random selection of time-series (triangle dots), each with 90% of the data set used for calibration (training) and 10% for validation. Each cross-validation was repeated five times, illustrated by five identical points. 3.2 Validation The cross-validation results ( Figure 5 ) indicate that the random time-series validation approach yields significantly better outcomes compared to the random time-points approach (see Supplementary Material). For both cross-validation approaches and across all four experimental data sets, we applied the top four models to perform cross-validation, with each approach repeated five times. Using NRMSE, the variability across these five repetitions ( Figure 5 ) demonstrates that no model consistently outperforms the others, regardless of the data set. For the validation for completely unknown mixtures (see Supplementary Material), the results shown in Figure 6 indicate that models without NN components (blue bar) are the most effective, except in the cases of the flupyradifurone-spiromesifen pair in data set 2 and the diflufenican-metribuzin pair in data set 4. Download figure Open in new tab Fig 6. The horizontal bars represent the NRMSE values obtained for the four models (distinguished by different colors) when predicting entirely new mixtures composed of combinations of active ingredients not present in the experimental data sets. 3.3 Prediction The predictions use parameter values and their probability distributions inferred during the calibration phase performed under a Bayesian framework. These distributions are then used as inputs into the model to predict survival at various times, mixture proportions, and for new mixtures. An effect prediction and the associated uncertainties for different time points and a mixture of varying proportions of two chemicals are shown in Figure 7 for the example of deltamethrin and tralomethrin. Uncertainty is quantified by the extent of credible intervals resulting from the propagated uncertainties of parameter estimates obtained through Bayesian inference. We used the 95% credible intervals (CI) of the predicted survival rate as an uncertainty measure. The width of this interval serves directly as an uncertainty measure because the survival interval is bounded in [0, 1]. Download figure Open in new tab Fig 7. In the top panel, heat maps display the predicted survival probabilities for a binary mixture across different time points, using the neural network model nn ReLU n exp . Light blue lines indicate survival probability isoclines, where the shapes reflect patterns of interaction. In the lower panel, heat maps illustrate the uncertainty of the predictions, calculated as the range of their 95% credible intervals. In Figure 7 (middle and lower panels), white indicates low uncertainty (95% CI near 0), while purple represents high uncertainty (95% CI near 1). As expected, the uncertainty heat map at time 0 is entirely white. For later time points, white areas surround the points or occupy the upper right corners of heat maps where the survival probability drops to 0. The purple regions, where uncertainty peaks, indicate critical scenarios in which minor adjustments in the relative proportions of the compounds lead to substantial shifts in survival probabilities. This observation is particularly valuable, as it demonstrates that the uncertainty interval (95% CI here) serves as an effective metric for the quality of the prediction. In addition, complex models may occasionally produce predictions disconnected from the input data. This is especially relevant for hybrid models combining mechanistic components and neural networks. Due to the constraints imposed by the mechanistic model component, outputs are generally kept within a plausible range, making errors (sometimes termed “hallucinations”) harder to identify. This uncertainty measure enables the identification of such “hallucination” regions, or at least regions where further experimental data would improve predictions by reducing uncertainty. 3.3.1 Deviation in prediction The final results of this study, illustrated in Figure 8 , focus on the deviation from the CA and IA hypotheses. For each prediction ( Figure 7 ), we calculated, for each data point, a deviation relative to the expected results under the CA or IA hypotheses. We then calculated the mean deviation for all points, weighted by an uncertainty penalty (see Equation (11) ). To mitigate the influence of high uncertainty, we excluded regions where uncertainty exceeded 50% ( i . e ., where the 95% credible interval spanned more than 50%). Download figure Open in new tab Fig 8. Classification of interaction patterns based on deviations from the specified interaction hypotheses, either Concentration Addition (left panel) or Independent Action (right panel), for each mixture combining two active ingredients. Labels of the y-axis are red when no data exists in the dataset for these mixtures and green when some data points exists in our database. Missing point appear when 95% uncertainty covers more than 50% of the survival probability. Data points indicate the deviation calculations for all four aggregation models, and each line on the graphs represents one of the pairs tested. A value of zero on the x -axis indicates no deviation, signifying that the CA or IA hypothesis holds. In contrast, negative or positive values denote a potential for antagonistic or synergistic effects, respectively. This deviation metric ranges from 0 to 1, where 0 denotes no deviation, and -1 or 1 are the highest (asymptotically) deviations. For effect quantification, we selected thresholds of 0.1 for a weakly probable effect and 0.2 for a more certain effect. Across all predictions, there is no evidence of antagonistic effects and only a few instances of synergistic effects. The CA and IA hypotheses hold in general across all pairs. Two pairs, clothianidin-thiodicarb and clothianidin-beta-cyfluthrin, deviate from the CA hypothesis, suggesting a potential for synergistic effect ( Figure 8 ). For the IA hypothesis, synergistic effects were observed for the pairs trifloxystrobin-bixafen, flurtamone-aclonifen, and deltamethrin-tralomethrin, but with a large uncertainty in the prediction for all three couples. 4 Discussion This study pursued three primary objectives by developing and applying a hybrid model that integrates mechanistic toxicokinetic-toxicodynamic (TKTD) components with neural network (NN) elements, implemented within a Bayesian framework. Our first objective was to demonstrate the practical applicability of this model with standard toxicology data sets. We achieved this by applying it to 99 data sets on fish survival exposed to various types of pesticides, including herbicides, fungicides, and insecticides. Secondly, our goal was to evaluate the predictive performance of the model for new mixtures, including those with varied active ingredient ratios and entirely novel mixtures. By conducting extensive calibration-validation exercises across different subsampling methods, we observed that incorporating NN elements within the TKTD framework improved calibration quality, but did not consistently enhance predictive accuracy for untested mixtures. Finally, we quantified potential synergistic or antagonistic interactions within these mixtures by assessing deviations from the Concentration Addition (CA) and Independent Action (IA) hypotheses. Our results showed that most of the compound combinations adhered to either the CA or the IA assumptions, in agreement with previous findings [ Belden et al., 2007 , Cedergreen, 2014 , Martin et al., 2021 ]. This linearity in the aggregation of active substances for the vast majority of cases studied in this paper explains the low added value of neural networks compared to classical models. Hence, datasets are needed that include substances that induce synergy or antagony, and ideally also include variable chemical exposures to fully evaluate the added value of the neural network approach. 4.1 Mixture toxicity with Neural-TKTD Modeling the aggregation of active ingredients from standardized toxicity tests presents substantial challenges due to limited data relative to the vast number of possible mixtures. To address these challenges, we developed a hybrid framework that combines mechanistic (ODE) and non-mechanistic (NN) modeling approaches, leveraging the strengths of each methodology [ Chen et al., 2018 , Lu et al., 2021 ] to enhance the accuracy and robustness of the predictions. Recent work in pharmacology (PKPD modeling) has demonstrated that integrating pharmacological principles into NN architectures (Neural-PKPD) enables simulations of responses to previously untested dosing regimens [ Lu et al., 2021 ]. For cross-validation, we tested different data partitioning strategies by splitting data sets into calibration (training) and validation subsets. We found that selecting the entire time series for validation, rather than randomly selecting individual time points, provided the best results. Despite this, model performance across validation subsets was largely equivalent. Calibration results suggest that NN exhibit a degree of overfitting, likely due to their high flexibility, which enables close fits to training data, but may capture noise rather than true underlying distributions. A promising outcome of this study is the straightforward protocol developed to build the model ( Figure 3 ) across data sets with varying mixtures and variable exposure profiles, using the number of active ingredients to set layer sizes and employing the ReLU activation function. Although the NN architecture is relatively simple by modern standards, it effectively meets our objectives. Using WAIC as a goodness-of-fit criterion, we found that more complex NN architectures (with additional layers and parameters) were less effective in capturing experimental data patterns. Increasingly, pharmacological studies use hybrid models to capture processes lacking explicit equations. For instance, Gruber et al. [2024] and Führer et al. [2024] implemented NN upstream of pharmacokinetics (PK) to incorporate molecular and mixture-specific parameters as PK inputs, enhancing drug candidate discovery. Given the parallels between TK and PK, future work will extend this approach within the TKTD framework. 4.2 Inference of model parameters For parameter estimation, we employed Bayesian inference, which incorporates prior knowledge and propagates uncertainties [ Baudrot and Charles, 2019 ], allowing for straightforward parameter updates as new data become available. Bayesian inference has recently been used to discern the interactions between covariates and PK parameters, while simultaneously assessing various levels of random effects and measuring uncertainty [ Elmokadem et al., 2024 ]. In artificial data sets, hybrid models with a single NN layer outperformed traditional linear models in detecting potential synergistic and antagonistic effects. As anticipated, linear models excelled in cases of additive aggregation, as they require fewer parameters for similar fit quality. However, in real data sets, hybrid models enhanced with NN provided superior calibration outputs compared to linear models ( Figure 4 ), reflecting the well-known flexibility of NN to fit diverse functions. The enhanced performance of NN on real data sets implies that certain compound mixtures do not follow simple additive mechanisms. Calibration-validation procedures identified non-additive interactions for specific com-pound pairs, with Figure 6 showing non-linear effects for some couples. However, note that in some of these combinations the tested concentration range did not allow strong conclusions, hence, a dataset with more pronounced interaction effects is needed to test the full potential of the model. 4.3 Predictability An additional question concerned whether the effects of new mixtures could be predicted from the marginal toxicities of the active ingredients. TKTD models, with or without neural networks, enable simulation of scenarios lacking empirical data, including varied exposure profiles, varying time points of assessments, assessing different active ingredient ratios, and potentially entirely new active ingredient combinations if they behave similarly to already tested compounds. However, particularly the latter assumption needs to be tested on empirical data. Previous models for evaluating synergism or antagonism in binary mixtures have relied on isobologram-based statistical tests [ Sørensen et al., 2007 ]. Our TKTD approach provides an advantage by enabling predictions across time points through interpolation or extrapolation, offering insights into mixture behaviour over time rather than being limited to specific experimental time points. This dynamic capability improves predictive power by capturing time-dependent interaction patterns. Using Bayesian inference for calibration allowed uncertainty propagation throughout the predictions. This approach effectively manages uncertainty, yielding low-uncertainty predictions where data support is robust, and avoids unrealistic predictions, such as increased survival at lethal concentrations. This capability confers two key benefits: (1) quantifying the reliability of the prediction in varying scenarios and (2) identifying data gaps, which can inform the cost-effective experiment design by prioritizing experiments with high information value given the significant costs associated with experimental toxicology. 4.4 Deviation from additivity hypotheses To test the concentration addition (CA) and independent action (IA) hypotheses, we calculated deviations of the model predictions from expectations under these hypotheses ( Figure 8 ). Deviations from these hypotheses indicate a potential for non-linear aggregation of active ingredients, which is critical for targeted risk assessment. The results suggest that the additivity hypotheses are generally valid, with synergistic effects appearing in the rare instances where deviations are observed. Antagonistic effects were absent, in line with previous studies on interactions of non-additive compounds, where synergistic effects predominate [ Cedergreen, 2014 , Martin et al., 2021 ]. There is no mechanistic foundation that antagonism is less common than synergism, but there is a bias: by design, synergistic mechanisms raise concerns for safety purposes (e.g., cocktail effects [ Altenburger et al., 2013 ]), which is not the case for antagonism. As a consequence, antagonism is likely to be less scrutinized in scientific studies. There is no mechanistic foundation that antagonism is less common than synergism, nor is it found less frequently than synergy when reviewing published data [ Cedergreen et al., 2008 ]. This is despite there being a bias in the literature of looking for synergistic interactions as synergistic mechanisms raise concerns for safety purposes if using a reference model to predict worst case effects, which is not the case for antagonism [ Martin et al., 2021 ] (e.g., cocktail effects [ Altenburger et al., 2013 ]). As a consequence, antagonism is likely to be less scrutinized in scientific studies. Deviations from the CA hypothesis were observed for two pairs, clothianidin thiodicarb and clothianidin-beta-cyfluthrin, which would suggest synergistic effects. However, the predictions for both combinations lack consistency (see Supplementary Material): e.g., the models predict a decrease in effect with increases in doses in both cases, accompanied by high uncertainty and lack of data for verification. Hence, such results are typically situations of hallucination highlighted by the propagation of uncertainties. Under the IA hypothesis, strong synergism was observed in the trifloxystrob bixafen, deltamethrin-tralomethrin, and flurtamone-aclonifen pairs. However, this is not uncommon and will always occur when concentration-response curves have steep slopes (¿ 1 for the log-logistic concentration-response model) [ Drescher and Boedeker, 1995 , Cedergreen et al., 2008 ] which is the main reason for CA being used as the preferred reference model in a risk assessment context. The discrepancies that we observed in the deviation between CA and IA models are in line with their underlying assumptions. In the CA model, all substances are summed after estimating the potency, while for the IA model, the responses given in proportions of the control are multiplied. At low concentrations, the mixture of substances below their individual threshold (NEC) still contributes to the mixture effect in CA while fully ignored in IA, while for binary mortality responses using IA, non-lethal effects are ignored. The difference between CA and IA predictions depends on the slope of the dose-response curve. In addition, based on their mathematical definitions, over all concentration ranges, the combined effect in IA tends generally to be lower than what CA would predict for the same concentrations. CA is more conservative compared to IA for safety purposes. Also, the detection of a synergy under IA (a stronger effect observed in data compared to the model prediction) can still be additivity (or even antagonism) under CA [ Drescher and Boedeker, 1995 , Altenburger et al., 2013 ]. We also observed a greater variability in deviation in the IA model compared to the CA models. This is largely because the IA model combines individual responses multiplicatively (or via a function that is non-linear) whereas concentration addition is a linear summation. As a consequence, any variability or uncertainty in the estimation of each individual effect (e.g. due to experimental error, biological variability, or model uncertainty) are much more amplified under IA compared to CA. These results, which illustrate a low deviation from the additive model, are due to the fact that 90% of the data sets validates the additive hypothesis [ Cedergreen, 2014 ]. Thus, from the present study, if we have shown that theoretically, from artificial data sets, we can retrieve their synergistic or antagonistic characteristics, we demonstrate that the method does not create false positives. 4.5 Conclusion Our framework presents a robust and accessible solution for researchers and practitioners in ecotoxicology, providing a powerful tool for the study and prediction of toxicity results in complex mixtures. By integrating the mechanistic components of TKTD with neural networks within a Bayesian framework, this model improves our ability to understand and quantify interactions among active ingredients, providing information on potential synergistic or antagonistic effects with greater confidence. Bayesian inference enables rigorous parameter estimation and effective uncertainty propagation, essential for reliable predictions and data-informed decision-making. Designed for usability, this advanced framework is both adaptable and practical, promoting its widespread adoption for applications in the evaluation and management of ecotoxicological risks. 5 Data availability All scripts and data are available in the Supplementary Material. The study report for empirical data is available for non-commercial use on request via email: cropscience-transparency{at}bayer.com. Please indicate in your request the test report numbers as listed in data tables in the Supplementary Material. 6 Declaration of interest A.G. is employed by Bayer AG, which produces and sells agrochemicals including those investigated in this study. Bayer AG applies TKTD models to assess the effects of chemicals on target and non-target species. The contributions of N.C., V.B., T.K., S.C. to the study were partially funded by Bayer AG (052-2023 Research Agreement between UCBL-Bayer-UCPH-Qonfluens). 7 Credit authorship contribution statement Conceptualisation : VB, SC ; Data Curation : AG, VB, TK ; Formal Analysis : VB, SC ; Funding Acquisition : AG, SC ; Methodology : VB, SC ; Project Administration : SC, AG ; Software : VB, TK ; Supervision : AG, SC ; Interpretation & Validation : NC, AG ; Visualisation : VB ; Writing – Original Draft Preparation : VB, SC ; Writing – Review & Editing : VB, SC, NC, AG. References ↵ R. Altenburger , T. Backhaus , W. Boedeker , M. Faust , and M. Scholze . Simplifying complexity: Mixture toxicity assessment in the last 20 years simplifying complexity: Mixture toxicity assessment in the last 20 years . Environmental Toxicology & Chemistry , 32 ( 8 ), 2013 . ↵ R. Ashauer , I. O’Connor , and B. I. Escher . Toxic mixtures in time: The sequence makes the poison . Environmental science & technology , 51 ( 5 ): 3084 – 3092 , 2017 . OpenUrl CrossRef PubMed ↵ S. Bart , T. Jager , A. Robinson , E. Lahive , D. J. Spurgeon , and R. Ashauer . Predicting mixture effects over time with toxicokinetic–toxicodynamic models (guts): Assumptions, experimental testing, and predictive power . Environmental Science & Technology , 55 ( 4 ): 2430 – 2439 , 2021 . OpenUrl CrossRef PubMed ↵ V. Baudrot and S. Charles . Recommendations to address uncertainties in environmental risk assessment using toxicokinetic-toxicodynamic models . Scientific reports , 9 ( 1 ): 11432 , 2019 . OpenUrl CrossRef PubMed ↵ J. B. Belden , R. J. Gilliom , and M. J. Lydy . How well can we predict the toxicity of pesticide mixtures to aquatic life? Integrated Environmental Assessment and Management: An International Journal , 3 ( 3 ): 364 – 372 , 2007 . OpenUrl ↵ M. C. Berenbaum . The expected effect of a combination of agents: the general solution . Journal of theoretical biology , 114 ( 3 ): 413 – 431 , 1985 . OpenUrl CrossRef PubMed Web of Science ↵ E. S. Bernhardt , E. J. Rosi , and M. O. Gessner . Synthetic chemicals as agents of global change . Frontiers in Ecology and the Environment , 15 ( 2 ): 84 – 90 , 2017 . OpenUrl CrossRef ↵ J. Bezanson , A. Edelman , S. Karpinski , and V. B. Shah . Julia: A fresh approach to numerical computing . SIAM review , 59 ( 1 ): 65 – 98 , 2017 . OpenUrl CrossRef ↵ N. Cedergreen . Quantifying synergy: a systematic review of mixture toxicity studies within environmental toxicology . PloS one , 9 ( 5 ): e96580 , 2014 . OpenUrl CrossRef PubMed ↵ N. Cedergreen , A. M. Christensen , A. Kamper , P. Kudsk , S. K. Mathiassen , J. C. Streibig , and H. Sørensen . A review of independent action compared to concentration addition as reference models for mixtures of compounds with different molecular target sites . Environmental toxicology and chemistry , 27 ( 7 ): 1621 – 1632 , 2008 . OpenUrl CrossRef PubMed ↵ N. Cedergreen , K. Dalhoff , D. Li , M. Gottardi , and A. C. Kretschmann . Can toxicokinetic and toxicodynamic modeling be used to understand and predict synergistic interactions between chemicals? Environmental science & technology , 51 ( 24 ): 14379 – 14389 , 2017 . OpenUrl CrossRef PubMed ↵ N. Cedergreen , G. Bellisai , L. Herrero-Nogareda , E. Boesen , and K. Dalhoff . Using tktd models in combination with in vivo enzyme inhibition assays to investigate the mechanisms behind synergistic interactions across two species . Environmental Science & Technology , 55 ( 20 ): 13990 – 13999 , 2021 . OpenUrl CrossRef PubMed ↵ R. T. Chen , Y. Rubanova , J. Bettencourt , and D. K. Duvenaud . Neural ordinary differential equations . Advances in neural information processing systems , 31 , 2018 . ↵ G. Cybenko . Approximation by superpositions of a sigmoidal function . Mathematics of control, signals and systems , 2 ( 4 ): 303 – 314 , 1989 . OpenUrl CrossRef ↵ K. Drescher and W. Boedeker . Assessment of the combined effects of substances: The relationship between concentration addition and independent action . Biometrics , pages 716 – 730 , 1995 . ↵ P. Efsa , C. Ockleford , P. Adriaanse , P. Berny , T. Brock , S. Duquesne , S. Grilli , A. F. Hernandez-Jerez , S. H. Bennekou , M. Klein , et al. Scientific opinion on the state of the art of toxicokinetic/toxicodynamic (tktd) effect models for regulatory risk assessment of pesticides for aquatic organisms. panel on plant protection products and their residues (ppr) . EFSA journal , 16 ( 8 ): e05377 , 2018 . OpenUrl ↵ A. Elmokadem , M. Wiens , T. Knab , K. Utsey , S. P. Callisto , and D. Kirouac . Hierarchical deep compartment modeling: A workflow to leverage machine learning and bayesian inference for hierarchical pharmacometric modeling . Clinical and Translational Science , 17 ( 10 ): e70045 , 2024 . OpenUrl CrossRef ↵ European Commission . Eu pesticides database , 2024 . URL https://ec.europa.eu/food/plants/pesticides/eu-pesticides-database . Accessed: 2024-06-07 . ↵ M. Faust , R. Altenburger , T. Backhaus , W. Bödeker , M. Scholze , and L. H. Grimme . Predictive Assessment of the Aquatic Toxicity of Multiple Chemical Mixtures . Journal of Environmental Quality , 29 ( 4 ): 1063 – 1068 , 2000 . doi: 10.2134/jeq2000.00472425002900040005x . OpenUrl CrossRef ↵ C. Fritsch , B. Appenzeller , L. Burkart , M. Coeurdassier , R. Scheifler , F. Raoul , V. Driget , T. Powolny , C. Gagnaison , D. Rieffel , et al. Pervasive exposure of wild small mammals to legacy and currently used pesticide mixtures in arable landscapes . Scientific Reports , 12 ( 1 ): 15904 , 2022 . OpenUrl CrossRef PubMed ↵ F. Führer , A. Gruber , H. Diedam , A. H. Göller , S. Menz , and S. Schneckener . A deep neural network: mechanistic hybrid model to predict pharmacokinetics in rat . Journal of computer-aided molecular design , 38 ( 1 ): 7 , 2024 . OpenUrl CrossRef PubMed ↵ H. Ge , K. Xu , and Z. Ghahramani . Turing: a language for flexible probabilistic inference . In International Conference on Artificial Intelligence and Statistics, AISTATS 2018, 9-11 April 2018, Playa Blanca, Lanzarote, Canary Islands, Spain , pages 1682 – 1690 , 2018 . URL http://proceedings.mlr.press/v84/ge18b.html . ↵ J. Glassey and M. Von Stosch . Hybrid modeling in process industries . CRC Press , 2018 . ↵ A. Gruber , F. Führer , S. Menz , H. Diedam , A. H. Göller , and S. Schneckener . Prediction of human pharmacokinetics from chemical structure: combining mechanistic modeling with machine learning . Journal of Pharmaceutical Sciences , 113 ( 1 ): 55 – 63 , 2024 . OpenUrl CrossRef ↵ T. Jager and R. Ashauer . Modelling survival under chemical stress. A comprehensive guide to the GUTS framework . Leanpub, 2023 . ISBN 978-1-9999705-0-5 . URL https://leanpub.com/guts_book . ↵ T. Jager , C. Albert , T. G. Preuss , and R. Ashauer . General unified threshold model of survival-a toxicokinetic-toxicodynamic framework for ecotoxicology . Environmental science & technology , 45 ( 7 ): 2529 – 2540 , 2011 . OpenUrl CrossRef PubMed ↵ L. Kaikkonen , T. Parviainen , M. Rahikainen , L. Uusitalo , and A. Lehikoinen . Bayesian networks in environmental risk assessment: A review . Integrated environmental assessment and management , 17 ( 1 ): 62 – 78 , 2021 . OpenUrl CrossRef ↵ S. Y. Lee . Using bayesian statistics in confirmatory clinical trials in the regulatory setting: a tutorial review . BMC Medical Research Methodology , 24 ( 1 ): 110 , 2024 . OpenUrl CrossRef PubMed ↵ J. Lu , B. Bender , J. Y. Jin , and Y. Guan . Deep learning prediction of patient response time course from early data via neural-pharmacokinetic/pharmacodynamic modelling . Nature machine intelligence , 3 ( 8 ): 696 – 704 , 2021 . OpenUrl CrossRef ↵ O. Martin , M. Scholze , S. Ermler , J. McPhie , S. K. Bopp , A. Kienzler , N. Parissis , and A. Kortenkamp . Ten years of research on synergisms and antagonisms in chemical mixtures: A systematic review and quantitative reappraisal of mixture studies . Environment international , 146 : 106206 , 2021 . OpenUrl CrossRef PubMed ↵ OECD . Test no. 203: Fish, acute toxicity test , 2019 . URL https://www.oecd-ilibrary.org/content/publication/9789264069961-en . ↵ A. Pal . Lux: Explicit Parameterization of Deep Neural Networks in Julia , April 2023 . URL doi: 10.5281/zenodo.7808904 . If you use this software, please cite it as below. OpenUrl CrossRef ↵ M. Plummer . Jags: Just another gibbs sampler , 2004 . ↵ M. Plummer . rjags: Bayesian Graphical Models using MCMC , 2023 . R package version 4 – 15 . ↵ R Core Team . R: A Language and Environment for Statistical Computing . R Foundation for Statistical Computing , Vienna, Austria , 2024 . ↵ C. Rackauckas and Q. Nie . Differentialequations.jl–a performant and feature-rich ecosystem for solving differential equations in julia . Journal of Open Research Software , 5 ( 1 ): 15 , 2017 . OpenUrl CrossRef ↵ M. Raissi , P. Perdikaris , and G. E. Karniadakis . Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations . Journal of Computational physics , 378 : 686 – 707 , 2019 . OpenUrl CrossRef ↵ V. C. Schreiner , E. Szoöcs , A. K. Bhowmik , M. G. Vijver , and R. B. Schäfer . Pesticide mixtures in streams of several european countries and the usa . Science of the Total Environment , 573 : 680 – 689 , 2016 . OpenUrl CrossRef PubMed ↵ K. Soetaert , T. Petzoldt , and R. W. Setzer . Solving differential equations in R: Package deSolve . Journal of Statistical Software , 33 ( 9 ): 1 – 25 , 2010 . doi: 10.18637/jss.v033.i09 . OpenUrl CrossRef ↵ H. Sørensen , N. Cedergreen , I. M. Skovgaard , and J. C. Streibig . An isobole-based statistical model and test for synergism/antagonism in binary mixture toxicity experiments . Environmental and Ecological Statistics , 14 : 383 – 397 , 2007 . OpenUrl CrossRef ↵ S. Watanabe . A widely applicable bayesian information criterion . Journal of Machine Learning Research , 14 ( 1 ): 867 – 897 , 2013 . OpenUrl View the discussion thread. Back to top Previous Next Posted March 17, 2025. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following A Bayesian Neural Ordinary Differential Equations Framework to Study the Effects of Chemical Mixtures on Survival Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share A Bayesian Neural Ordinary Differential Equations Framework to Study the Effects of Chemical Mixtures on Survival Virgile Baudrot , Nina Cedergreen , Thomas Kleiber , André Gergs , Sandrine Charles bioRxiv 2025.03.14.643272; doi: https://doi.org/10.1101/2025.03.14.643272 Share This Article: Copy Citation Tools A Bayesian Neural Ordinary Differential Equations Framework to Study the Effects of Chemical Mixtures on Survival Virgile Baudrot , Nina Cedergreen , Thomas Kleiber , André Gergs , Sandrine Charles bioRxiv 2025.03.14.643272; doi: https://doi.org/10.1101/2025.03.14.643272 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Pharmacology and Toxicology Subject Areas All Articles Animal Behavior and Cognition (7636) Biochemistry (17705) Bioengineering (13899) Bioinformatics (41967) Biophysics (21460) Cancer Biology (18600) Cell Biology (25526) Clinical Trials (138) Developmental Biology (13384) Ecology (19909) Epidemiology (2067) Evolutionary Biology (24326) Genetics (15613) Genomics (22512) Immunology (17740) Microbiology (40423) Molecular Biology (17193) Neuroscience (88645) Paleontology (667) Pathology (2835) Pharmacology and Toxicology (4825) Physiology (7647) Plant Biology (15159) Scientific Communication and Education (2046) Synthetic Biology (4302) Systems Biology (9825) Zoology (2271)