Identifying Site-Specific Superoxide and Hydrogen Peroxide Production Rates From the Mitochondrial Electron Transport System Using a Computational Strategy

Abstract Mitochondrial reactive oxygen species (ROS) play important roles in cellular signaling; however, certain pathological conditions such as ischemia/reperfusion (I/R) injury disrupt ROS homeostasis and contribute to cell death. A major impediment to developing therapeutic measures against oxidative stress-induced cellular damage is the lack of a quantitative framework to identify the specific sources and regulatory mechanisms of mitochondrial ROS production. We developed a thermodynamically consistent, mass-and-charge balanced, kinetic model of mitochondrial ROS homeostasis focused on redox sites of electron transport chain complexes I, II, and III. The model was calibrated and corroborated using comprehensive data sets relevant to ROS homeostasis. The model predicts that complex I ROS production dominates other sources under conditions favoring a high membrane potential with elevated nicotinamide adenine dinucleotide (NADH) and ubiquinol (QH2) levels. In general, complex I contributes to significant levels of ROS production under pathological conditions, while complexes II and III are responsible for basal levels of ROS production, especially when QH2 levels are elevated. The model also reveals that hydrogen peroxide production by complex I underlies the non-linear relationship between ROS emission and O2 at low O2 concentrations. Lastly, the model highlights the need to quantify scavenging system activity under different conditions to establish a complete picture of mitochondrial ROS homeostasis. In summary, we describe the individual contributions of the electron transport system complex redox sites to total ROS emission in mitochondria respiring under various combinations of NADH- and Q-linked respiratory fuels under varying workloads.


Introduction
Reactive oxygen species (ROS) are by-products of cellular respiration, but they are now recognized as important signaling molecules. [1][2][3][4] Under pathological conditions such as ischemia/reperfusion (I/R) injury, elevated ROS levels contribute to cell death. 5 This type of injury occurs after a period of partial or complete loss of tissue blood flow (ischemia) followed by the restoration of normal flow (reperfusion). Although reperfusion is necessary to salvage ischemic tissue, it leads to significant tissue damage in a ROS-dependent manner. 5 This unfortunate and common event is also associated with clinical situations such as organ transplantation 6 and hypovolemic shock. 7 Though I/R injury can affect any organ, metabolically active tissue such as brain and myocardium are the most sensitive to this injury. Despite enormous efforts to ameliorate the detrimental effects of reperfusion-dependent oxidative stress, pharmacological interventions have not produced effective clinical treatment options. 5,8,9 A more complete understanding of redox site-specific ROS emission in I/R injury may lead to the design of efficient pharmacologic interventions.
In most mammalian cells, mitochondria are the primary source of ROS. [10][11][12] As such, mitochondrial ROS homeostasis is essential to maintaining mitochondrial and cellular physiology. In essence, ROS production and elimination act as countering forces which determine the mitochondrial net ROS emission. The mechanisms underlying ROS production have been reviewed extensively elsewhere, 13,14 and only key details are provided here. The electron transport system (ETS) complexes I and III are accepted as the major mitochondrial ROS producers; however, their respective contributions to total ROS emission vary according to the bioenergetic state of the organelle. 10 Matrix and non-ETS inner membrane proteins can also produce ROS, but their contribution is limited. [15][16][17] Under physiological conditions, electron flux through the ETS favors forward electron transport (FET). In this operating mode, pyruvate (P) from glucose metabolism is oxidized through the Krebs cycle to generate nicotinamide adenine dinucleotide (NADH), which enters the ETS at complex I. When succinate (S) is oxidized, electrons enter the ETS at complex II. Despite the different entry points, electrons converge at the quinone (Q) pool before they are passed onto cytochrome c and ultimately molecular oxygen (O 2 ). This is the normal flow of electrons during energy transduction. During certain pathological states, ETS circuitry is rewired due to a change in substrate availability. For example, during anoxia, the forward electron flow through complex I is supported by reversal of complex II while the activities of complexes III and IV dramatically fall near zero. 18,19 Reorganization of electron flux causes succinate and ubiquinol (QH 2 ) to accumulate as fumarate reduction at complex II replaces O 2 reduction at complex IV. Upon reperfusion, O 2 is reintroduced into the system, and electron flow normalizes. With high levels of QH 2 , thermodynamic and kinetic processes favor ETS operation in the reverse electron transport (RET) mode. In this mode, complex I enters a near equilibrium state with nearly all its redox centers fully reduced and produces ROS at significantly elevated rates. Thus, ROS production is enhanced under RET compared to FET with the former playing a significant role in I/R-induced oxidative stress.
The majority of O 2 is fully reduced to water at complex IV; however, one-or two-electron reduction of O 2 by an ETS redox center upstream of complex IV produces superoxide or hydrogen peroxide, respectively. In complex I, the redox centers include a flavin mononucleotide (FMN) at the NADH oxidase site (site I F ), 20 a semiquinone (SQ) at the Q reductase site (site I Q ) 21 and a chain of iron-sulfur (Fe-S) clusters that rapidly relay electrons between the two. 22,23 Likewise, complex II harbors a flavin adenine dinucleotide (FAD) at the succinate oxidase site (II F ), a Q reductase site (II Q ), and a chain of Fe-S clusters. 24,25 The mammalian complex III redox centers include cytochrome c 1 , a highand low-potential b type heme, the Rieske iron-sulfur protein (ISP), and two quinone binding sites (Q N and Q P ). 26 Superoxide arising from complex III is generally accepted to originate from the SQ at the quinone binding site proximal to the intermembrane space (Q p ). 27 The topic of site-specific ROS production has attracted experimental researchers for many years. Both pharmacological and genetic approaches have been utilized to dissect the contribution of redox centers to total mitochondrial ROS. While complexes I and III are the major ROS producers, the contribution of specific redox centers at physiologically relevant metabolic states remains elusive. For example, both site I F [28][29][30] and I Q 31 are proposed as major ROS-producing sites of complex I. Complex II was widely accepted as a negligible source of ETS ROS. 32 However, several more recent studies report that complex II produces significant amounts of ROS under appropriate conditions. 33,34 Such discrepancies in experimental data arise from the caveats inherent to experimental studies and have contributed to unsuccessful efforts to target oxidative stress in clinical settings. In particular, the use of chemical inhibitors and genetic models inevitably alter the native electron distribution in unpredictable patterns, leading to epiphenomena and conflicting experimental results. Differences in experimental conditions including species, tissues, developmental stage, etc. introduce additional confounding factors that are often neglected and make comparing results from different studies difficult. Lastly, the pursuit is further impeded by the lack of a robust method to distinguish species-specific ROS without interfering with other mitochondrial and extra-mitochondrial processes.
Computational modeling is a useful platform to address these challenges. A computational model serves as a quantitative framework to analyze mitochondrial bioenergetic data from a single, unified perspective. Several models of mitochondria exist at varying levels of complexity, focus, and approach. [35][36][37][38][39][40][41][42] Some investigate the kinetics underlying ROS production by individual ETS complexes [38][39][40][43][44][45] while others integrate mitochondrial ROS production and elimination in a substratespecific context. 41,42,46 However, to our knowledge, none consistently reproduces a wide range of experimental data using a unified, coherent, and consistent framework. Thus, the origin of superoxide and hydrogen peroxide from the ETS remains elusive from a computational or experimental perspective.
We herein developed, analyzed, and corroborated a comprehensive model that simulates mitochondrial ROS homeostasis under a range of bioenergetic conditions. This model is focused on ROS from the ETS and is referred to as the ETS-ROS model. Modules relevant to ROS originating from the ETS were individually developed, parameterized, and corroborated against a variety of data sets. [38][39][40] The model includes the primary ROS producers in mitochondria with their biochemical reactions constrained thermodynamically and balanced with respect to mass and charge. These modules were then incorporated into a single integrated framework, producing a model that is coherent and operates within biophysical constraints using a single, consistent set of parameters. The ETS-ROS model reproduces experimental data of not only net ROS production but also other mitochondrial bioenergetic variables. Using this approach, we have identified the regulation of ROS production by substrate metabolism via alterations to the mitochondrial membrane potential along with the NAD and Q pool redox states. The model also reveals that the scavenging system is saturable with a K M for H 2 O 2 in the nM range. Therefore, under pathological conditions such as I/R injury, cellular oxidative stress is a result of ROS production overwhelming ROS elimination.

General Approach to Modeling and Processes Included
Our modeling approach is modular and contains enzyme kinetic and ROS production models of complexes I, 38 II, 39 and III 40 previously published by our group. These models are biophysically detailed and thermodynamically consistent. They were individually calibrated and corroborated against a wide range of kinetic data sets. In this ETS-ROS model, we opted to explicitly model part of the Krebs cycle while lumping other parts. The explicit parts of the model consist of the enzymatic reactions that support succinate, fumarate, and malate transport and oxidation. Pyruvate transport is explicitly modeled; however, pyruvate dehydrogenase and the upper branch of the Krebs cycle (citrate synthase through succinyl-CoA synthetase) were modeled as a lumped reaction. This lumped reaction was sufficient to explain the experimental data. The detailed reactions for the lumped portion of the model are excluded to stay faithful to the parsimonious principle. In addition, our in-house experimental data reveal a moderate level of malic enzyme (ME) activity in mitochondria isolated from guinea pig ventricular cardiomyocytes ( Figure S1). Several other studies also support the expression of ME in cardiac tissue of guinea pigs, 47,48 and we found this reaction necessary to fit the model to the calibration data sets. Other biochemical processes such as spontaneous oxaloacetate decarboxylation 49 may also help explain these data. We will test this and other ideas when the model is updated to include more detailed Krebs cycle reactions. A summary of processes that are explicitly modeled are shown in Figure 1. Notations of the redox centers included in the text are described in Table 1.

Model Simulations
The model was numerically simulated using MATLAB (R2020b). The parameter optimization was performed on a Dell desktop PC (64-bit operating system and x64-based processor Intel R core TM i7-7700 CPU @3.60GHz and 16 GB RAM) using the Parallel Computing Toolbox. A parallelized simulated annealing algorithm was first used to globally search for feasible parameters, which were then refined using a local, gradient-based optimization algorithm. When fitting to the calibration data, the difference between model output and data were normalized with the respetive uncertainty in the data and minimzed using the standard least squares method.

Data Sets
The intricate relationship between substrate utilization and ROS homeostasis necessitates that model calibration includes experimental data describing both NADH-and QH 2 -linked substrate supported respiration and ROS production. The calibration data set includes the oxygen consumption rates (J O2 , nmol/mg/min), net hydrogen peroxide emission rates (J H2O2 , pmol/mg/min) and NADH reduction (%) under NADH-and QH 2 -linked supported respiration in both FET and RET modes. FET occurs when respiration is supported by pyruvate/L-malate (P/M) and S in the presence of rotenone (S/R). RET occurs under QH 2 -linked supported respiration when rotenone is absent (S). The calibrated model is subsequently challenged against novel data sets, those which were not used for model calibration. The first corroboration data set includes the J O2 and J H2O2 data obtained when respiration is supported by both NADH-and QH 2 -linked substrates (P/M/S). The second data set recapitulates the monotonic relationship of J H2O2 on oxygen concentrations in FET mode. The list of all the data sets utilized in the construction, calibration and corroboration of the model is summarized in Table 2.
The majority of the experimental data used for model calibration were from a previous publication of ours. 50 In this study, we measured the oxygen consumption rate (J O2 , nmol mg -1 min -1 ) and net hydrogen peroxide emission rate (J H2O2 , pmol mg -1 min -1 ) of mitochondria respiring on P/M, S and S/R. Briefly, all experiments were performed using a standard KClbased mitochondrial suspension buffer containing 130 mM KCl, 5 mM K 2 HPO 4 , 20 mM MOPS, 1 mM MgCl 2 , 1 mM EGTA, and 0.1% w/v BSA (pH of 7.1 at 37 • C). For general methods and protocols regarding mitochondrial isolation and bioenergetic experiments, see our prior works. 50,51 When used, the final concentrations of mitochondria, substrates, and inhibitors were 0.1 mg/mL mitochondria, 5 mM pyruvate (P), 1 mM L-malate (M), 10 mM succinate (S), 1 μM rotenone (R), and 500 μM adenosine diphosphate (ADP). The data were analyzed using the software MATLAB R .

Model Parameters
The fixed model parameters are obtained from thermodynamic and kinetic data obtained from prior work. 18,[38][39][40] The adjustable parameter set consists of 14 parameters related to both ROS kinetics and substrate metabolism. We chose the smallest set of adjustable parameters necessary to simulate the calibration data sets. The ROS kinetic parameters are related to mitochondrial production and clearance of superoxide and hydrogen peroxide. The capacity of mitochondria to remove hydrogen peroxide is interchangeably referred to as the scavenging capacity in this paper. The scavenging system is modeled in a Michaelis-Menten like fashion-saturable and dependent on the hydrogen peroxide concentration. The ETS parameters include the enzyme contents for complex I, II, and III, the estimation of the semiquinone stability constant in complex III, and the apparent rotenone binding constant to complex I. The remaining adjustable parameters are concerned with metabolite transport and substrate utilization (Table 3). Normalized local sen-  For details regarding the approximate location of each site on the complexes, see Figure 1. Site locations are not 100% accurate and drawn for clarity. Subscript indicates the form of the redox cofactor responsible for ROS production. In complex II, site [3Fe-4S] is located at the quinone binding site and thus referred to as IIQ in Figure 1.
sitivity coefficients are computed to inform the contribution of each adjustable parameter to the model output in a similar manner previously described. [38][39][40] For details on how the analysis was performed, see Section 3 of the supplement. These values are calculated from model outputs coincident with the data and associated dynamics. The sensitivity analysis reveals that the complex I content is the most sensitive parameter followed by complex III content and the H 2 O 2 affinity of the matrix scavenging system. These results are not unexpected considering the tight association between respiratory state and ROS emission.

Isolated Mitochondria Preparation
We supplemented prior data of ours 50 with additional data to further challenge the model. For these studies, we used Hartley guinea pigs weighing between 350 and 450 g (4-6 wk). Animal care and handling conformed to the National Institutes of Health's Guild for the Care and Use of Laboratory Animals and was approved by Michigan State University's Institutional Animal Care and Use Committee. Isolated mitochondria were isolated from ventricular cardiomyocytes as described in our prior work. 51,52 In brief, the guinea pigs were euthanized after being anesthetized with 5% isoflurane and tested to be unresponsive to noxious stimuli. The heart was perfused in chest with icecold cardioplegic solution until no blood was observed in the coronary arteries and cardiac veins. The cardioplegic solution consisted of 25 mM KCl, 100 mM NaCl, 10 mM dextrose, 25 mM MOPS, and 1 mM EGTA at pH = 7.15. The heart was then excised and washed with ice-cold isolation buffer consisting of 200 mM mannitol, 50 mM sucrose, 5 mM K 2 HPO4, and 0.1% w/v bovine serum albumin (BSA) at pH = 7.15. Ventricular tissue was minced in ice-cold isolation buffer until small pieces remained (∼mm 3 ).
The homogenate was then transferred to a 50 mL canonical tube containing 0.5 U/mL protease (Bacillus licheniformis) in 25 mL isolation buffer. A handheld Omni tissue homogenizer was run at 18 000 rpm for 20 sec to liberate mitochondria from cellular components. Mitochondria were then recovered by gradient centrifugation in isolation buffer at 4 • C. Mitochondrial protein was quantified using the bicinchoninic acid assay and BioTek Synergy H1 plate reader. The mitochondrial stock solution was diluted to a working concentration of 40 mg/mL. We routinely run quality control steps to test for mitochondrial quality that include respiratory control ratios (RCRs) (maximum oxphosstate rate/leak-state rate), and we only use mitochondria with control ratios above 16. All reagents were from Sigma unless otherwise noted.

Oxygen Consumption and Net Hydrogen Peroxide Emission Rates
The Oroboros Oxygraph (O2k) was used to simultaneously measure J O2 (nmol mg -1 min -1 ) and J H2O2 (pmol mg -1 min -1 ), described in detail in our previous publications. 50 Briefly, J H2O2 values were measured using the Amplex UltraRed/Horseradish Peroxidase/Superoxide Dismutase assay. The reduction of hydrogen peroxide is coupled with the oxidation of Amplex Ultra-Red to resorufin. Resorufin fluorescence is converted to hydrogen peroxide concentration by using a hydrogen peroxide calibration curve obtained each experiment day from freshly prepared hydrogen peroxide standards. A post-hoc analysis of our previous J H2O2 data sets using S as the substrate reveals that Sdependent J H2O2 is highly sensitive to the RCR. Higher J H2O2 values correspond to RCR (energized with P/M and in the presence of 4 μM Ca 2+ ) greater than 18. Thus, we repeated J H2O2 measurements using mitochondria with an RCR above 18. Mitochondria were energized with L-M, S or the P/M/S combination (leak state). An ADP bolus was added to stimulate oxidative phosphorylation (oxphos state) after 5 min under L-M supported respiration, 2.5 min under S-and 2 min under P/M/S-supported respiration. These time periods were chosen to allow measurement of steady-state J O2 while maintaining adequate oxygen supply for oxidative phosphorylation turnover.

NADH Fluorescence Measurement
The fluorescence of NADH was monitored on an Olis DM-245 spectrofluorometer. Fluorescence was monitored with λ excitation = 355 nm (8 nm bandpass filter) and λ emission = 450 nm (13 nm bandpass filter). Baseline fluorescence was monitored for 1 min followed by mitochondria addition. Substrates were added after 5 min of equilibration. A bolus of ADP stimulated oxidative phosphorylation. Fluorescence minima were obtained by uncoupling mitochondria with trifluoromethoxy carbonylcyanide phenylhydrazone (FCCP) to produce maximally oxidized nicotinamide pools. Maximal fluorescence was obtained in the presence of rotenone to produce maximally reduced nicotinamide pools.

Model Calibration
The experimental fidelity of the ETS-ROS model is critically related to its ability to predict the in vivo response of the ETS under a wide variety of physiological and pathophysiological scenarios. Including both FET and RET modes of ROS production data into model calibration was essential to identify the parameters associated with site-specific ROS fluxes under these diverse experimental conditions. In the absence of exogenous ADP, the substrate combinations of P/M and S/R induce FET whereas RET was favored with S as the sole substrate. We used the P/M/S data, which consist of both FET and RET processes, to corroborate the model.

Oxygen Consumption and Net Hydrogen Peroxide Emission Rates
Outputs from the calibrated model were compared to experimentally determined J O2 (nmol/mg/min) and J H2O2 (pmol/mg/min), as shown in Figure 2 and summarized in Table 4.
In general, the model captures the dynamic profile of O 2 consumption as well as H 2 O 2 emission rates in a range of energy demand conditions. In the leak state, respiration rates are low and H 2 O 2 emission rates are high. This is general for all substrate conditions and is due to membrane polarization and highly reduced mitochondrial redox couples (eg, high NADH/NAD + ). When succinate is present, the QH 2 /Q ratio is large and H 2 O 2 emission rates reach even higher values. The highly reduced redox centers in the ETS proteins thermodynamically and kinetically favor elevated ROS emission rates. When ADP is added, mitochondria enter the oxphos state and dramatically increase their respiratory rate. In this state, redox couples become more oxidized and thus produce ROS at a lower rate compared to leak state. This relationship is also consistent across substrate conditions. Relative to the S/R condition, the absence of rotenone for the S condition did not affect leakstate J O2 but slightly depressed oxphos-state J O2 . This results from NAD pool oxidation leading to oxaloacetate accumulation and complex II inhibition due to malate dehydrogenase turnover. 50,53 Importantly, the H 2 O 2 emission rate is dependent on O 2 levels (will be discussed in detail below and briefly here). In these simulations, we modeled the buffer O 2 dynamics as we will show later. Our initial buffer O 2 concentration was 188 μM in a 2 mL volume of uniformly dispersed isolated mitochondria. In our conditions, a respiratory rate of 100 nmol O 2 /mg/min corresponds a rate of 10 μM/min which will deplete all the O 2 in 18.8 min. These simulation results show the model not only reproduces the respiratory dynamics associated with various substrates and mitochondrial work rates, but it also simulates  Experimental data are obtained using mitochondria isolated from the ventricular cardiomyocytes of guinea pigs, as described in our previous work. 50 As shown, the model outputs are within the experimental ranges for both JO2 and JH2O2 in both FET and RET modes.

NADH Reduction States
An essential bioenergetic variable associated with H 2 O 2 emission is the NAD redox state. The changes to the NAD + pool redox state vary with substrate delivery, ETS flux, and phosphorylation potential, reflecting the overall redox balance of the mitochondrial matrix. As shown in Figure 3, the model is capable of faithfully reproducing experimentally observed NADH dynamics of mitochondria respiring in the leak and oxphos states with P/M or S/R as substrates. Succinate-supported respiration with rotenone yielded a more reduced NADH pool in both leak and oxphos states compared to the P/M substrates. With P/M as the substrates, turnover of NADH at complex I favors a more oxidized NADH pool in leak (80%) and oxphos (30%) states relative to S/R (87.5% and 90%, respectively). The rapid NADH reduction dynamics seen with S/R are a result of residual ME and site I Q activity. We used 1 μM rotenone in these experiments, and the reported apparent dissociation constant of rotenone binding at the I Q site is 20 nM. 54 In that study, the I 50 value was determined using submitochondrial particles with a CoQ 10 analogue. To properly identify the rotenone dissociation constant, kinetic and free radical studies using intact mitochondria with the native Q molecule used by complex I are necessary. Thus, we opted to fit this parameter with our calibration data sets. The resultant model simulations predict that complex I operates at 2% of its maximum capacity. That said, assuming simple inhibitor kinetics and the above stated I 50 value for rotenone, residual complex I activity is expected to be ∼2% which fits with our model simulations. With the J O2 , J H2O2 , and now the NADH calibrations, we next sought to challenge our model with an additional substrate combination without adjusting any model parameter other than setting the initial conditions (substrate concentrations).

Oxygen Consumption, Net Hydrogen Peroxide Production Rates, and Oxygen Concentration Dynamics Under Both NADH-and QH 2 -linked Substrate Metabolism
The final adjustable parameter set ( Table 2) was used to simulate the model and predict the J O2 , J H2O2 , and [O 2 ] dynamics when mitochondria were energized with a combination of NADH-and QH 2 -linked substrates (P/M/S). Under this condition, electrons entering the ETS at complexes I and II converge at the Q pool. 55 Such conditions favor the reduction of the NAD + pool, the Q pool, and a high membrane potential, resulting in ROS generation at complex I through RET. As shown in Figure 4, the ETS-ROS model outputs are consistent with experimental data collected from the P/M/S condition. Experimentally determined leak-state J O2 (nmol/mg/min) are 129 ± 5 at 2 min and 121 ± 12 at 4 min. The oxphos-state J O2 is 729 ± 91 nmol/mg/min. The corresponding simulated J O2 (nmol/mg/min) are 116 at 2min leak state, 114 at 4-min leak state, and 654 during oxphos ( Figure 4A). Experimentally measured J H2O2 (pmol/mg/min) are 4304 ± 499 at 2-min leak state, 3843 ± 713 at 4-min leak state, and 82 ± 7 during oxphos. The corresponding model predictions are 4449, 3100, and 42 pmol/mg/min ( Figure 4B). In general, J O2 and J H2O2 were significantly elevated compared to the calibration substrate mixtures, which favor electron entry at either complex I or II. With the P/M/S substrate combination, both complex I and complex II are fully engaged to deliver electrons into the Q pool during oxphos. This will increase the respiration rate with the high substrate delivery favoring reduction of mitochondrial redox couples. The model predicts the resting and active respiratory states with very good accuracy. The small underestimation of ROS emission in the oxphos state is due to minor sources we did not include in this model.

Monotonic Relationship Between J H2O2 on [O 2 ]
Although J O2 is essentially independent of [O 2 ] until it drops below 1 μM, the monotonic relationship between J H2O2 and [O 2 ] has been described by several groups. 50 Figure 5B). The rapid drop in hydrogen peroxide production as [O 2 ] approaches zero is due to a decreased availability of the reduced I FMN . In our model, flavin sites must be unoccupied before they are able to react with O 2 . This occurs because of the following: first, the reduced I FMN preferentially binds to NAD + over NADH. 38,58 Second, the kinetic pressures between NADH production (via residual ME activity) and NADH consumption (via ROS production reactions) results in a gradual but complete reduction of the NAD pool as the [O 2 ] required for ROS production approaches zero. Consequently, the net J H2O2 vs [O 2 ] is non-linear with the S/R combination which we will now explicitly show in the next section. A similar mechanism likely exists for the P/M condition, but our simplified, lumped Krebs cycle expression does not capture this phenomenon well. We anticipate that we will reproduce this relationship when the model is extended to include better a representation of the enzymatic processes in the Krebs cycle.

Model Predictions
The high level of consistency between model simulations and experimental data demonstrate that our model is structurally sound and appropriately calibrated. We next used the ETS-ROS model to (1) investigate the linear and non-linear relationship of [O 2 ] on J H2O2 , (2) predict the effects of substrate utilization on site-specific and species-specific contributions to mitochondrial ROS production during different respiratory states and ETS flux directions, (3) assess the degree of mitochondrial scavenging activity and % utilization during the experimental conditions in the calibration and corroboration data sets, and (4) elucidate the impact on the rate of ROS emission of the redox states of the NAD and Q pools, along with the mitochondrial membrane potential. These model predictions provide valuable mechanistic information about mitochondrial redox homeostasis that are experimentally inaccessible.

Hydrogen Peroxide Production By Site I F Underlies the Kinetics of Net ROS Emission Rates At Low [O 2 ]
Preliminary model analysis of the monotonic dependence between net J H2O2 and [O 2 ] reveals that, when mitochondria are energized by succinate and RET is inhibited with rotenone, H 2 O 2 production is non-linearly sensitive to changes in [O 2 ] ( Figure 5B). We previously explained this phenomenon with an   empirical model which was unable to identify redox center contributions to this behavior. 50 Using our computational ETS-ROS model, we dissected the individual processes underlying the net J H2O2 measurements as shown in Figure 6. Here, the total ROS (J ROS ) is defined with respect to 2-electron equivalents and equals the sum of hydrogen peroxide (J H2O2 ) plus half of the superoxide flux (J SO /2). The scavenging capacity shown in right y-axis of panels shown in Figure 6A and 6D reflect the fraction scavenging capacity utilization. In the P/M condition, the scavenging system activity is only a few % of the maximum rate of H 2 O 2 clearance. In contrast, the scavenging system activity increases to 15% when mitochondria are energized with S/R. This increase in scavenging activity is expected as mitochondria respiring in the S/R condition produce more ROS compared to the P/M condition ( Figure 6A and 6D). This goes up even further when S or P/M/S are used as respiratory fuels, as we will show later. In the presence of P/M, a linear dependence between total J ROS and [O 2 ] was observed, and H 2 O 2 production is low ( Figure 6A). However, J ROS exhibit a non-linear dependence on [O 2 ] under S/R-supported respiration ( Figure 6D). The primary reason for this difference is the H 2 O 2 producing activity of the I F site.
Regardless of the substrates, most H 2 O 2 originates from the fully reduced flavin of complex I, which produces the non-linear The model further demonstrates that the redox states of the NADH and the Q pools are highly sensitive to experimental conditions (C and F). Rotenone inhibits QH2 oxidation at site IQ, which indirectly prevents effective NADH oxidation at site IF. The NADH pool is, thus, significantly reduced in the S/R simulation, creating a condition conducive for H2O2 production. relationship between J ROS and [O 2 ] for S/R supported respiration ( Figure 6B and 6E). The combination of highly reduced NAD + and Q pools redox states explains the J H2O2 profile associated with each substrate combinations ( Figure 6C and 6F). In the P/M simulation, the NADH pool is less reduced, and the Q pool is significantly oxidized. In contrast, rotenone in the S/R simulation prevents quinol oxidation at site I Q and NADH oxidation at site I F , maintaining the NADH and Q pools in mostly reduced states. As seen in Figure 6F, the NADH level at high [O 2 ] is 90% reduced and rises to 100% as [O 2 ] approaches zero. At high [O 2 ], most of the reduced flavin is bound up with NAD + despite it being only 10% of the NAD + pool. This is because NAD + binds with about 10-fold higher affinity than NADH to the reduced flavin on complex I. However, as O 2 is depleted, the NAD + pool becomes progressively reduced freeing the binding site near the reduced flavin to maintain H 2 O 2 production at elevated rates until anoxia ( Figure 5B, inset).

Site-Specific Superoxide and Hydrogen Peroxide Production Rates
As shown in Figures 7 and 8, the site-specific contributions to total superoxide and hydrogen peroxide production are determined by substrate combinations, respiratory states, and ETS flux directions. The site-specific superoxide production rates are also shown in Table 5 whereas site-specific hydrogen peroxide production rates are included in the text below.
During leak state with substrate combinations favoring RET (S and P/M/S), ROS is produced primarily at sites I Q and I FMNH2 ( Figure 7A). The rates from these sites exceed 25 times that of the next highest site. When operating in FET mode, these two sites are joined by III Qp as the top 3 ROS producers. In the S/R condition, sites I Q and I FMNH2 are replaced by both CII sites as top ROS producers when their individual contributions are combined. Even so, I FMNH2 still produces ROS in this condition because the FMN on complex I is unobstructed in the presence of rotenone. Both complex II sites increase their individual ROS production rates due to the presence of a highly reduced Q pool. During oxphos, bioenergetic conditions are quite different and the mitochondrial redox couples become more oxidized, in general. Figure 7B shows that the ROS production rates in the oxphos state from all sites are dramatically reduced with a similar profile relative to the leak state. Notably, sites I Q and I FMNH2 remain the sites with the highest capacity to produce ROS. Despite high ETS turnover, complex I redox centers are maintained in a reduced state, producing high ROS emission rates. In the oxphos state, the NAD + pool becomes more oxidized, but the Q pool gets more reduced. This is due to the activation of parts of the tricarboxylic acid (TCA) cycle via metabolite and allosteric modulatory changes. 18 As a result, complex II and III sites do not significantly drop below leak-state rates despite a large drop in membrane potential and NADH levels. The high rates of sites I Q and I FMNH2 in S and P/M/S conditions during the oxphos state is somewhat surprising. During leak state, the FMN of complex I is operating in a near-equilibrium state under highly reducing conditions. The I Q site feeds electrons into the complex to support ROS production at both sites. In contrast, the sites at complex I are displaced from equilibrium during oxphos and contribute to electron flow down to O 2 . Yet, in this state, these sites still produce ROS at a rate of around 5 times or more when compared to the other sites. The sites that produce H 2 O 2 show a similar pattern.
As direct H 2 O 2 production by the respiratory chain is a 2electron redox reaction, only the I FMNH2 and II FADH2 are considered as ETS sources of H 2 O 2 in the model. Regardless of respiratory states and substrate sources, almost all H 2 O 2 produced by mitochondria originate from the I FMNH2 site (Figure 8). Site II FADH2 produces very little H 2 O 2 in the presence of fumarate. 39 Transitioning to oxphos from leak state leads to a more dramatic reduction J H2O2 relative to J SO . This is because J H2O2 relies more heavily on NADH rather than QH 2 or the membrane potential. 38,39 In the S/R condition, the mitochondrial NAD pool is more reduced relative to the P/M condition due to the (i) presence of ME activity and (ii) inhibition of complex I. Under S and P/M/S conditions, the model predicts that the hydrogen peroxide production rate from site I FMNH2 is similar. Consequently, most of the decrease in total hydrogen peroxide after the transition to oxphos results from decreased production by this site. However, phosphorylating mitochondria energized by P/M/S produce   more hydrogen peroxide from the I FMNH2 site compared to S alone. In the P/M/S condition, pyruvate metabolism produces NADH through the TCA cycle which results in a more reduced NAD + pool. Consequently, the I FMNH2 site is more reduced under P/M/S-supported respiration, enhancing hydrogen peroxide arising from this site.

Effects of Substrate Utilization and Electron Transport Mode on Scavenging Activity
The activity of the scavenging systems' response to changes in ROS production remains an enigma for the same reasons that the origins of mitochondrial ROS under native conditions are inconclusive. Therefore, due to the lack of quantitative data, we modeled the scavenging system assuming that it is saturable and responsive to changes in matrix H 2 O 2 production. To do so, we used a simple, lumped empirical Michaelis-Menten expression shown in Section 3 of the supplement. This expression consists of a V max and a K M term for H2O2 affinity which were fit to the data presented herein. Because regenerating the scavenging system relies on NAD(P)H, ROS scavenging has a limited capacity. 10 Thus, the responsive nature of the scavenging system to matrix H 2 O 2 is necessary to keep net ROS production in the physiological range. In this model, these assumptions were adequate as the model outputs in the calibration and corroboration stages are consistent with experimental data. Thus, we used the model to predict the rate of ROS scavenging in the experimental conditions simulated in Figure 1. The model reveals that scavenging activity is lower in FET mode, which is unsurprising since it is consistent with lower ROS production. In RET mode, the scavenging activity is operating near maximum and ROS emission rates increase nearly 40-fold. Currently, the kinetic properties of the intact scavenging system remain a deep mystery. This is, in part, due to the complex biochemistry occurring in the mitochondrial matrix. For example, there is compelling evidence for a ROS-related substrate channeling system embedded in the matrix, 59 but the molecular details of this phenomena are still obscure. Thus, changes in scavenging capacity or H 2 O 2 sensitivity can result in elevated ROS emission rates that will lead to cellular oxidative stress.

Primary Determinants of ROS Production
The corroborated model is further used to explore the effects of the mitochondrial membrane potential along with the NADH and Q pool redox state on J SO and J H2O2 from individual complexes ( Figure 9). In this simulation, we simulated the model until reaching the steady state with respect to mitochondrial state variables under the P/M condition. We clamped substrates at saturating levels and oxygen at a physiological level of 20 μM 10 in this simulation. We then perturbed the membrane potential, NAD + redox state, and Q pool redox state while keeping the other state variables fixed and computed the resultant free radical fluxes. This simulation exemplifies the utility of model in testing hypotheses that are not experimentally feasible. The model predicts a maximum rate of ROS production occurs when the membrane potential is high and the NAD + and Q pools are extremely reduced (Figure 9). In general, at high membrane potentials, mitochondrial ROS production is more sensitive to NADH levels than QH 2 levels. When NADH levels begin to exceed 90%, ROS production exponentially increases. That said, QH 2 levels still significantly affect total ROS produciton showing saturation at more oxidized levels. Of the H 2 O 2 producing sites, complex I produces significantly more ROS compared to complex II. Its maxima lie along the 100% NADH and 220 mV membrane potential plane where QH 2 is maximally reduced. In contrast to the H 2 O 2 producings sites, the rates of the SO producing sites are distributed differently. Complex I dominates in RET mode and contributes little to SO production under more oxidized conditions. Complex II produces the least amount of SO under these conditions and is independent of NADH. As expected from our prior study, complex II SO production is very sensitive to QH 2 levels. In contrast, complex III produces a moderate amount of SO and contributes to the basal level of SO production by the mitochondrial respiratory chain. As with the other ETS pumps, this complex produces most SO when the membrane potential, NADH levels, and QH 2 levels are high.
The ETS-ROS model was used to predict the bioenergetic response of mitochondria to an increase in adenosine triphosphate (ATP) demand for the substrate combinations use for calibration and corroboration. These results shown in Figure 10 reveal the respiratory dynamics, ROS emission profile, and dynamics of other bioenergetic variables encoded by the model follow expected behavior. With an increase in ATP demand, respiration monotonically increases until the max capacity of the ETS and oxphos-related processes are attained. At high demand, the ETS pumps are outcompeted by the adenine nucleotide translocase, inorganic phosphate carrier, and ATP synthase. Consequently, the membrane potential decreases in a Figure 10. Select model outputs demonstrating the changes in key bioenergetic variables during changes in ATP demand. Model simulation protocols were like those used for simulating the experimental conditions presented in the accompanying article. The only major differences are the inclusion of 10 mM ATP in the extramitochondrial space as an initial condition, a more physiological O2 concentration at 20 μM, and external metabolites, buffer pH, and O2 were clamped. The model was then run out to a steady state at the ATPase Vmax's as shown on the x-axis. These simulation results are consistent with what we know about the impact of ATP-demand on respiration (A), ROS output (B), membrane potential (C), and extra-mitochondrial, steady-state ADP levels (D) for the given conditions. The decrease in the respiratory rate with succinate as ATP demand increases is due to OAA inhibition. titratable manner. As a corollary, mitochondrial ROS emission is decreased. Due to the relative membrane potential independent kinetics of succinate dehydrogenase, the S/R condition results in the fastest ROS emission rate when ATP demand is high. What is also revealed by these simulations is the ability of mitochondria to defend against ATP demand when substrate supply is ample and diverse. This notion is reflected in the membrane potential and extra-mitochondrial ADP concentrations. The P/M/S combination resulted in the highest membrane potential and the lowest ADP concentration profile among all combinations. A strong dependence on the membrane potential for ROS output is reflected in the 5-fold drop of ROS emission rates when the membrane potential decreases only about 10% (∼20 mV). Figures S2 and S3 show Figure 10 results plotted against ATPase rate instead of Vmax. The plots look quantitatively different because our ATPase expression (given in Section 3 of the supplement) encodes an ADP-dependent negative feedback mechanism we developed to represent external ATPases in our experiments. 18 In addition, we explored the role succinate plays on enhancing energy production and ROS emission at sub-saturating pyruvate (100 μM) and malate concentrations (0.5 mM) in the presence of moderate ATP demand (∼50% of max). The results shown Figure S4 reveal that as succinate levels are increased, there is a corresponding increase in respiration, ROS emission, membrane potential, and % reduced redox pools (NADH, QH 2 , and cyt c) with a concomitant decrease in extra-mitochondrial ADP. The lower ADP levels are a direct consequence of increasing the membrane potential that enables the generation of stronger phosphorylation potentials (∼[ATP]/[ADP]/[Pi]). Each of these results demonstrate the model's ability to characterize bioenergetic changes associated with energy demand and support its general utility. That said, these results should be interpreted with care until the model is validated after the addition of other major matrix enzymes and transporters involved in energy metabolism.

Model Predictions of Site-Specific ROS From the ETS Complexes
Site-specific ROS production has inspired many experimental studies and led to a wealth of experimental data. However, the data are limited to experimental conditions wherein inhibitors are present. The use of different systems, experimental conditions, and techniques has also led to discrepancies among these studies. Simpler systems, such as purified enzymes, have the advantage that fewer variables need to be accounted for in interpreting experimental results. But the conclusions may not be applicable when these results are integrated into biochemical reaction networks. On the contrary, a more intact system affords an environment more similar to in vivo conditions but introduces more hidden variables. Perhaps the biggest hidden variable in all mitochondrial studies is the dynamic redox state of the Q pool. There is currently no experimental method that can be used to quantify the real-time dynamic changes of the Q pool without sample destruction. Additionally, factors that affect thermodynamics such as temperature, buffer ionic  D). ADP addition results in maximal oxidative phosphorylation (oxphos state). During oxphos, the scavenging activity decreases regardless of the electron source and transport mode. Under RET mode, the scavenging activity is significantly increased but unable to keep the net JH2O2 at comparable levels to under FET. The model predictions of active scavenging, thus, together with the changes in JH2O2 and membrane potential support the assumptions that the scavenging activity is saturable and responds to changes in ROS production. strength, and buffer pH can all regulate enzyme activities but are often neglected. These changes can lead to large changes in pathway fluxes and study outcomes. A computational approach can be applied to fully realize unified analyses at high temporal and spatial resolutions. Thus, different experimental conditions are analyzed using a unified framework to reconcile contradictory reports. This quantitative framework not only checks for the internal consistency of data but also enables extrapolation to experimentally untestable space for hypothesis testing or generation.
While detailed computer models of mitochondrial metabolism have been around since the 1970s, [60][61][62][63][64][65][66] there are only a few that integrate ROS homeostasis with substrate metabolism. 18,35,41,42,46,67 For an overview of these models, the review by Mazat et al. is an excellent source. 68 These models each have their own unique strengths and focuses, but the current ETS-ROS model is distinguished from these prior models in several important ways. First, complexes I, II, and III are included as sources of ROS while prior models only included I and/or III. As this study demonstrates, ROS from complex II is not insignificant, and is produced at high rates under the right circumstances. So, to quantitatively ascribe ROS production rates to each complex, all three major ROS producing sites are required. Second, the ETS-ROS model includes mechanisms required to simulate the dynamic range between low (leak) and high (oxphos) flux regimes. Not all prior models possess this ability without requiring parameter or structural changes. In the ETS-ROS model, the metabolite concentration profiles control electron flow through the ETS. This type of behavior is desirable when extrapolating up to more complex systems. And third, the core elements of the ETS-ROS model are coded up to ensure thermodynamic consistency. This feature is shared by some models and keeps model behavior within the feasible biochemical reaction space, and it strengthens the utility of multi-scale models. In any case, each of these models are extremely useful in aiding our ability to understand the general features of mitochondrial ROS metabolism.
Using the ETS-ROS model, we identify the ROS species and their origin (Figures 8 and 11) under the different metabolic conditions of our experimental data ( Figure 2). We found that the III Qp site, site I FMNH2 , and site I Q as a semiquinone are the primary sources of superoxide in non-phosphorylating mitochondria energized with P/M (Figure 7 and Table 5). In nonphosphorylating mitochondria, site I FMNH .
gives rise to a considerable amount of superoxide under QH 2 -supported respiration when RET is inhibited (S/R). These model predictions reconcile experimental results regarding complex I from different studies. For example, Galkin and Brandt 30 and Grivennikova and Vinogradov 29 independently concluded that both the fully reduced flavin (FMNH 2 ) and the flavin radical (FMNH .-) can give rise to superoxide. Kussmaul and Hirst also found that site I F is a primary source of superoxide but only when fully reduced. 28 In contrast, Lambert and Brand reported that site I Q is the predominant source of superoxide originating from Complex I. 21 The Ohnishi group found that both sites I Q and I F are ROS sources when they exist in the radical forms. 31 Thus, computational modeling shows that the differences among experimental data are likely due to different experimental conditions and that even seemingly contradictory results may be true under the right conditions.
Like complex I, the identity and origin of ROS produced by complex II remain unresolved. Until recently, it was thought to be a negligible source of mitochondrial ROS. Quinlan et al. found that site II F produces comparable amounts of ROS to site I Q in the presence of fatty acids. 20 The notion that site II F produces most ROS from complex II is shared by Siebels and Drose. 34 However, others contend that significant amounts of ROS originate from the site II Q 69 and the ISC near the Q site. 70 We used the model to identify conditions that favor ROS production by complex II and the contributions of its redox centers. In particular, the model identifies site II [3Fe-4S] as the primary source of superoxide arising from complex II. In FET, superoxide arising from site II FADH becomes more significant only under QH 2supported respiration (S/R). Thus, our model supports the experimental results that the ISC near the Q site and site II F give rise to most of complex II's ROS. Moreover, the model predicts that FET in non-phosphorylating mitochondria favors superoxide production by complex III. The bc1 complex also contributes significantly to total superoxide flux in phosphorylating mitochondria although its contribution is more pronounced in FET mode. Unfortunately, all these predictions remain difficult to experimentally test due to the lack of a robust, quantitative, and straight-forward method to directly measure them without interfering with other processes. An interesting theory relevant to this study is the concept of a partitioned or segmented Q pool. In the early days of the bioenergetics field, evidence against Q pool segmentation was provided with the aid of computer modeling. [71][72][73][74] These earlier studies are supported by more recent experimental evidence. 75 However, with the advent of new molecular and genetic tools, there has been a renewed interest in this theory. [76][77][78][79] These new data compel further research into this phenomenon. Based on the available evidence at this time, an integrated model is likely required to support evidence for or against this idea. It is quite possible that Q pool segmentation is a unique phenomenon associated with specific physiological circumstances that control cristae and inner membrane arrangement or that specific cells and tissues are wired with sufficient differences to observe this condition. Currently, the ETS-ROS model is unable to test this theory in its present form, but we believe the ETS-ROS model can be extended to test this theory.

Next Steps for Expanding the ETS-ROS Model to Include More ROS-Related Pathways
In developing this model, we focused on the known primary sources of mitochondrial ROS, the ETS complexes. Electron flavoprotein, 80-82 dihydrolipoamide dehydrogenase, 15,16,83 and other mitochondrial dehydrogenases 84,85 are additional sources of ROS that may become significant in different physiological conditions and tissues. In the next evolution of the ETS-ROS model, we will expand the model to include 2 major modifications. The first is to expand the model to explicitly include the full TCA cycle with additional enzymatic processes uncovered during the development of this model (eg, ME). This requires additional experiments to collect data on metabolite profiles under the various conditions. Second, we will develop models of beta-oxidation that are able to simulate both kinetic turnover and ROS production from electron flavoprotein and other associated ROS processes. This will, again, require additional experiments focused on fatty acid metabolism under various conditions spanning across the anticipated physiological domain.
Lastly, we will revisit the scavenging system to include glutathione and thioredoxin scavenging enzymes in addition to NAD(P) transhydrogenase. Despite the high level of consistency between model simulations and experimental data ( , this approach is incapable of partitioning the roles of the glutathione and thioredoxin systems involved in the ROS detoxifying pathway. Addressing these aspects require additional experimental data under the prevailing conditions to properly calibrate the model. Individual kinetic models of each of the four primary enzymes responsible for most of the H 2 O 2 scavenging are already developed, but they remain untested in an integrated model. [86][87][88] In all cases, each future generation of the model will be more refined and enable testing of hypotheses that are experimentally infeasible.

Conclusions
In conclusion, we have developed, analyzed, and corroborated a model of mitochondrial ROS that identifies the species-specific ROS originating from the ETS redox centers under varying respiratory states and electron flux directions. Each ROS-producing module is constructed in a thermodynamically faithful manner and tested prior to being unified in a single platform. [38][39][40] The biophysical details included in the model are supported by experimental data and further refined by the process of sensitivity analysis. Being modular, parsimonious, and consistent, this modeling approach enables the individual modules within the model to work harmoniously with each other in a thermodynamically faithful manner. This approach also facilitates the inclusion of other mitochondrial physiology processes that are related to ROS homeostasis, as mentioned above, when data become available. These qualities distinguish this model from existing models which also attempt to recapitulate mitochondrial ROS homeostasis 41,42 as ours is capable of consistently reproducing experimental data on various aspects of mitochondrial physiology and making predictions that are physiologically relevant.

Supplementary Material
Supplementary material is available at the APS Function online.

Author Contributions
Concept and design: JNB and QVD. Acquisition of data: JNB, QVD, MJD, JOSR, and YL. Analysis and interpretation of data: JNB, QVD, and YL. Development of computer code: JNB and QVD. Drafting or revising the article: JNB, QVD, and YL.

Data Availability
All data, model equations, and model codes described in this manuscript are given in the main body or supplement.