Reactive Oxygen Species Production by Forward and Reverse Electron Fluxes in the Mitochondrial Respiratory Chain

Reactive oxygen species (ROS) produced in the mitochondrial respiratory chain (RC) are primary signals that modulate cellular adaptation to environment, and are also destructive factors that damage cells under the conditions of hypoxia/reoxygenation relevant for various systemic diseases or transplantation. The important role of ROS in cell survival requires detailed investigation of mechanism and determinants of ROS production. To perform such an investigation we extended our rule-based model of complex III in order to account for electron transport in the whole RC coupled to proton translocation, transmembrane electrochemical potential generation, TCA cycle reactions, and substrate transport to mitochondria. It fits respiratory electron fluxes measured in rat brain mitochondria fueled by succinate or pyruvate and malate, and the dynamics of NAD+ reduction by reverse electron transport from succinate through complex I. The fitting of measured characteristics gave an insight into the mechanism of underlying processes governing the formation of free radicals that can transfer an unpaired electron to oxygen-producing superoxide and thus can initiate the generation of ROS. Our analysis revealed an association of ROS production with levels of specific radicals of individual electron transporters and their combinations in species of complexes I and III. It was found that the phenomenon of bistability, revealed previously as a property of complex III, remains valid for the whole RC. The conditions for switching to a state with a high content of free radicals in complex III were predicted based on theoretical analysis and were confirmed experimentally. These findings provide a new insight into the mechanisms of ROS production in RC.


Introduction
Reactive oxygen species (ROS) are side products of electron transport in the mitochondrial respiratory chain, the principal component of energy transformation in mitochondria. ROS generation starts with the formation of a superoxide radical (O 2 2 ) as a result of interaction between molecular oxygen and free radicals, e.g. semiquinone (Q 2 ): O 2 +Q 2 RO 2 2 +Q [1]. This extremely active compound can be deactivated in cells, mainly through superoxide dismutase [2]. However, H 2 O 2 formed in this process can interact with various intracellular compounds to produce ROS. ROS production serves as a metabolic signal [3][4][5]. However, when released in excess under certain stress conditions such as hypoxia/reoxygenation, ROS can also directly damage cells [6]. This destructive function of the electron transport chain represents the main problem in organ transplantation [7] and in many systemic diseases, as diverse as Parkinson disease [8] and diabetes [9]. The problem can be so great that in some organisms disruption of the electron transport chain can be a positive factor in increasing lifetime [10].
Although electron transport and coupled ROS production have been the focus of intensive research, important details are still not understood. There is currently debate regarding the relative contribution of various sites of the respiratory chain to overall ROS production [11,12] and the factors that may alter this contribution [13]. The use of specific inhibitors can localize the sites of ROS production, but their contribution under normal and stress conditions without inhibitors in vivo is not clear. It is generally accepted that electron transport from succinate through complex I to NAD + , the phenomenon known as reverse electron transport [14], is important for respiration and ROS production [15,16]. However, the mechanism of ROS production as a result of electron transfer upstream in the respiratory chain is not understood. Some details of the general mechanism of electron transport, such as the interaction of complex I with quinones that results in translocation of four protons through the membrane and reduction of one ubiquinone molecule per two electrons transported, remain the subject of discussions [17,18]. Answering these questions will help in understanding the mechanisms of electron transport and coupled ROS production, and will be useful for advances in transplantology and therapy.
The solution to such problems requires not only improvements in experimental techniques and new experiments, but also modification of methods for theoretical analysis. Specifically, kinetic modeling, which is an efficient method for investigating complex systems, still needs to be adopted for the mitochondrial respiratory chain. In fact, kinetic modeling in its classical form has been used for analysis of mitochondrial respiration. However, even the most detailed models [19] could consider only simplified scenarios. Huge number of differential equations is necessary to describe the behavior of respiratory complexes, so an automated procedure is required for their construction. Previously we developed a rule-based methodology for the automated construction of large systems of differential equations for analysis of 13 C isotope tracing experiments in metabolic flux analysis [20][21][22]. We extended this methodology to the mathematical description of multienzyme complexes, specifically mitochondrial respiratory complex III based on a Q-cycle mechanism [23]. A detailed description of complex III operation revealed that in a certain range of parameters complex III has the property of bistability, where two different steady states exist for the same parameters and the system can reach one or the other, depending on its initial state. Perturbations, such as fluctuations in succinate concentrations or temporal hypoxia, can switch the system from low to high ROS producing steady state. Such a switch explains the damaging increase in ROS production on reoxygenation after hypoxia.
The prediction of bistability for the mitochondrial respiratory chain was based on analysis of the Q-cycle mechanism for complex III. The contribution of other parts of the respiratory chain and linked processes that provide substrates must affect the properties of the respiratory chain. To study mitochondrial respiration as a whole, we extended the model of complex III [23]. The extended model includes the following elements: a detailed mathematical description of complex I; the stoichiometry of electron transport and proton translocation by the respiratory chain; the transmembrane potential; proton leak; oxidative phosphorylation; the TCA cycle that produces NADH and succinate as substrates for complexes I and III; and the transport of TCA cycle metabolites. The objective of this extension to the whole respiratory chain and linked processes was to create a tool for analysis of the basic behavior of the respiratory chain, in particular under conditions defining different fluxes in the forward and reverse directions. The ultimate aim was to reveal characteristics that have not been measured, such as the content of various free radicals, and thus to provide an insight into the relationship between states of the respiratory chain operation and ROS production. Figure 1 shows the components of the respiratory chain connected in the extended model. Respiratory complexes I to III Figure 1. Scheme for mitochondrial respiration and linked processes simulated in the model. Two reactions lead from pyruvate to succinate and further transformation to oxaloacetate reduce NAD + to NADH. The latter is used by complex I to generate a transmembrane electrochemical proton potential (Dm H+ ) and reduce ubiquinone (Q) to ubiquinol (QH 2 ), oxidation of which by complex III also contributes to Dm H+ . Complex III reduces cytochrome c, oxidation of which by complex IV and reduction of molecular oxygen to H 2 O is also coupled to Dm H+ generation. Oxidation of succinate to fumarate by complex II is coupled to the reduction of ubiquinone and thus fuels complex III. The product of electron transport, Dm H+ , is consumed for ATP synthesis. doi:10.1371/journal.pcbi.1001115.g001

Author Summary
Respiration at the level of mitochondria is considered as delivery of electrons and protons from NADH or succinate to oxygen through a set of transporters constituting the respiratory chain (RC). Mitochondrial respiration, dealing with transfer of unpaired electrons, may produce reactive oxygen species (ROS) such as O 2 2 and subsequently H 2 O 2 as side products. ROS are chemically very active and can cause oxidative damage to cellular components. The production of ROS, normally low, can increase under stress to the levels incompatible with cell survival; thus, understanding the ways of ROS production in the RC represents a vital task in research. We used mathematical modeling to analyze experiments with isolated brain mitochondria aimed to study relations between electron transport and ROS production. Elsewhere we reported that mitochondrial complex III can operate in two distinct steady states at the same microenvironmental conditions, producing either low or high levels of ROS. Here, this property of bistability was confirmed for the whole RC. The associations between measured ROS production and computed individual free radical levels in complexes I and III were established. The discovered phenomenon of bistability is important as a basis for new strategies in organ transplantation and therapy.
are components of the electron transport chain connected through ubiquinone. Complex III is linked to complex IV through reduction/oxidation of cytochrome c. NADH, which is a substrate for complex I, is produced in the TCA cycle. Since the total concentration of NAD + and NADH is conserved, NADH consumption, which fuels electron transport in the respiratory chain, defines the levels of NAD + , which is a substrate for several reactions in the TCA cycle. In this way, the extended model links electron transport with central energy metabolism, in particular with the reactions of the TCA cycle.

Determination of parameters by fitting experimental data
As described in the Methods section, the model of the respiratory chain and linked substrate transport and TCA cycle reactions contains 51 parameters. Out of 22 parameters of complex III, six ratios for forward and reverse rate constants were expressed through midpoint potentials. The order of magnitude for the rate constants for forward electron transport reactions in complex III can be estimated based on previous studies [19]. A qualitative reproduction of measured triphasic dynamics of cytochrome b H reduction by succinate in isolated cytochrome bc1 complex [24], as described in Text S1 and Figures S1, S2, S3 and S4, provides some restrictions for rate constants for binding/ dissociation of complex III with ubiquinone species.
The rates of respiration in the presence of ADP (state 3) or an uncoupler characterize the maximal capacity of the respiratory chain. In the absence of ADP (state 4), the respiration rate is characterized by proton leaks, which must be compensated by respiration. According to our measurements, the respiration rate is 480640 and 170630 ng atom O/min/mg protein in the uncoupled and in state 4 in succinate-fueled mitochondria, and 410630 and 80620 ng atom O/min/mg protein in mitochondria fueled by pyruvate and malate, respectively.
If mitochondria fueled by succinate do not expend the energy of the transmembrane electrochemical potential on ATP synthesis (state 4), succinate oxidation results in fast reduction of intramitochondrial NAD + . In the presence of rotenone, an inhibitor of electron transport in complex I, NAD + reduction is characterized by NAD + -dependent reactions of the TCA cycle and in particular the forward respiratory flux resulting from succinate oxidation. In the absence of rotenone, reverse electron transport [14] also participates in NAD + reduction, which makes the process much faster (Figure 2A). These data define the rate constants for reverse electron transport.
While succinate fuels complex III through succinate dehydrogenase, the oxidation of malate and pyruvate in the TCA cycle fuels complex I by reducing NAD + to NADH. Respiration under such conditions defines the characteristics of complex I. To evaluate the model parameters, we used a procedure that simulates all the different types of data listed above for the same set of parameters. The ratio of forward and reverse constants defined by a known midpoint potential or dissociation constant was kept fixed, and the conditions of substrate supply or membrane permeability for protons were changed in accordance with experimental conditions. The procedure fitted all the data by changing the free parameters within the order of magnitude indicated in [19], summarizing and minimizing the deviations in several calculations that simulated measurements. Minimization was performed using a standard stochastic procedure in the global space of parameters as described in Methods.
The best fit reproduces well the dynamics of NAD + reduction measured in brain mitochondria in the presence and absence of rotenone using the same set of parameters ( Figure 2A). The insets in Figure 2A show respiration rates and DY in the presence and absence of rotenone. These characteristics remain practically the same in both conditions. Without rotenone inhibition reversible electron flow through complex I, which fits the experimental data shown in Figure 2A, is directed to NAD + reduction (is negative) only during a short period of time ( Figure 2B), although ROS are constantly produced for a much longer time under such conditions [15,25]. Reverse electron flow is believed to induce excessive ROS production, but evidently these two processes are not correlated.
Rotenone essentially changes the dynamics of NADH measured before succinate addition. It is slightly oxidized by the RC in the absence of rotenone, but slowly reduced in its presence. This reduction is a result of oxidation of internal substrates while electron flow through the RC is blocked. We found that the metabolites of TCA cycle cannot be substrates that provide NADH reduction, because oxidation of TCA cycle metabolites results in much faster initial reduction of NADH. If the parameters of TCA reactions are changed to slow down and reproduce the initial dynamics of NADH, maximal respiration rate with pyruvate becomes inconsistent with experimental data (not shown). Rather, slow oxidation of other metabolites, probably aminoacids or lipids, contributes to NADH reduction. The simulation of such slow oxidation did not prevent NADH oxidation in the absence of rotenone, and reproduced NADH reduction in its presence.
While, in the absence of rotenone, succinate induced much faster NADH reduction due to reverse electron transport, the steady state levels are lower than in the presence of rotenone. The steady state levels are defined by NADH production and consumption in respiration. Rotenone blocks the consumption, therefore NADH levels further increase when rotenone is added after succinate. The model parameters were adjusted without considering subsequent NADH increase and the reproduction of this phenomenon validates the model.
The model reproduces measured maximal and state 4 respiratory electron flows for succinate-fueled mitochondria, as well as for mitochondria fueled by pyruvate/malate ( Figure 2C). The change in DY in the same simulations qualitatively corresponds to known changes measured under such conditions ( Figure 2D). The parameters for simulations shown in Figure 2 are listed in Table 1 (column indicated as best fit).
These simulations of measured data provide an insight into important hidden characteristics, such as the capacity of ROS production. ROS are produced by the respiratory chain as a consequence of one-electron transfer directly to oxygen from free radicals of electron transporters such as the semiquinone radical (SQ) at the Qo site in complex III [26][27][28] or FMNH [29], SQ bound to complex I [30], or N2 centers [30] in complex I. Simulating the experimental data as presented in Figure 2 the model at the same time simulates the dynamics of these free radicals.

Qualitative analysis of associations between the overall ROS production and individual radicals
The model describes various states of respiratory complexes formed in the process of electron transport, including those containing free radicals. Such radicals could be responsible for passing unpaired electrons to oxygen thus forming superoxide radicals and other forms of ROS. The contributions of various radicals to ROS production remain unknown; to clarify it we compared measured ROS production and the levels of various free radicals predicted by the model for the same conditions. A similar change in radical content and measured ROS production indicates qualitative accordance between the model and the described process and thus validates the model.
It is generally known that inhibition of reverse electron transport by rotenone decreases ROS production in succinate-fueled brain mitochondria [15,28]. In our measurements, ROS accumulation was inhibited immediately after rotenone addition ( Figure 3A). The model predicts that the SQ content at site Qo in complex III in succinate-fueled mitochondria is practically unchanged by the presence of rotenone ( Figure 3B) and this remains valid for simulations with any set of parameters describing the data well. This is the reason for the coincidence of intervals for SQ at Qo for the first two types of simulation shown in Table 1. Thus, the ROSstimulating role of reverse electron transport and the ROSinhibitory effect of rotenone cannot be explained at the level of complex III. Apparently, reverse electron flow mainly affects complex I by increasing the concentrations of free radicals able to pass electrons to oxygen.
The model predicts that rotenone essentially decreases initial levels of SQ bound on site Qn ( Figure 3C), FMNH ( Figure 3D), and the content of reduced N2 centers (inset). After an initial decrease, levels of SQ and FMNH increase, in agreement with the acceleration of ROS production measured after initial inhibition induced by rotenone ( Figure 3A). The reason for accumulation of free radicals and acceleration of ROS production is the production of malate from succinate, which then reduces NAD + in malate dehydrogenase reactions. This supply of substrate for complex I increases ROS production in rotenone inhibited mitochondria. The increase in the rate constant for malate-succinate exchange eliminates a slow increase in free radical content when rotenone is present (dashed curve in Figure 3C and D). It should be noted that the acceleration of ROS accumulation is not always observed experimentally and this agrees with the predicted disappearance of this slow component after acceleration of malate-succinate exchange. Such similarity of experimental and simulated behavior supports the mechanism accepted for its simulation and in this way validates the model.
The fact that the species of complex I can contain more than one radical makes it more difficult to understand the contribution of each site. In particular, the species 1101001 (the positions of digits correspond to Qp-Qp-Qn-Qn-N2-FMN-FMN), which contain SQ and FMNH radicals, slowly accumulate after inhibition by rotenone. This accumulation defines the dynamics . State 4 respiration (st4) was simulated using the parameters described in Methods. The action of rotenone was simulated by setting k fI5 = k rI5 = 0 (rate constants for N2-ubiquinone interactions). The maximal respiration rate (2u, uncoupled) was simulated by setting a high proton leak (k lk = 50000 s 21 ). Arrows denote the time of additions of mitochondria (Mt, usually at time 0), succinate (suc), and rotenone (rot) if they ate not present in the medium before mitochondria. doi:10.1371/journal.pcbi.1001115.g002 of SQ and FMNH, whereas 1101100 defines the fast component in levels of Qn-bound SQ and reduced N2 (inset in Figure 3D). It is possible that only one of coupled radicals makes the major contribution to ROS production, but in this case the levels of other radicals would also correlate with ROS production. On the other hand, radicals situated inside the same species could interact, so that the specie as a whole produce superoxide. In the considered example the behavior of the whole ensemble of radicals in complex I agrees with the observed effect of rotenone, and this validates the model.
Overall, according to the model predictions, rotenone hardly affects SQ levels in complex III, but initially it significantly decreases the levels of free radicals produced in complex I; this is the reason for the decrease in ROS production induced by rotenone in succinate-fueled mitochondria. The model also explains the subsequent increase in ROS production as a result of the formation of malate in rotenone-inhibited mitochondria.
Rotenone induces a large increase in ROS production in pyruvate/malate-fueled mitochondria ( Figure 4A). The corresponding simulations show that rotenone greatly increases the levels of FMNH and SQ at Qn site, but decreases the levels of reduced N2 ( Figure 4B). Since the changes in N2 disagree with measured ROS production, probably N2 center does not make essential contribution in ROS production under the considered conditions. The same species (1101001) that defined the slow component in the increase in free radicals now change faster and defines the main part of the response to rotenone. Species 1101100 also constitute an essential part of the total radical content, but their levels decrease in response to rotenone in accordance with the decrease of total N2 radical levels.
Stimulation of electron transport by addition of ADP or an uncoupler such as FCCP to succinate-fueled mitochondria results in a decrease in ROS production ( Figure 5A). This generally known phenomenon [15] validates the model prediction that the levels of all free radicals decrease when electron transport is stimulated by addition of ADP or an uncoupler.
Mitochondria fueled by pyruvate/malate also produce less ROS when electron transport is stimulated by an uncoupler ( Figure 6A). Such measurements also validate the model, which predicts a decrease in the levels of free radicals ( Figure 6B).
At high succinate concentrations, brain mitochondria produce much more ROS than those fueled by pyruvate ( Figure 7A). The model also predicts higher levels of free radicals in complex III, as well as in complex I, for mitochondria fueled by succinate ( Figure 7B).
Thus, the study of associations between measured ROS production and predicted radical levels in RC revealed qualitative consistency of measurements with all types of radicals and therefore validated the model, or showed a way of discrimination between possible sites of ROS production, and even between possible ROS producing species. However, in the latter case, a special, quantitative study is needed, which currently is beyond of the scope of presented study.

Prediction of bistability for the whole respiratory chain
It has been predicted that the Q-cycle mechanism of complex III can in principle induce bistable behavior [23]. The whole These intervals were calculated for each parameter separately among the sets of parameters that give x 2 below than a fixed threshold [31]. The sets of parameters were found using a stochastic optimization algorithm (Simulated annealing) that minimized the deviation from measured dynamics of NAD+ reduction in the presence and absence of rotenone, and uncoupled and state 4 respiration rates. Only forward rate constants are shown assuming that the reverse constants change proportionally keeping constant the ratio in accordance with respective DE m or dissociation constant. respiratory chain considered here, with the parameters that fit the experimental data, also has two different steady states for the same parameters. Figure 8A shows that the SQ content at site Qo in complex III could persist at different values, depending on whether the respiratory chain is initially reduced or oxidized. Figure 8B shows how the steady states for free radicals of complexes I and III change with the external succinate concentration for a set of parameters that reproduces the experimental data described above.
With increase in succinate concentration at some point the system switches to the state with the highest levels of semiquinone radicals at Qo site of complex III. The difference here from the similar curve in Figure 7 is that pyruvate is present in addition to succinate. Once the system is switched to the state of highest SQ content at Qo, it remains in this state even if the succinate concentration decreases back to low values. Thus, if the system is initially in an oxidized state, the steady state SQ levels at Qo depend on the succinate concentration, in accordance with the blue curve in Figure 8B. If the system is initially in a reduced state, it remains in this state until succinate concentrations decrease to the micromolar range. Since complex III is directly connected to complex I through a common substrate (ubiquinone), the bistable behavior of complex III induces bistability in complex I. However, when complex III enters the state with high SQ levels at Qo, SQ levels at Qn decrease ( Figure 8B), as well as the levels of other free radicals in complex I (not shown). In some range of succinate concentrations total amount of radicals in the two presented steady states can be similar, but this does not necessary means similar ROS production in the two states since the probability of ROS production can be different for various radicals.
Thus, bistable behavior remains valid for the extended model of the RC with proton translocation and transmembrane potential (DY) generation, and with parameters defined by fitting the experimental data and validated by qualitatively similar predicted and measured ROS production. The model predicts also that a pulse of succinate is associated with decrease of DY. Such counterintuitive decrease of DY induced by increase of substrate for respiration is shown in Figure 9A. The value of DY decrease, induced by the same pulse of succinate, can be different, depending, for instance, on membrane leak, as illustrated by two curves in Figure 9A. Measurements of DY using safranine O fluorescence revealed that the mean DY at low succinate (0.2 mM) is greater than at high succinate (2 mM) ( Figure 9B), thus validating the paradoxical prediction of the model.
The succinate threshold for a switch to the reduced state depends on the parameters of pyruvate transport and TCA reactions; here we do not investigate the quantitative details with respect to bistability, but emphasize only the qualitative similarity of predicted and measured behavior. With regards to the considered in the previous sections normal ''working'' steady state, the predicted levels of free radicals are robust with respect to the model parameters, as the next section shows.

Sensitivity to parameters and robustness of the model predictions
The sensitivity of simulations to variations in model parameters is shown in Table S1 for each type of experimental data presented in Figure 2 (dynamics of NAD + reduction, maximal and state 4 respiratory fluxes). The sensitivity is also listed for simulated levels of free radicals shown in Figure 3. The results indicate that significant changes in some parameters hardly affect the simulations (e.g. k qp_FS ). Evidently, the data do not restrict the parameter values and they could not be defined unambiguously. However, changes in these parameters within the range, for which fitting remains good, do not affect the predictions in terms of free radical levels. The parameters shown in red highly affect the simulations. However, it is possible that different combinations of such parameters could fit the measured data equally well because of mutual compensatory changes. In this case, despite the high sensitivity, the parameters can have a wide range of values for which a good fit is obtained. Confidence intervals rather than sensitivity are used to characterize the robustness of parameter determination. Different sets in the global space of parameters that fit the experimental data could be identified using our stochastic algorithm for minimization of the objective function x 2 (sum of squares of deviations from measured data normalized by standard deviations). The algorithm identified confidence intervals for parameters based on fixed thresholds of x 2 [31]. Table 1 shows the 99% confidence intervals for the free parameters. The ranges for which the values give a good fit to the data are large. Thus, even though the measurements cover various modes of respiratory chain operation, the data do not restrict the parameters sufficiently to define them unambiguously. Various sets over a wide range of parameters can describe the data equally well. However, the situation is different for free radical levels predicted for the simulated experimental conditions. Table 1 lists intervals for predicted free radical levels simulated using the parameters sets that fit the data with x 2 that is below the threshold. The confidence intervals for free radical levels are generally much narrower, so the predicted values are more robust. Although the intervals for SQ at Qo sites in succinate-fueled mitochondria are relatively large, they are clearly almost the same for both conditions (with or without rotenone). This result agrees with data indicating that the SQ content at Qo practically shows no dependence on the presence of rotenone ( Figure 3B). The levels of all free radicals in complex I under the conditions for the first two simulations are very robust, despite the high parameter variability. If the parameters give a good fit, the model predicts similar levels of complex I radicals. Although the intervals are relatively large under the third condition (pyruvate/malate supply), it is evident  that they are much lower than the intervals for the condition of succinate supply, as well as the levels of radicals in complex III.

Discussion
To construct a detailed mathematical model that accounts for all redox states formed during electron and proton transport in complexes III and I, we used our rule-based methodology for automated construction of large systems of ODE [23]. This model further extends our methodology previously used to model the distribution of 13 C isotopes in central metabolism [20][21][22], development of which occasionally coincided in time with that of similar rule-based methodology for signal transduction pathways [32,33]. For the study of mitochondrial processes our methodology gives a deep insight into the mechanics of respiration and ROS production. Here, rule-based algorithms for mathematical description of mitochondrial respiration coupled to proton translocation and DY formation was linked to a classical kinetic model that accounts reactions of the TCA cycle, which provides succinate and NADH as substrates for respiration and substrate transport in mitochondria.
After fixing the ratios of forward and reverse rate constants for electron transport reactions, free parameters were defined by fitting of forward and reverse electron flows measured under various conditions. High variability of parameters with a good fit to experimental data precluded definition of their values. However, the levels of free radicals calculated in the model showed much less variability. Different sets of parameters with a good fit to experimental data define very similar patterns for free radicals formed in complexes I and III. Thus, the analysis gives a valid insight into the mechanism of respiration and ROS production, even without precise evaluation of the model parameters.
A substantial body of experimental data on mitochondrial ROS production cannot be satisfactorily explained within the current experimentally based paradigm. Some of these results were obscure, such as acceleration of succinate-driven ROS production after initial inhibition by rotenone ( Figure 3). Others, such as a lower membrane potential in mitochondria fuelled by higher succinate concentration (Figure 9), were even counterintuitive. Calculation for mitochondrial constituents not measurable by current techniques represents a powerful tool for mechanistic explanation of accumulated data and for directing experimental research to test model predictions.
A body of evidence indicate that either FMNH [29,30], or SQ bound to Qn sites of complex I [30], or reduced N2 centers [30,34,35] may be a major contributor to ROS production, depending on the tissue, substrate, energy demand and oxygen tension [36,37]. The simulations revealed correlations between measured ROS production rates and levels calculated for each type of free radical. In this first step of the study we did not assume  any explicit link between any specific radical and ROS, but qualitatively compared all of them, taken separately, with measured ROS production. However, the method, which we use, opens a direction for future studies of quantitative contribution of various radicals of electron transporters, and even specific species of complex I and III, into total ROS production.
The similarity between changes in the ROS production rate and in the levels of specific free radicals validates the model and also provides an insight into the mechanism of ROS production. Rotenone inhibition of ROS production in succinate-fueled mitochondria correlated with the free radicals formed in complex I, but not in complex III. Evidently, under the given conditions, reverse electron transport must contribute to free radical formation in complex I, although the net flux reducing NAD + through complex I exists for only a very limited period of time.
In accordance with our previous study that revealed bistability for complex III [23], the extended model confirms the existence of two steady states for the same set of parameters. In one of these states (oxidized), mitochondria can develop a maximal rate of respiration, DY, and a capacity for ATP synthesis. This is the usual working state. In the presence of pyruvate high succinate concentrations can induce a switch of respiration to the reduced steady state, where lack of electron acceptors strongly restricts electron flow. The levels of free radicals in complex III greatly increase in this state, but decrease in complex I, in contrast. The switch to a more reduced state results in DY decrease. Indeed, we observed a DY decrease in isolated mitochondria in conjunction with an increase of succinate concentrations in the presence of pyruvate.
Q-cycle mechanism of complex III operation assumes bifurcation of electron flow at Qo site: one electron goes to Rieske center and further to complex IV, and another one reduces cytochrome b. This bifurcation of electron flow underlies the bifurcation between the two steady states. If in some moment the rate of first electron transition to Rieske center is higher than that for the second electron (because cytochrome b is reduced), semiquinones at Qo accumulate, thus preventing Qo liberation, binding and oxidation new molecules of ubiquinol, and thus limiting electron flow. In the case, shown in Figure 9A, greater proton leak resulted in greater transient discrepancy between the two electron flows at the point of bifurcation, which ultimately leaded to more significant inhibition of respiration and deeper descent of DY.
The decrease of DY, in the case shown in Figure 9B, is relatively small, however, in living cardiomyocytes a much higher decrease of DY can be observed, accompanied by high ROS production, and associated with mitochondrial permeability transition (MPT) [38]. Although the presented study does not touch a possible link between bistability and MPT, it puts forward some hypotheses, which verification in future can essentially clarify the mechanism of MPT. There are at least two phenomena, which do not find appropriate explanation in terms of current state of knowledge. First, the increase of permeability in this process does not increase electron flux and proton recirculation, as in case of uncouplers. Second, ROS production is high despite the decrease of DY. If we assume that the switch into the reduced state precedes MPT, both phenomena would find a natural explanation. This hypothesis, although not proved yet, opens avenues for deeper investigation of the MPT mechanisms. The presented new insight into mitochondrial respiration was possible due to the application of novel methodology of modeling that allowed a detailed mathematical description of mitochondrial respiration. The phenomenon of bistability, predicted based on this methodology, has a potential to be a basis of new paradigm for the mechanism of ROS production, which will initiate new research with outcome on academic and practical levels.

Methods
The file executable in Linux, which runs the simulations, and the C++ code of the program could be downloaded free from http://www.bq.ub.es/bioqint/ros_model/plcb2010.cpp.tar.gz.

Electron transport in complex III reflected in the model
The model of complex III described elsewhere [23] was used as a part of the extended model presented here. For each reaction two values, forward (K f ) and reverse (K r ) rate constants, were used as parameters. The order of magnitude of K f was set based on [19] and then K r was determined as described in [23] using midpoint electrochemical potentials, which determinations was variable and allowed refinement by fitting the data presented in ''Results''. Table 2 summarizes the reactions and values of parameters for complex III that simulate the data.

Electron transport in complex I reflected in the model
The overall process catalyzed by complex I is oxidation of NADH coupled with ubiquinone reduction and pumping 4 protons from negative to positive side of the membrane: NADHzH z z4H z n zQ'NAD z zQH 2 z4H z p ðI:0Þ This is a complex process that involves electron transport through a chain of intermediates coupled with proton translocations through inner mitochondrial membrane. The structure and mechanism of catalysis of complex I is reviewed in [39] and the data from this review are used for the construction of model of complex I. It is assumed that proton translocation is a result of Q reduction (with proton binding) at the negative side and its oxidation (and proton release) at the positive side. If several protons are translocated per one electron, then this electron must pass several cycles of Q reduction and oxidation. Such mechanism, similar to that accepted for complex III, called Q-cycle, was suggested for complex I (see e.g. [17]). We constructed a model based on electron cycling that is in accordance with the measured stoichiometry of proton translocations per one electron passed through the chain.
The initial step of such transport is the oxidation of NADH coupled with the reduction of FMN; further, electrons from FMN pass through a relay of eight different iron-sulfur (Fe-S) containing centers [40], which possibly form a relay for electron transport from FMN to the last Fe-S center N2 (see e.g. the review [40,41]). The Fe-S centers have similar midpoint potential close to that for FMN (E,2350 mV) with an exception of N2, which is much more positive (2150 mV, [40]). In this model the relay of Fe-S centers is simplified, so that electrons pass from FMN directly to the N2 center, which can interact with quinones. In this way, twoelectron transporter FMN and one-electron transporter N2 form the core of complex I, N2-FMN-FMN (referred as core).
The mechanism of interaction of N2 center with ubiquinone that results in the translocation of four protons from matrix to cytosol and one ubiquinol synthesized is not fully understood. Here we implemented in the model a proposed mechanism, which we consider as a working hypothesis that could be checked by the analysis of model behavior. In this way the model could serve as a tool for checking different possible mechanisms.
According to the EPR data [42,43] there are two ubiquinonebinding sites; bound ubiquinones possess different EPR characteristics, one of them is fast-and another is slow-relaxing. The former one bound in oxidized form in the proximity of N2, in Qn site, could be reduced by N2 and bind protons taking them from negative (matrix side of the membrane) (indicated as Qn below). The other one, bound in the reduced form in Qp site, situated in the proximity of Qn, can interact with Qn-bound semiquinone releasing protons to the positive (cytosolic) side of the membrane. This interaction of two quinones in fact is in agreement with the idea outlined in [39]  Reverse rate constants marked (K r ) are calculated from the respective forward rate constants and midpoint potentials (DE) as described in [23]. Although the reactions of electron transport are shown in simplified form as bimolecular, in fact they are performed (and simulated in the model) as transitions between the states of the whole complex (monomolecular). The units of rate constants are described in the legend of The proposed mechanism of N2-ubiquinone interactions, which we implemented in the model, is shown in Figure 10, and the individual reaction steps are described in the legend. It satisfies the known stoichiometry of proton translocation and ubiquinone reduction (four protons translocated and one ubiquinol synthesized per two electrons taken from NADH).

Reduction of oxidized FMN by NADH. In traditional form this equation is expressed as
NADHzFMNzH z 'FMNH 2 zNAD z ðI:0:0Þ Any of the forms of complex I with reduced FMN can receive two electrons from NADH, however, subsequent transitions require the interaction of three centers, N2, Qn and Qp. Therefore effective outcome produces only the reduction of FMN in the specie qnpc with ubiquinone bound at Qn and ubiquinol bound at Qp, which is reflected by binary number 1100000 corresponding to decimal 96. The reduction of FMN results in the production of redox state 1100011 (decimal 99): The forward and reverse reaction rates for this transformation are expressed in accordance with mass action law: v fI0~kfI0 |C 1100000 |NADH|H z ; v rI0~krI0 |C 1100011 |NAD z ðI:0:1Þ Here, as described for the complex III, ''0'' designates oxidized and ''1'' reduced states. The ratio of rate constants from I.0.1 could be found from the known redox potentials. Equilibrium constant for this reaction as a function of midpoint electrochemical potentials could be found from the condition of equality of electrochemical potentials at equilibrium: Figure 10. Interactions between N2 centers and quinones in complex I. Numbers above or below a species indicate the redox state of the complex as a combination of electron transporters. The last two digits indicate the presence (1) or absence (0) of two valence electrons of FMN (not shown graphically). The third digit from the right denotes the state of the N2 center, the next two digits from the right indicate the presence or absence of two valence electrons of Q/Q 2 /QH2 at the n-site. The next two digits from the right indicate the valence electrons of Q/Q 2 /QH2 at the psite. Numbers 0-8 above arrows denote individual reactions. 0, FMN reduction by NADH; 1, electron transition from FMN to the N2 center; 2, electron transition from reduced N2 to n-site ubiquinone. This interaction results in electron transfer from p-side ubiquinol to n-side semiquinone, which is coupled to binding of two protons taken from the matrix side and release of two protons to the intermembrane space. 3, ubiquinol thus produced is released and p-site semiquinone changes its position, releasing the p-site, which binds the released ubiquinol; 4, n-site semiquinone takes an electron from p-site ubiquinol and forms ubiquinol, taking two protons from the matrix, while the p-site semiquinone formed releases two protons to the pside of the membrane. 6, ubiquinol formed at the n-site dissociates and semiquinone bound at the p-site changes its location, binding to the n-site. 7, the non-paired electron of N2 is captured by n-site semiquinone, which subsequently takes two protons from the matrix and is converted to ubiquinol. 8, release of n-site bound ubiquinol, and binding of ubiquinol at the p-site and ubiquinone at the n-site. doi:10.1371/journal.pcbi.1001115.g010 ðI:0:2Þ taking into account that the difference between midpoint potentials for NADH (2320 mV) and FMN (2340 mV) [40] is DE m = 220 mV. 1. Reduction of the N2 center by FMN (step 1 in Figure 7): First electron of FMNH2, which by convention occupied second position from the right in binary representation, passes to N2 converting 0 into 1 in the third position from the right: The relationship between forward and reverse rate constants could be defined similar to (I.7). For the first transition at equilibrium In binary form:

1100101'1101001
The semiquinone Q 2 n is very active [17], so it reacts with QH 2 bound at p-site: In binary form: This reaction is symmetrical: p-side quinol and n-side semiquinone give p-side semiquinone and n-side quinol. The distance between the two quinone binding sites can be estimated as follows. Fastrelaxing semiquinone (bound to n-side oriented proton well) situated at the distance of ,12 Å from N2, slow-relaxing semiquinone (bound to p-side oriented proton well) situated at the distance of ,30 Å from N2 [42]. The distance between the bound quinones could be around 18 Å , which makes possible the interaction between them, taking into account the high energy of electron coming from FMN to Qn-bound quinone through N2 center.
The assumption of such interaction fulfills the known stoichiometry of translocation of four protons and overall reduction of one ubiquinone coupled with the transport of two electrons through complex I. We grouped together these two reactions: Overall in this reaction the oxidation of N2 center is coupled with the reduction of Q n .
RT nF |ln( FMNH|Fe 2z FMNH 2 |Fe 3z ) Since Here is considered that Em for N2 center is 2150 mV [40] and Em for ubiqunone one-electron reduction is 245 mV [42]. 3. Dissociation of QH 2n at n-site, transition of p-site SQ p to the n-site and binding of dissociated QH 2n at p-site.
In this step the three reactions are combined: dissociation of QH 2 formed at n-site, change of position of p-site bound semiquinone, and binding QH 2 at p-site. Overall in binary form:

1011001'1110001
ðI:3Þ and the forward and reverse rates are: v rI3~krI3 |C 1011001 v rI3~krI3 |C 1110001 4. Second electron (from radical FMN 2 , which by convention occupied the right position) passes to N2 converting 0 into 1 in the third position from the right: The transition of second electron characterized by the same DE m as accepted in (I.1), but it is not related with proton binding or release, therefore the right hand side value of (I.1.3) equals to the ratio k f /k r .
5. Reduction of N2 by FMN in step 4 induces the interaction of n-site semiquinone with p-site quinol resulted in the production of n-site quinol and p-site semiquinone coupled with the translocation of two protons: The dissociation of n-site quinol produced and the change of position of p-site semiquinone: 1011100'10100zQH 2 v fI6~kfI6 |C 1011100 v rI6~krI6 |C 10100 |QH 2 ðI:6Þ 6. The reduction of n-site semiquinone by N2 coupled with the binding of two protons: In equilibrium Since K eq~k f k r~( and the ratio of forward and reverse rate constants could be expressed as taking into account that Em for ubiquinol reduction in 263 mV [42] and Em for reduced N2 center oxidation is 2150 mV [40]. 7. QH2 dissociates, Q binds at n-site and QH 2 binds at p-site, overall: 11000zQ'1100000 v fI8~kfI8 |Q|C 11000 ; v rI8~krI8 |C 1100000 ðI:8Þ Above, the ratios of forward and reverse rate constants for the redox reactions of complex I are defined and summarized in Table 3. The particular values were defined from fitting the experimental data presented in ''Results'' using these ratios as restrictions. In some cases the fitting required different value of midpoint potential. This may indicate the differences in the operation of complex I in situ and under the specific conditions of midpoint potentials determination. Recognizing the importance of this subject we leave its studying for future because it deserves a separate consideration.
Although the mathematical description of complex I and complex III are similar, they differ in the strictness of rules for electron transport and proton translocation. For complex III the transition between two transporters allowed for any states of other transporters. This assumes participation of all 400 redox forms in electron transport. For complex I the rules accepted in the model allow participation in electron transport only several selected form. This illustrates the flexibility of methodology applied.
Transmembrane electrochemical gradient of H + Proton binding to ubiquinone at the matrix side of the membrane and their dissociation from ubiquinol to the intermembrane space results in the translocation of protons and arising the transmembrane gradient of H + concentration and electric potential. As described above for complex I, the reactions (I.2), (I.5) and (I.7) reduce ubiquinone each time taking two protons from the matrix. In complex III the reduction of ubiquinone at Qi site by reduced cytochrome b H is coupled with binding two protons taken from the matrix. The rate of this process (v 35 ) is calculated as described in [23]. The electron flow (v O ) through complex IV results in the reduction of oxygen with the uptake of two protons from the matrix and additional two protons are transferred from the matrix to cytosol. Proton leak (v lk ) and ATP synthesis (v syn ) return the protons transferred back to the matrix: v lk is leak of protons through the membrane: Here bo, bi, Vo, Vi are the buffer capacity and volume of outer and inner intracellular space with regards to mitochondria respectively. The differential equation for electric potential difference (y) used the same terms as that for proton concentration, but multiplied by a coefficient, which transforms the flux of ions into the change of electric potential: where F is Faraday number (96000 cu/mol or 0.96?1024 cu/ nmol), C is electric capacity of the membrane (2?1024 F/mg of protein, as computed based on [44].

Connection of respiration with central energy metabolism
Substrates for respiration, i.e. NADH and succinate are produced in TCA cycle inside mitochondria and in the model the connection of this part of intracellular metabolism with respiration through these common metabolites is taken into account by the simulations of following reactions.
Since the emphasis of work described here is the operation of respiratory chain, the reactions of TCA cycle were simulated in simplified form, as linear function of each substrate. Such expressions assume that the substrate concentrations are far from saturation, which should be true for the most cases. In this case the usual hyperbolic dependence of enzymatic reactions on substrate concentrations is close to the linear dependence. On the other hand, this simplification allows to avoid such unfavorable situation, when choosing inappropriate Km makes reactions artificially insensitive to substrate changes. Therefore we used such assumption as a first approximation, which could be easily corrected with obtaining more information about the properties of system.
Pyruvate transport and transformation to acetil coenzyme A: here the conversion of pyruvate into acetyl coenzyme A, linked with NAD + RNADH transformation, is included in the same reaction. The reactions converting citrate into succinate were joined together, taking into account that NAD + is used in these reactions: v tca~ktca |C NAD |C cit ; ðT:3Þ Then succinate is transformed into fumarate in succinate dehydrogenase reaction, which reduces Q taking two protons from matrix: v SDH~V SDH |Q|suc (K m(Q) zQ)|(K m(suc) zsuc) ðT:4Þ Here the total content of reduced and oxidized ubiquinone is conserved at the levels defined by [45] (6 nmol/mg prot).
Succinate not only could be produced in TCA cycle but also transported from outside of mitochondria in exchange to fumarate or malate (which are lumped in one pool in the present version of the model):  The differential equation, which describes the evolution of NADH takes into account the stoichiometry of its production in TCA cycle and consumption by complex I (reaction (I.0) described above.
The total concentration of NAD and NADH is assumed to be constant, so that NAD+, which defines the rates of TCA cycle reactions is computed as CNAD = CNADt2CNADH. The reactions linked with electron transport and respective parameters are summarized in Table 4. In total, without the reactions of pyruvate transport and ATP synthase, which were switched out in accordance with experiments analyzed, this module contains 11 parameters.
As the presented equations show, although the expressions for reaction rates are simplified, the stoichiometry of succinate and NADH production and succinate transport is reflected precisely in the model and this was the most important for the presented step of study of the link between central metabolism and ROS production by electron transport chain and the role of reverse electron transport in this process.

An algorithm for fitting experimental data
The whole model contains (2226)+(1827)+11 = 51 parameter (22 for complex III, 18 for complex 1, and 11 for the rest of reactions simulated). The six parameters of complex III and seven parameters of complex I are defined by the known values of midpoint potential. The other parameters were validated by fitting experimental data. To fit the experimental data our modification of Simulating Annealing algorithm was implemented in the way similar to that in [22]. The specificity of this algorithm was defined by the specificity of experimental data. The dynamics of NAD + reduction was measured under the two different conditions, in the presence and absence of rotenone, an inhibitor of reduction/ oxidation of quinones in complex I. The presence of rotenone was simulated by decreasing to zero the rate constants of step 5 in the reactions performed by complex 1 (k f15 = k r15 = 0), and the two conditions were fitted simultaneously for the same values of all other parameters. The procedure consisted of minimization of x 2 , normalized sum of squares of deviations from experimental data.
x 2 was calculated based on two simulations, first, normal conditions and, second, the presence of rotenone (k f15 = k r15 = 0, and all other parameters as in the first simulation).
The fitting algorithm made the following actions: 1. made the stochastic perturbation of given set of parameters (Vmax for the reactions of TCA cycle and substrate transport through the membrane) 2. performed coordinate descent, taking the parameters one by one and changing them in the direction, which decreased x 2 3. after reaching the local minimum of x 2 the program saved the respective set of parameters 4. returned back to step 1.
The cycles of perturbations and coordinate descent repeated thousands times and saved sets of parameters were analyzed: program read the saved sets with corresponding values of x 2 , defined the best fit (absolute minimum of x 2 ), the set of parameters, corresponding to the best fit, and defined confidence intervals for the parameters using the criterion of Dx 2 [31].

Experimental procedures
All procedures involving animals were approved by Children's Hospital of Pittsburgh and were in compliance with ''Principles of Laboratory Animal Care'' and the current laws of the United States.
Brain mitochondria were isolated from the cortex of adult Wistar rats. After removal, tissue was minced and homogenized in ice-cold isolation buffer I (IB I) which contained: 225 mM mannitol, 75 mM sucrose, 5 mM HEPES buffer (pH adjusted to 7.3 with KOH), 0.1 mg/ml fatty acid free BSA, 1 mM tetrapotassium EDTA and 12% Percoll. The homogenate thus obtained was carefully layered on the top of a discontinuous gradient of Percoll (24% and 42%) prepared using the same buffer. The preparation was then centrifuged at 16,0006g for 10 min. The fraction containing the mitochondria located between 42% and 24% Percoll was carefully withdrawn by a syringe and washed from Percoll twice by pelleting in IB I. The resulting mitochondrial suspension was diluted in isolation buffer II (IB II), which was similar to IB I, except for the concentration of EDTA (0.1 mM) and lack of albumin, and spun down at 12,0006g for 10 min. The deposit of mitochondria was homogenized in IB II at a final protein concentration of ,20 mg/ml and stored on ice until use. The protein concentration in the mitochondrial samples was determined using a Protein Assay kit (Pierce Chemical Company, Rockford IL) according to the manufacture's instructions. Mitochondria prepared in this way were active for at least 5-6 hours, as determined by their ability to maintain a stable transmembrane potential in the presence of oxidizable substrates. Fluorescence measurements were performed in a stirred cuvette mounted in a Shimatzu RF-5301 spectrofluorimeter maintained at 37uC. Mitochondria (0.2 mg/ml of protein) were added to 1.5 ml of the basic incubation medium that contained: 125 mM KCl; 2 mM KH 2 PO 4 ; 2 mM MgCl 2 ; 10 mM Tris;10 mM HEPES (pH 7.0); 100 mM EGTA; and oxidizable substrates as indicated in a particular experiment. Concentration of rotenone, when indicated, was 1 mkM.
Fluorescence of NAD(P)H was measured at excitation/ emission wavelengths 365 nm (slit 5 nm)/463 nm (slit 10 nm), respectively. To quantify the measurements a calibration curve was constructed using standard concentrations of commercial NADH.
Hydrogen peroxide was measured by fluorescence of Amplex red (2 mM), which increased upon oxidation to resorufin in the presence of 1 U/ml of horseradish peroxidase (HRP) as previously described [15]. Measurements were carried out at excitation/emission wavelengths of 560 nm (slit 1.5 nm)/590 nm (slit 3 nm), respectively. Amounts of H 2 O 2 released by mitochondria were estimated by constructing calibration curves using known H 2 O 2 concentrations in the standard incubation buffer together with Amplex red and HRP, but without mitochondria.
Mitochondrial transmembrane potential, DY m , was estimated using fluorescence quenching of the cationic dye safranine O. Since polarized mitochondria have a negative charge inside, positively charged molecules of safranine O are accumulated inside the matrix; increase in dye concentration inside the matrix leads to fluorescence quenching, thus, a decrease in fluorescence corresponds to an increase of membrane potential. The excitation wavelength was 495 nm (slit 3 nm) and emission 586 nm (slit 5 nm), and the dye concentration used was 2.5 mM [15].
Mitochondrial respiration rates were measured by an Oroboros High Resolution Respirometer (Innsbruck, Austria) in a stirred 2 mL chamber at 37uC in the same incubation media as indicated above. Oxygen sensor was calibrated at each experiment according to the manufacture's instructions. Calculations of respiratory rates were performed by software built into the instrument. Figure S1 The scheme of reactions performed by complex III as it is generally accepted (considered in Selivanov et al, 2009). One of two electrons taken from ubiquinol (QH2), which releases its two protons into the intermembrane space, recycles through cytochromes bh and bl reducing another quinone. The other electron continues its way to oxygen through cytochromes c1 and c and complex IV. Complexes I and II provide QH2. The reactions 0-12 are described in detail in the text. Table S1 Sensitivity of simulation of mitochondrial respiration with regards to the parameters. First column gives the list of parameters, next four columns give the relative change of respectively dynamics of NAD+ reduction in the absence and presence of rotenone, uncoupled respiration fueled by succinate, and pyruvate/malate. Next four columns give the relative change of SQ at Qo site of complex III, in the same four simulations as above, then, relative change of SQ at Qn site of complex I, then FMNH, and finally, reduced N2 centers. The highest changes marked black.