Systems modeling of chemical-induced mitochondrial toxicity: pilot integration with multiple assays data from three cell-lines
1 Introduction
Mitochondria are main bioenergetic organelles. They are subjected to chemical-induced toxicity, resulting in adverse organ outcomes like neurodegenerative diseases (EFSA WG EPI1 Members et al. 2018; Stel, Bennekou, et al. 2020). With molecular initiating events of inhibition to complex I-IV of the electron transport chain (ETC)，a couple of adverse outcome pathways (AOP) were developed like OECD-endorsed AOP:3 and another AOP:273, see the adverse outcome pathways (AOPN) in Figure 1.
Previous work in quantitative AOP from one case study on renal tox (Zgheib et al. 2019) highlighted better flexibility and extrapolation of quantitative systems modeling approach compared to concentration-response and Bayesian network approaches. Regarding mitochondrial tox in liver systems, MitoSym™ from DILISym® showed a systems model of oxidative phosphorylation processes using readout from in vitro effect assay like seahorse oxygen consumption assay (Y. Yang et al. 2015). Focusing on ~20 complex inhibitors and uncouplers, our previous modeling (H. Yang et al. 2021) showed a reduced but sytems models with a coherent structure can be derived and calibrated well to mitochondrial membrane potential (proxy of KE of mitochondrial dysfunction). This reduced model was also benchmarked with a detailed biophysical model (Beard 2005), showing satisfactory approximation to ion leaky perturbed by uncouplers. However this model was only validated using key event like membrane potential and only validated in HepG2 cell-lines. To quantify the mitochondrial tox AOP with inhibition to complex I being the MIE, a more recent paper utilized non-mechanistic systems modeling approach: concentration-response fitting to three cell-lines (LUHMES, HepG2, and RPTEC/TERT1). The non-systems approaches can capture data rich compounds like rotenone and deguelin. However, when validating the model with another 8 compounds (using LUHMES-specific assays), for half of the compounds, activities of neuron degeneration were underestimated. We speculated that the unsatisfactory behavior could arise from limited capacity of a purely empirical concentration-response regression approach with missing mechanistic understanding of key events. Furthermore, it is well-noted that systems modeling is more capable in extrapolation than other two: concentration-response curve regression and Bayesian nets, see Table 4 in a comparative study (Zgheib et al. 2019).
OUr present work explores systems modeling to capture temporal concentration-response relationship of key events from tissue-specific cell-lines and later key events like cell death or neuron degeneration, which is the proxy to the adversity at a tissue level. We presented the AOPN starting from MIE i.e. binding to complex I to the later key event like cell death or neuron degeneration, see Figure 2. The main objective of our present modeling work is to demonstrate that a coherent structure of quantitative AOP model can be calibrated from multiple cell-lines.
2 Methods
2.1 in-vitro assay data from three cell-lines
We extracted and curated the experimental data from the previous study using three cell-lines (Tebby et al. 2022).
In total, ten compounds with suspected inhibition to complex I were included in the experimental study in (Tebby et al. 2022). The compounds were Antimycin A, Azoxystrobin, Capsaicin, Carboxine, Cyazofamid, Fenamidone, Fenazaquin, Fenfuram, Fenpyroximate, Flutolanil, Hydramethylnon, Kresoxim-methyl, Mepronil, Picoxystrobin, Pyraclostrobin, Pyridaben, Pyrimidifen, Tebufenpyrad, Thifluzamide, and Trifloxystrobin. However, the availability of measurements about MIE and KEs differs per compound per cell type. For LUHMES, all 10 compounds were included, and rotenone and deguelin had most rich datasets with various key events. For other two cell types, HepG2 and RPTEC/TERT1, only rotenone and deguelin were included. Per readout, multiple concernetation levels were utilized in the in-vitro testing.
As noted in the oxidative phosphorylation processes covered MIE and very first key events but later key event like cell death only associated with one in-vitro assay: resazurin viability assay for HepG2 and RPTEC/TERT1, readout of neurite area for LHUMES. Proteostasis assay were only applicable to LHUMES. Hence in our pilot modeling study, we only considered representative MIE and/or another KE related to oxidative phosphorylation.
2.2 Quantitative systems modeling
We considered pharmacokinetics following a simple decay, while initial concentration depends on the administrated concentration of the compound. For compound PK, compounds had their specific pharmacokinetics parameters like decay rate constant \(\alpha\), Michaelis–Menten kinetics parameters \(V_{max}\) and \(Km\), which maps the applied concentration (\(C_{ex}\)) to target-site concentration \(C_{in} =\frac{V_{max}}{Km + C_{ex}}\).
To quantify the AOP, our previous model (H. Yang et al. 2021) focused on two variables in oxidative phosphorylation: oxygen level and mitochondrial membrane potential, providing a mechanistic yet streamlined framework upon ETC inhibitors and uncouplers. Here, as extension to impaired and cell death/neuronal degeneration, we modelled later effect/key events following transit compartment modeling proposed for pharmacodynamics field (Krzyzanski 2011; Friberg et al. 2002).
For the entire set of ODEs, we built models for three different cell types and all with the same model structure.In the qAOP we have 10 parameters per cell-type independent of compounds. These parameters characterize the key event relationship. We also considered the steady state constraints to the four states: oxygen level, mitochondrial membrane potential, proteostasis activity, and cell viability. The normalized activity \(r\) defined ratio between time-varying variable \(x(t)\) and corresponding initial condition \(x_0\) at their steady states.
2.3 Model simulation and calibration with in-vitro assay data
We performed model calibration for data measured from multiple assays for individual cell-lines. We also noted that the MIE/KE assays were conducted at different timepoint upon exposure. From LUHMES with 10 compounds, we had complex I activity at 30 mins, proteostasis and neurite area both at 24 hour. From RPTEC/TERT1 with 10 compounds, we had mito. membrane potential and viability both at 24 hour.
In model calibrate/optimization, we applied multiple-starting value strategy using latin-hyper sampling technique (Raue et al. 2013). For each starting-value, we performed maximal likelihood estimation using sensitivity equation technique Raue et al. (2013).
We utilized our OneHealthStudio platform to develop and calibrate the models. OneHealthStudio employs CVODES solver from SUNDIALS to solve the differential equations of qAOP model. Besides, OneHealthStudio incorporates reliable and AI-empowered parameter inference technique, which ensure this performance and scalability even for large-scale physiologically-based kinetics modeling (PBPK)-qAOP model, see also details in another case study of PBPK/PD modeling siRNA treatment in model-informed drug development with OneHealthStudio.
Detailed about the model parameterization and estimates will be released into dedicated supplemental material. We are deploying the model into shinyapps.io. Please message us for early access to the model and simulator.
3 Results
We presented the modeling and simulation outcome from our quantitative systems modeling of mitochondrial toxicity network for three cell-lines. For each cell-line, We first explored the alignment of prediction of calibrated models against in-vitro assay data. Then we dived into more detailed model-based concentration-response.
3.1 Model fitting performance with three cell-lines
In general model fitting align well with data, see Figure 3. As note for HepG2 and LHUMES celllines, We saw large variability of observed activities when predict values are close to 100 %. We suspected that for large variability per se from data, we will have a deep dive in the next section when checking concentration-response data.
3.2 Model-based concentration-response relationship
We presented detailed concentration-response relationship simulated from the calibrated models of individual cell-lines. We extrapolated the concentration levels \(10^{-1}\) fold lower and 10-fold higher to the bounds of applied concentrations.
Concentration-response prediction for LHUMES
From calibration with data from 10 compounds, we showed the model-predicted concentration-response curves against experimental data, see Figure 4. The experimental data suggested some compounds can cause neuronal degeneration like rotenone but not carboxine. Such compound-specfic contract was captured in our model predictions.
OneHealthStudio.
Concentration-response prediction for HepG2
For both rotenone and deguline, the model could capture well with all concentration-response curves for three assays, also show in lumped plot in panel B of Figure 3. In model-based extrapolation, in higher concentration like 100 \(\mu M\) of rotenone (Figure 5) could also greater reduction in viability than that deguelin (Figure 6).
Concentration-response prediction for RPTEC/TERT1
For rotenone and degulin, we observed significantly large variability in the readout like viability, see Figure 7 for rotenone. The overall simulation matched the observation, but we also suspected wonder the cause cell death for RPTEC/TERT1 cells, given the large variability readout in the viability assay.
4 Discussion and perspective
This study adapted our previous quantitative systems model of oxidative phosphorylation (H. Yang et al. 2021) and further incorporated later key events like impaired proteostasis and enhanced cell death or neuronal degradation. With a coherent model structure, our systems models can capture experimental concentration response relationship measured from various assays from three different cell-lines, namely LUHMES, HepG2, and RPTEC/TERT1. This pilot study demonstrates that systems modeling of adverse outcome pathways can provide chemical-agnostic hazard prediction for 10 different compounds and for three cells from three tissues (neuronal, liver and kidney). We have been processing the parameter estimates from calibration model with three cell types. As a next step, we aim to provide an further refined model with robust estimates.
In terms of mitochondrial perturbation, one compound could inhibit multiple complex activity or act as uncouplers see example compounds in (Stel, Carta, et al. 2020). For such compounds, one can potentially test mechanism-based systems models from this study as well as previous models Y. Yang et al. (2015). On the other side, to evaluate hazard of mixture of compounds, mechanistic modeling could also provide a testable soluation when multiple MIEs are activated by a mixture of chemical stressors.
Xia et al. (2018) pointed out co-activation of mitochondrial toxicities and (anti-)oxidative stress response using a broader panel of assays from a rat-specific testing battery. As noted，The AOP:273 as well as (Tebby et al. 2022) have scoped with key events like oxidative stress with increase of reactive oxygen species (ROS). Further extension with existing sytems models of stress pathways could gain more compherhensive understanding of species/tissue specific adversity (Kuijper et al. 2017; Huan Yang et al. 2020; Zhang et al. 2010; Burgers et al. 2025).
Another way to extend quantitative systems modeling of AOP can also consider the possible distinct temporal characterization in key event relationship, as variability of time scales could present in different species as well as different tissue-specific cell(lines) (Shamir et al. 2016). Omics data (Meijer et al. 2025) and other High-throughput screening data (Xia et al. 2018; Rezaei Adariani et al. 2024) can provide valuable information about key event activation, and they should be considered in such extension of quantitative systems modeling. The emerging technique like cell painting at single-mitochondrion/cell resolution would provide more mechanistic and phenotypic information about chemical insult to mitochondria. As an essential step to model the data, machine-learning-aided imaging analysis approaches are important to categorize mitochondria into possible different phenotypes (Stel et al. 2023) . In practice. key event activities are often measured at different time points, our systems modeling approach can provide further flexibility to integrate data measured at different time points. Our pilot work above indeed incorporate the key event measurements: Complex I activity at 15-30 minutes while cell viability and membrane potential measured at 24 hours. In further incorporating other key event measurements, our systems modeling can provide a basis template for flexible and coherent extension and data integration.
As the mechanistic systems modeling can have better chance to model the pharmaco/toxicodynamics of the compounds as well as easy integration with pharmacokinetics modeling, one promising integration can be using whole-body physiologically-based kinetics modeling to extrapolate external exposure to internal exposure (i.e. tissue-specific target-sites). Previous systems modeling study (H. Yang et al. 2021) suggested complex interplay between pharmacokinetics and pharmacodynamics, linking possible target-site exposure to more realistic hazard characterization with concentration-response assessment. We could leverage such existing in-vitro and in-silico resource in future study.
Another note is on the learning paths/opportunity from cross-sectors on mitochondrial toxicity. The original design/purposes of the quantitative systems modeling focus on agrochemical like rotenone and other pesticides with parkinson liability. As a knowledge sharing tool, adverse outcome pathways bring different mechanistic understanding of mitochondrial from other chemicals like DILI drugs (Vlasveld et al. 2024), same class of drug like (Tolosa et al. 2015) and general chemicals (Xia et al. 2018).