Modeling of Mitochondria Bioenergetics Using a Composable Chemiosmotic Energy Transduction Rate Law: Theory and Experimental Validation

Mitochondrial bioenergetic processes are central to the production of cellular energy, and a decrease in the expression or activity of enzyme complexes responsible for these processes can result in energetic deficit that correlates with many metabolic diseases and aging. Unfortunately, existing computational models of mitochondrial bioenergetics either lack relevant kinetic descriptions of the enzyme complexes, or incorporate mechanisms too specific to a particular mitochondrial system and are thus incapable of capturing the heterogeneity associated with these complexes across different systems and system states. Here we introduce a new composable rate equation, the chemiosmotic rate law, that expresses the flux of a prototypical energy transduction complex as a function of: the saturation kinetics of the electron donor and acceptor substrates; the redox transfer potential between the complex and the substrates; and the steady-state thermodynamic force-to-flux relationship of the overall electro-chemical reaction. Modeling of bioenergetics with this rate law has several advantages: (1) it minimizes the use of arbitrary free parameters while featuring biochemically relevant parameters that can be obtained through progress curves of common enzyme kinetics protocols; (2) it is modular and can adapt to various enzyme complex arrangements for both in vivo and in vitro systems via transformation of its rate and equilibrium constants; (3) it provides a clear association between the sensitivity of the parameters of the individual complexes and the sensitivity of the system's steady-state. To validate our approach, we conduct in vitro measurements of ETC complex I, III, and IV activities using rat heart homogenates, and construct an estimation procedure for the parameter values directly from these measurements. In addition, we show the theoretical connections of our approach to the existing models, and compare the predictive accuracy of the rate law with our experimentally fitted parameters to those of existing models. Finally, we present a complete perturbation study of these parameters to reveal how they can significantly and differentially influence global flux and operational thresholds, suggesting that this modeling approach could help enable the comparative analysis of mitochondria from different systems and pathological states. The procedures and results are available in Mathematica notebooks at http://www.igb.uci.edu/tools/sb/mitochondria-modeling.html.


Introduction
Throughout the mitochondria inner membrane are many energy-transducing protein complexes that help transform the chemical energy from the cell's metabolic intake into various useful forms of energy for the cell. Some of these complexes use the free-energy extracted froms reduction-oxidation (redox) reactions to transport proton across the membrane and establish a proton gradient, while others use this proton gradient in combination with the membrane potential, the proton-motive-force (pmf), to drive otherwise energetically unfavorable processes such as ATP synthesis and assorted transporters of other ions and/or molecules. This use of the pmf as an intermediate driving force in the overall conversion of energy is the essence of the chemiosmotic theory [1], and the flow of energy between these chemiosmotic complexes constitutes the core of mitochondrial bioenergetics [2]. The chemiosmotic complexes consist of complex I, III, IV, and V of the oxidative phosphorlyation (OXPHOS) pathway, and are partly encoded by the mitochondria DNA (mtDNA). Genetic variation or mutations in the mtDNA can alter the protein structures of the complexes, which can then affect their functional output, the bioenergetics of the system, and ultimately the health of the organism. In particular, a mtDNA mutation in a polypeptide of an electron-transport-chain (ETC) complex may cause its enzyme machinery to become less efficient in its energy transduction. An increase in the slippage of complex I and III of the ETC [3][4][5] can lead to an increase in the production of the respiration byproduct, reactive oxygen species (ROS), which can further damage the mtDNA and create a vicious feed-forward loop of energetic decline. When the total damage to the OXPHOS surpasses a functional threshold whereby it can no longer fulfill the energetic requirements of the cell, the cell may undergo apoptosis (programmed cell death) to remove itself from the population. The consequence of such an energetic decline is potentially grave, as over time, when enough cells are lost, the organism would begin to lose the functions of its organs, which might be manifested as either the normal progression of aging, or more seriously as the onset of major metabolic and degenerative diseases such as diabetes, Alzheimer, Parkinson, as well as cancer [6,7].
Interest in the roles that mitochondria play in mammalian health and disease has grown markedly over the past two decades, resulting in an abundance of genetic [7][8][9][10], structural [11,12], biochemical [13][14][15], and pathological [16,17] studies on mitochondria systems. An effective integration of the heterogeneous data coming from these studies is the main focus of the emerging field of systems biology [18]. However, the development of one of its key ingredients-the kinetic modeling of mitochondria bioenergetics-has not kept pace with the rest of the field and could potentially become a bottleneck.
Models of mitochondrial bioenergetics range from the top-level network-constraint variety [19], all the way down to the low-level molecular dynamics simulations [20]. Nevertheless, ODE based deterministic approaches still provide the best balance of dynamic descriptions and computational tractability across several length and time scales [21]. Of the deterministic models, the simplest and most direct approach is to approximate the respiration flux through the whole mitochondria by using a single empirical oxygen consumption equation [22,23]. To introduce a general thermo-kinetic approach, Jin and Bethke describe a respiration rate law that encompasses the overall electro-chemical reaction of the respiratory chain, and is applicable to both mitochondrial and bacterial respiration models [24,25]. Single equation approaches such as these, allow easy assimilation of cellular respiration into higher scale models, but lack the level of detail required to understand the contributions from the individual components.
At the next level of detail, the respiration process is divided into its principle components in the OXPHOS pathway, each with its own kinetic description. There are currently three main components-based OXPHOS models that serve as the basis for other larger and more extensive physiological models: the Yugi and Tomita model [26,27], the Korzeniewski model [28], and the Beard model [29][30][31][32]. These approaches have been used in the studies of in vivo cardiac energy metabolism [33] (Beard model) [34] (Korzeniewski model), dynamic OXPHOS respiration simulation [35] (Beard model), volume dynamics of mitochondrial bioenergetics [36](Beard model), mitochondrial fatty acid b -oxidation network [37](Yugi and Tomita model), the modeling of the ETC in purple non-sulfur bacteria [38] (Korzeniewski model), etc. Although the three OXPHOS models have been useful in studying several aspects of bioenergetics and physiology, they are limited by their choices of mechanism schemes. In particular, the Yugi model assembles a large array of detailed kinetic descriptions derived from specific enzyme binding mechanisms in past literature, e.g. the Ping-Pong bi bi mechanism of complex I [39] etc., but treats them as separate and independent ''reactors'' that do not incorporate the thermodynamic constraints necessary to characterize the dependence of their energy transduction processes to the chemiosmotic forces of the system. The Korzeniewski model primarily uses empirical data-driven relationships in both its reaction equations and systems properties, and it incorporates linear thermodynamic constraints on its OXPHOS components based on their free-energy profiles. However, the validity of its linear approximation is limited to near-equilibrium conditions. The Beard model inherits many components from the Korzeniewski model, but it extends the thermodynamic constraints to non-linear and far-from-equilibrium regions, and it explicitly treats the membrane potential and proton gradient of the system as separate state variables. However, the reaction rate equations in the Beard model lack detailed ''kinetic descriptions of enzyme activity'' [35]. In other words, these reaction rate equations do not intrinsically account for the kinetic properties of the ETC complexes since they were not derived in view of the internal mechanisms of the complexes. Instead, kinetic parameters are incorporated mainly through phenomenological control factors, which are introduced to compensate for specific modulations shown in experimental data sets.
The internal mechanisms of the complexes have been modeled with explicit elementary reaction steps to track intermediate metabolic species or reaction byproducts such as the generation of ROS [3,20]. However, at this scale, model validation becomes exceedingly difficult with large uncertainty as the individual reaction rates are not observable with the experimental technology currently available. For models that do incorporate detailed mechanistic schemes, there is the additional danger of overfitting in that mitochondria across different systems display a high degree of variability in their components [15], and consequently if a detailed model is based on a particular system, it may not generalize well to other systems. In addition, kinetic properties that are determined under a certain set of conditions (e.g. an in vitro laboratory setup) may not be transformable into another set of condtions (e.g. in vivo system). Thus, if a model of mitochondrial bioenergetics is to be both extensive and flexible, it must be able to adapt to the mechanism of a particular system, and be applicable to various system conditions. Furthermore, the dynamics of a system depend on the system's component descriptions, yet a system is often more than just the sum of its components. Thus, to establish a functional relationship between the dynamics of a system and its components, it is also necessary to allow for potential unobserved properties to emerge from the synergistic network interactions that vary from system to system. Only then can a model effectively capture the system's response to a perturbation across multiple scales.
In this paper, we present a new modeling approach for mitochondrial bioenergetics that addresses the problems of crossmechanism adaptability, cross-conditional applicability, and crossscale analysis through a novel composable kinetic rate equation that we term the chemiosmotic rate law. The rate law incorporates three configurable modulating factors that, via steady-state and rapid-equilibrium approximations, separately encapsulate the binding kinetics, the redox transfer potential, and the non-linear thermodynamic force-to-flux relationship of a prototypical energy transduction complex in the OXPHOS pathway. The separation of the factors allows one to selectively configure the mechanistic scheme of the complex, while the six biochemically relevant kinetic parameters of the rate law allow one to completely specify its kinetic properties at a particular reference condition. The kinetic parameters are estimated from data obtained experimentally using simple extensions to the standardized in vitro assay protocols for the various ETC complexes in the OXPHOS pathway. The resulting in vitro reference dynamics can be transformed into dynamics of an in vivo system through changes in the system's constituent thermodynamic forces, the reaction equilibrium constant, and the rate constants in a systematic way. One particularly useful type of transformation is a ''slippage'', or a loss in the thermodynamic efficiency of the complex's energy transduction process, which can be used to generally express a perturbation on the complex's reaction rate.
To validate the rate law, we conduct six parallel sets of assays on rat muscle homogenates, two for each complex. For each complex, we obtain kinetic parameter estimates from one of the two sets of assays, and then compare the progress curve predicted by our model to the experimental progress curve from the other set of assays. In addition, we formally establish the theoretical connections between our approach and other models, and compare their predictive accuracies. Finally, by using the optimized parameters for complex I, III, and IV, we perform a complete perturbation study of these parameters to reveal how they can significantly and differentially influence global flux and operational thresholds, suggesting that this modeling approach could help enable the comparative analysis of mitochondria from different systems and pathological states.

Methods
In this section, we first describe the derivation of the new chemiosmotic rate law based on the common properties of the ETC complexes. We then describe the experimental protocols that can be used to obtain the various biochemical parameters appearing in the rate law.

ETC Complex Rate Law Derivation
ETC Complexes and their Common Properties. The ETC consists of four integral membrane protein complexes that facilitate a sequence of electron transfer steps: (1) Electrons from citric-acid cycle products NADHzH z and FADH 2 enter the complex I (NADH-ubiquinone oxidoreductase) and complex II (succinateubiquinone oxidoreductase) respectively; (2) both complexes then catalyze the reduction-oxidation transfer of electrons to CoQ (ubiquinone) to produce CoQH 2 (ubiquinol); (3) subsequently, CoQH 2 is bound by complex III (ubiquinol-cytochrome c oxidoreductase), which transfers the electrons to cytochrome c; and (4) cytochrome c binds to complex IV (cytochrome c-O 2 oxidoreductase), which transfers the electrons to oxygen, the final electron acceptor [40] ( Figure 1A). The ETC complexes can either move about freely by lateral diffusion in the plane of the membrane ( Figure 1B), or alternatively, the complexes can form aggregates or supercomplexes, ranging from small clusters of a few complexes to a complete ETC assembly [11] (Figure 1C), to enable direct electron channeling between the complexes.
Each type of ETC complex is structurally unique, has diverse catalytic and binding rates, and responds to different inhibitors. Much of the detailed mechanisms of each complex remain to be determined. Nevertheless, all ETC complexes and their supercomplex assemblies share the same general characteristic of binding to both a donor and an acceptor electron carrier (EC) species, and facilitating the flow of electrons through the complex by a series of intermediate internal redox steps [2]. In addition, all complexes except for complex II capture the energy from such an electron flow and couple it to the drive of protons against their electrochemical gradient across the mitochondria inner membrane, much like how a windpump (combination of a windmill and a water pump) captures the energy generated by air flow to move water against gravity. Both are examples of an energy transduction process in which free-energy from an energetically favorable flow of one type of particle drives an energetically unfavorable flow of another type of particle through a type of coupling catalytic machine. However, as with all types of coupling catalytic machines, they can deviate from their normal coupling efficiency (slip) when damaged or when operating at an abnormal turnover rate.
Electron-Proton Pump (E epp ) Representation. To provide quantitative descriptions for the ETC complexes, an abstract prototypical ETC complex model, the electron-proton pump (E epp ) complex, is introduced based on the aforementioned general properties of the ETC complexes. Externally, this E epp complex is embedded in a membrane that separates two compartments, and its function is to catalyze the energy transduction between an electron transfer reaction and a proton translocation reaction (Figure 2A). This overall transduction reaction is expressed by the chemical equation: where D ox and D red are the oxidized and reduced species of the donor-EC reactant (D), A ox and A red are the oxidized and reduced species of the acceptor-EC reactant (A), H z in and H z out are the proton inside and outside of the membrane, while n dr , n do , n ao , n ar , m in and m out represent their stoichiometric reaction coefficients respectively. The overall chemical reaction for each of the ETC complexes is tabulated in Table 1.
Assuming that the concentrations of chemical species and reactants are homogeneous inside each compartment, the dynamics of the overall reaction is governed by a set of timedependent ordinary differential equations from the law of conservation of mass [41]. The rate of change for the chemical concentrations is directly related to the turnover rate of the reaction, expressed in terms of the reaction flux J, by their respective stoichiometric coefficients: In addition, depending on the composition of D and A, for each turn of the reaction there is n net number of electron transferred between D to A, giving the implicit electron flux: The turnover rate depends on various catalytic properties of the enzyme that the E epp complex represents, thus a detailed quantitative description of the reaction flux J for the E epp complex must take into account its internal structure and molecular mechanisms. Internally, the E epp complex can be modeled using four general features: one donor EC binding site, one acceptor EC binding site, one high potential redox-center R H , and one low potential redox-center R L ( Figure 2B). At any instant, the EC binding sites can either be occupied or unoccupied, while the bound ECs and the internal redox-centers can either be reduced or oxidized, giving rise to a total of 64 microstates (2 6 distinct configurations) of the reaction system. Of these microstates, 36 are distinct configurations of the complex that contribute to the forward or reverse turnover rate ( Figure 3). Transitions between these discrete microstates arise from the infinitesimal thermal fluctuations in the system, and they can be partitioned into the following elementry reactions that make up the overall reaction: two EC substrate binding reactions that consist of the EC binding and unbinding transition events (D and A binding, Figure 2B); an electron exchange reaction between the donor EC and R H (represented by the reaction rate v D{RH ); an internal electron transfer reaction between R H and R L that can be coupled to the translocation of protons across the membrane (v TT ); and an electron exchange reaction between R L and the acceptor EC (v RL{A ) ( Figure 2D). An averaging of these microstates and elementry reactions over an ensemble of E epp complexes provides the necessary link to the more practical macroscopic states of the system governed by the principles of thermodynamics.
Micro-Macroscopic Thermodynamics. The laws of thermodynamics are phenomological in nature, but they provide a convenient and powerful method of relating experimental bulk properties of a system such as pressure, volume, temperature and composition, which allows one to obtain desirable information on a system even if explicit knowledge of the interactions within the system is not available. Statistical mechanics provides the link between the quantum mechanical molecular properties and the macroscopic properties of thermodynamics by predicting an appropriate thermodynamic function of a system from its molecular structure and intermolecular forces. In regards to the reaction flux J, a crucial prediction derivable from statistical mechanics is how the thermodynamic forces of a system (which can also be expressed in terms of the Gibbs free energy gradient) can affect the magnitude and direction of J at steady-state.
Macroscopically, in a closed-system where there is no exchange of material with the outside, all chemical and physical processes eventually reach a balance such that there is no net activity in the system. At this thermodynamic equilibrium, the ratio between the concentrations (or the difference between the chemical potentials) of the products and substrates of a chemical reaction reaches a constant value defined as the equilibrium constant for that reaction. For any other ratio, there exists a disequilibrium in the chemical potentials, which provides the thermodynamic force X to drive the reaction towards equilibrium. By the same token, a disequilbrium from an electrostatic potential difference across a membrane also contributes to X . Living systems are open-systems in which materials are constantly being transported into and out of the system, driving the potentials of the various reactants and species far from their equilibrium. Analysis of non-equilibrium conditions is more intricate, but one can extend much of the equilibrium analysis to steady-states in which the system reaches a stationary condition while the input and output of the system are maintained at the same constant rate. At the steady-state, all potentials settle at a ''local equilibrium'' in which X remains constant, and continually drives the reaction with a constant net flux J.
In the following subsections, we use a combination of statistical steady-state thermodynamics, macroscopic equilibrium thermodynamics, and kinetics theories to derive a quantitative expression of the reaction flux J of the E epp complex as a function of its thermodynamic force X and the substrate concentrations in the reaction system. The resultant rate law provides a unification of the key attributes from the existing approaches we described [24,26,28,29]. In particular, we retain in the framework of our rate law, the convenient modular approach of Jin and Bethke that partitions the different contributions to the reaction flux into separate factors that modulate on a maximum reaction velocity. The derivation proceeds first by focusing on the steady-state thermodynamics of the overall redox reaction in a post-binding ternary DE epp A complex (thermodynamic force function), then coupling the net steady-state flux to the average transition rate between the internal redox centers (redox state function), and finally superimposing the kinetics of the binding reactions that produces the active ternary DE epp A complex (saturation function). In addition, we extend the derivation to allow adaptation of the flux expression to different system settings by means of transformations of rate constants and/or equilibrium constants from a reference system. We then utilize the transformation framework to introduce the chemiosmotic free energy transduction mechanism of the E epp complex, as well as the slippage mechanism in the transduction process.
Thermodynamic Force Function T(X). Let us first consider the redox mechanism of the E epp complex after it is bound with both its donor and acceptor ECs to form the active ternary DE epp A complex by using the cyclic kinetic diagram in Figure 4. The cyclic diagram is a transition map for the 16 fully bound microstates of Figure 3 condensed into three stable operational states (E zz , E {z , and E z{ ), where contributions from D and A are implicit in the transfer of electron(s) to and from the E epp complex (transitions 1<2 and 3<1), while the complex shuttles the internal transfer of electron(s) between the source and the sink (transitions 2<3). J z is the forward cycle flux that traverses the three states in a forward reaction sense (counterclockwise), while J { is the reverse cycle flux that travels in the opposite sense (clockwise). In a large ensemble of E epp complexes, some complexes would be operating in the forward sense (J z ), while others would be operating in the reverse sense (J { ), and the observed net reaction flux J is the difference between these cycle At equilibrium J z and J { are in balance and the net flux produced is zero, whereas at a steady-state away from equilibrium, a net flux is produced by the total thermodynamic force X of the system. In Figure 4, this X is the redox force where F is the Faraday constant, n the number of electron(s) transferred, and is the redox potential difference for one electron between the electron carriers D and A. In Equation 6, DE m is the standard mid-point redox potential [2], R is the universal gas constant, T is the temperature, DY is the electrical potential difference across the membrane, while e D?A is a switch with the value of 1 if the electron is transferred from the negative N side to the positive P side of the membrane, 21 if transferred in the opposite direction, and 0 if the electron does not cross the membrane. For The three internal processes of the electron-proton pump: reaction v D{RH captures the transfer of electron(s) from reduced form of the donor electron carrier D red to the high potential internal redox center R H ; reaction v TT tunnels the electron(s) between R H and the low potential redox center R L and drives the proton translocation; reaction v RL{A allows the exit of the electron(s) from R L to the oxidized form of the acceptor A ox . D.
Chemical reaction equation representation of the three internal processes: the red bidirectional arrows indicate the rate limiting step, and the blue bidirectional arrows represents processes under fast equilibrium condition relative to the rate limiting step. doi:10.1371/journal.pone.0014820.g002 mitochondrial inner membrane, N is the matrix while P is the inter-membrane space. By applying detailed balance on the individual elementry transition rates between the states in a cyclic system, a flux-to-force relationship is found that equates the steady-state ratio of the forward and reverse cyclic flux to the exponential of the thermodynamic force [41][42][43]: Combining Equation 7 with the definition of J in Equation 4, one obtains the two equivalent equations and in which the net reaction flux can be expressed exclusively as a function of the forward flux J z and X , or as a function of the reverse flux J { and X . Although equivalent, we note that in the direction and limit of X ??, Equation 8 provides a more convenient and numerically stable expression of J, while Equation 9 is better used in the direction and limit of X ?{?. For regions in the domain of X close to the equilibrium, either equation is just as good and the choice depends on whether J z or J { is more obtainable. Here the thermodynamic force function is introduced as: which encapsulates the modulation on J by X , and simplifies the expression of J to: Redox State Function R(D r=o ,A o=r ). After establishing the flux-to-force relationship of the general redox reaction of the E epp complex, the next step is to find the expressions of cycle fluxes J z and J { . In the internal redox reactions defined earlier, the transfer-transport reaction v TT is a lump reaction of all intermediate electron transfer steps between the two outer-most redox reaction sites and also where the proton transport may occur. Therefore, comparing to v D{RH and v RL{A , v TT can be assumed to be the rate-limiting portion of the overall reaction, which, at the steady-state, can be coupled to the steady-state expression for J.
Microscopically, v TT represents the transitions between the two states E {z and E z{ (Figure 4 and 5). Analogous to the net reaction flux (Equation 4), the net transition flux between two microstates is the difference between its forward and reverse transition rates: where i and j are indexes referring to the two distinct neighboring states, p i and p j are the corresponding state probabilities, and a ij and a ji are the transition rate constants associated with the    transitions from i to j and from j to i respectively. For v TT , i~2 and j~3 and the state probability p 2 is the combined probability that R H is in its reduced state R red H while R L is in its oxidized state R ox L , a necessary condition for the forward reaction (see Figure 5). Similarly, the state probability p 3 is the combined probability of R ox H and R red L , the condition necessary for the reverse reaction. Taking into account that each E epp complex consists of one R H and one R L , the probability of a redox-center to be in either oxidized or reduced state over an ensemble of ternary DE epp A complexes is determined by dividing the concentrations in each state by the concentration ½DE epp A: which are subject to the constraints: The probability for the forward state can then be taken as and the reverse state as Substituting for the probabilities in the net transition flux (Equation 12) and redefine a 23 and a 32 as kinetic rate coefficients k z TT ½DE epp A and k { TT ½DE epp A, one arrives at the kinetic expression for v TT : The equilibrium constant for the reaction is found by noting that v TT~0 at the equilibrium: where the probabilities P(:) e are the state probabilities at the equilibrium. Equating the rate-limiting transition flux in Equation  17 to the steady-state expression of J in Equation 4 gives an updated expression of Equation 11: is the thermodynamic force in the range of the v TT reaction. Equation 19 now completely describes the steady-state flux of the E epp complex, but the actual state values for R H and R L could not be measured directly since they are states within the complex. Nevertheless, a correlation between the concentrations of the external ECs and the internal redox-centers could be made such that the internal state values can be inferred from measurements of external concentrations. Since v TT is the rate-limiting reaction, it follows that the kinetic constants for the exchange reactions v D{R H and v R L {A are much larger than k z TT and k { TT , and the exchange reactions may be assumed to achieve rapid-equilibrium. As described by Jin and Bethke, with this assumption the thermodynamic forces for the two boundary exchange reactions v D{R H and v R L {A can be approximated as zero: This allows one to approximate the total thermodynamic force spanning from D to A, as: and expand the reaction boundary of X TT in Equation 20 to include the chemical potentials of the external ECs. Furthermore, one can obtain from Equation 21, the equilibrium constants for which can be rewritten as the correlation between the internal concentration to the external concentrations: Using the algebraic constraints in Equation 13 and Equation 14, the relationships between the probabilities of each internal state and the concentrations of the external reactants are expressed as: Here shorthand notations for the ratios of the redox states are introduced: Substituting the relationships in Equation 26 and 27 for all the internal probability values in Equation 19 and Equation 20 allows the reaction flux to be expressed in terms of the external reactant species: where Noting the exponential form of K R,D and K R,A in Equation 24 and combining them with the midpoint potential DE mRH ,RL gives which is nearly the same as the definition of X redox in Equation 5 and 6 (identical if e D?A~0 or if the electron is assumed to traverse the membrane between R H and R L only).
Here the redox state function is introduced as: which encapsulates the modulation on J by the redox state ratios, and simplifies the expression of J to: Saturation Function S(D T ,A T ). Equation 32 describes the flux as a forward rate coefficient k z TT that is scaled by the concentration of the ternary complex ½DE epp A, and modulated by both a thermodynamic force function T and a redox state function R, both of which change with respect to changes in the donor and acceptor redox state ratios. However, it is important to note that the two functions describe only the redox electron transfer processes of the E epp complex, but not the binding kinetics of the complex with respect to the donor and acceptor ECs (D and A). Indeed, binding kinetics is a central focus of several kinetics based models of ETC complexes, such as the Yugi and Tomita model [26,27] or the more recent Chen and Beard model of complex I [44], because most enzyme regulation studies deal with the binding of substrates. However, the binding mechanism for each ETC complex may be unique and may also vary significantly across different mitochondrial systems, as suggested by the lack of consensus in the literature. Thus the strategy adopted here is to maintain the generality of the rate law so that it remains compatible with different commonly found mechanisms.
At steady-state, the electron exchange reactions between the ECs and the E epp complexes are assumed to have achieved a rapid-equilibrium. Furthermore, since the binding of the EC reactants to form the ternary DE epp A complexes precedes the exchange reactions, one can assume that the binding reactions must also be a part of the transient kinetics that takes place before the establishment of the steady-state. This enables us to superimpose the binding kinetics of EC reactants on Equation 32 by finding the fraction of the total E epp complex ensemble that forms DE epp A. An extended kinetic diagram in Figure 6 illustrates how the binding of the ECs separates the E epp microstates from wasteful cicles and transitions to the active cycles that contribute to the reaction flux. The more the E epp ensemble is ''saturated'' with (28) the reactants, the higher the probability that the E epp ensemble manifests itself in the DE epp A form. The extent of this saturation can be quantified by the ratio ½DE epp A=½E epp T , where Using the nomenclature of Cleland [45], the binding mechanism for a typical two-substrate ETC redox reaction can be classified by the reactant binding sequence (ordered sequential, random sequential, and ping-pong), and molecularity (uni, bi, ter, quad, etc) [46]. In an ordered sequential mechanism, substrates binding to and products release from the enzyme follow an exact order. In a random sequential mechanism, the order of binding between the two substrates or the order of release between the two products are random. In a ping-pong mechanism, one or more products must be released before all substrates can react.
Since the ping-pong mechanism requires two separate catalytic steps, whereas a major assumption for the E epp complex is the single rate-limiting catalytic step v TT , it is not incorporated at this time. The general scheme for reversable ordered sequential bi bi mechanism and reversable random sequential bi bi mechanism are EzD red  Since the relative difference in the oxidized and reduced species of either the donor or the acceptor EC are already accounted for in the redox state function R, one can combine the contribution of both redox species of an EC reactant into one state variable through the constraints: and for irreversible random sequential bi bi as found in [45,46]. In the situation where only one of the two substrate EC concentrations is varying while the other is held constant, one can essentially consider the effects of the individual variation separately as two separate and parallel reactions, each with a single substrateenzyme binding step At steady-state, such independent binding reactions can always be expressed using the Michaelis-Menten-type kinetics [41]. The effects of the two independent reactions can be combined multiplicatively to give where K S,D and K S,A are the Michaelis-Menten-like saturation parameters that are characteristic of the enzyme complex. Equation 43 can be rearranged to give a comparable form to the sequential mechanisms: Here a saturation function is introduced as:  to set X~0, the equilibrium constant of the reaction in the reference system is found to be which relates the ratio of the two maximum velocities to the midpoint redox potential of the reaction. System Transformation. One important point in using Equation 46 is that although Equations 8 and 9 are equivalent over the entire domain of X , they are only equivalent if the equality in Equation 7 is not ''perturbed''. If the composition of X changes such that an additional force X d acts on the system, then the equilibrium point of the flux ratio ought to change with respect to the new force, and correspondingly the left hand side of Equation 7 ought to change in the same amount, giving the new relationship: where of the E epp complex from a reference system to a new system with the added force. In general, for each additional force, an additional pair of exponential terms is applied to Equation 49. Thus, for a number i of additional forces, the J f of the final system is expressed in terms of the J of the original system by: Free Energy Transduction and Enzyme Slippage (c). The standard form of the chemiosmotic rate law in Equation 46 considers only the redox thermodynamic force X red of the reference system. Such a system exists when proton gradient is not available either because the system cannot compartmentalize protons (i.e. a continuous membrane is not present to keep proton concentrations apart), or because the proton concentrations exactly balance across the membrane through clamping, both of which are observable and controllable in in vitro experiments. Thus, the reference system represents a fundamental basis from which the rate law for many other in vitro or in vivo systems can be derived through the transformation framework of Equation 50. The incorporation of proton gradient can be viewed as such a transformation. The free energy transduction reaction in a chemiosmotic complex is driven by an overall thermodynamic force X chemio that is the sum of the two opposing forces: where X redox is the redox force defined in the reference system, X pmf is the proton motive force (pmf), F is the Faraday constant, and n and m are the stoichiometric values for the number of electron transferred and net proton translocated respectively. In the pmf, is the energy necessary to pump one proton across the membrane with respect to the proton gradient and the membrane potential. Schematically, the addition of X redox transforms the kinetic diagram of the redox reaction in Figure 4 to the diagram of the transduction reaction in Figure 7. Energy transduction occurs when the free energy of the electron carriers decreases by an amount X redox , and from this, an amount X pmf is used to increase the free-energy of the protons. At equilibrium, X redox is in balance with X pmf , but an imbalance between the two would produces a net thermodynamic drive. In accordance with Equation 49, X pmf is a negative perturbation force {X d upon the reference system, and the additional parameter f p is used to determine the fraction of the effect of the the perturbation force that is distributed on the reference forward Jz and reverse J{ fluxes. The transformed expression of J is then: Since the equilibrium constant can be expressed as the ratio of the rate constants, the equilibrium constant of the transformed reaction with respect to the reference reaction is then In general, a transformation such as the addition of the pmf would shift the equilibrium constant of the original reaction (Equation 54); however, if the sum of perturbations in Equation 50 affects the forward and reverse rate constants in an equal but opposite manner, the original equilibrium will be preserved. Another possibility is if the perturbation force is a function of the original thermodynamic force, then the equilibrium concentrations of the reactants will not changed (albeit the equilibrium constant would be modified). One such perturbation is the ''slippage'' in the free energy transduction process. Energy transduction processes are prone to slippages in which efficiency can be affected by several factors [47] such as increased proton leakage or the loss of electrons to form ROS [6]. As a simple illustration, the efficiency of the flux-force relationship in Equation  7 can be compromised if a short circuit occurs in the cyclic states of the free energy transduction process (Figure 8). Alternate enzyme transition cycles could diverge from the normal transduction path, dissipating portions of the free-energy acquired from the high energy substrate without performing the transduction on the secondary substrate. This decrease in the available free-energy, expressed in terms of a smaller thermodynamic force, would lessen the magnitude of the net transduction flux according to Equation 11. To describe this lost of thermodynamic force X slippage without explicitly expressing its content, a convenient c variable is introduced: such that c represents the percent of original thermodynamic force X original available after losing a fraction through the slippage transition path. c has the range between 0 to 1 as X slippage has a upper bound of X original . Setting X slippage as the perturbation force in Equation 49 but expressing X slippage in terms of c and X original with Equation 55 gives the slippaged corrected expression of J as

Methods for Determining the Kinetic Parameter Values
The standard form of the chemiosmotic rate law, assuming a parallel binding mechanism (Equation 43), contains a total of six kinetic parameters: apparent maximum forward and reverse velocities (V f ' max ,V r' max ), donor and acceptor reactant saturation constants (K S,D ,K S,A ), and donor and acceptor redox state constants (K R,D ,K R,A ) ( Figure 9). Through the equilibrium constant in Equation 47, the six parameters can be reduced to five as either V f ' max or V r' max can be expressed in terms of the other. One of the main advantages of this rate law is that all five of its basis kinetic parameters can be fully determined through enzyme kinetic studies (e.g. the consensus protocols of respiratory chain spectrophotometric assays for clinical diagnosis http://lbbma. univ-angers.fr/lbbma.php?id = 58). The principle technique for these assays involves the use of an ultraviolet/visible (UV/VIS) absorption spectrophotometer in which the time-course conversion of a redox substrate species to its product species by an ETC complex in a closed system (inside a curvette) is recorded to obtain the initial velocity (rate) of the reaction. Homogenated tissues and isolated mitochondria (where membranes are fragmented in both cases) are specifically used as the reference system since its pmf can be neglected due to the absence of an intact mitochondrial inner membrane. Thermodynamic constants from the literature and saturation reactant concentrations used in our experiments are listed in Table 2. For our studies, the protocol is extended to provide a complete time-series of the reaction until the substrate species is completely exhausted or the system reaches an equilibrium.
In the following sections, procedures to determine the five kinetic parameters from the time-series data are described. With the exception of the kinetics for the final electron acceptor in complex IV, where both oxygen and water are reactants that are open to the bulk concentration, these procedures apply to all kinetic parameters in the ETC complexes. All of the fitting and subsequent simulations of the rate equation are performed using Mathematica 89s NonlinearModelFit function, which produces leastsquares fits that are defined to minimize the quantity x 2~P i r i j j 2 , where the r i are residuals giving the difference between each original data point and its fitted value. The procedures and results are available in the form of Mathematica notebook files at http://www.igb.uci.edu/tools/sb/mitochondria-modeling.html.
Determining the Maximum Forward Velocity (V f ' max ) and the Saturation Parameters (K S,D andK S,A ). As indicated in the derivation of Equation 45, out of the three modulating function, the S saturation function affects the reaction flux first. Thus its two saturation parameters, K S,D and K S,A , can be obtained from the instantaneous initial velocity of the reaction when there are no products so that functions R and T have negligible contributions. Furthermore, since it is assumed that the two binding reactions in S are independent to each other, K S,D and K S,A can each be determined separately by varying the starting concentration of the corresponding EC substrate while saturating the EC substrate of the other parameter to minimize its contribution. Consequently, for each complex parameter determination, there are necessary two sets of time-course enzyme kinetics assays: (1) the D red set with ð53Þ ð57Þ variations in the initial concentration of the reduced donor substrate specie; and (2) the A ox set with variations in the initial concentration of oxidized acceptor substrate specie. To illustrate, Figure 10A shows the time-series of complex I with variations in the initial concentration of NADH, while Figure 10C shows the time-series of complex I with variations in the initial concentration of CoQH 2 . For each assay set, an initial velocity is approximated for each of the initial concentrations by calculating the change in concentration over a starting time period (dependent on the amount of enzymes used). The initial velocity-concentration pairs are then used as data points for calculating the residuals in the least-square fitting of the corresponding Michaelis-Menten-like factor in Equation 43 to obtain the associated saturation parameter ( Figure 10B and 10D). A final value of V f ' max is obtained by averaging the fitted initial velocity value from both D red and A ox variation sets.
Determining the Redox State Parameters (K R,D and K R,A ). Beyond the transient period of the instantaneous initial velocity, both R (Equation 31) and T (Equation 10) contribute significantly in the modulation of the reaction flux. The T function can be fully described by the time-dependent variables D r=o and A o=r and the thermodynamic constant DE mD,A (Equation 30 with pmf &0), leaving only the R function to be determined. R imposes on the flux, a hyperbolic dependence on the ratios D r=o and A o=r through the redox state parameters K R,D and K R,A , respectively. Both ratios change continuously as the reaction progresses, and their rate of change are related through the stoichiometric coefficients of their constituent species in the overall reaction (Equation 1 and 2); therefore, the effects of K R,D and K R,A are not separable in the reaction, and must be considered together.
Given all other parameter values are fixed or determined, the values for K R,D and K R,A can be determined by fitting the simulation output from Equation 46 to the time-series from the D red or A ox variation assay sets. Although each time-series from the two sets of assays differs in its total reactant concentrations, the values of D r=o and A o=r only depend on the concentration ratio of the respective redox-pairs at the specified time. Thus, ideally, every one of the time-series can supply the full value range of the ratios to obtain an estimate for K R,D and K R,A . However, because K R,D and K R,A must be determined simultaneously, and the factors containing them in the R function are symmetrical, the estimated parameter values from fitting a single time-series are not unique and might not provide the optimal representation of the complex reaction over various conditions. This is demonstrated in Figure 11A and 11B, where the best fit parameter values estimated individually from the time-series 50mM, 75mM, and 100mM of the NADH assays are shown to have a large variation in the parameter values, and the subsequent simulations of the rate law based on one of the three estimates are shown to have a poor agreement with the time-series from the 25mM, 50mM, 75mM, and 100mM CoQ assays (due to the level of measurement noise in the time-series, 10mM and 25mM assays for NADH and 12mM assay for CoQ are not included in this analysis). To obtain a better estimate of the parameters, multiple time-series are used simultaneously instead in a combined least-square fitting across various initial concentrations. Even though individually the error value may increase between the simulation time-series and the experimental time-series used for the fitting, the estimated K R,D and K R,A parameter values from the combined fitting produce simulations that agree with the complex's characteristics across a much greater range of conditions. This is shown in Figure 11C and 11D, where the same time-series from Figure 11A are simultaneously used in a combined fitting of the parameters, and the subsequent simulations of the rate law based on this combined fitting are then shown to give a much better agreement with the time-series from the CoQ assays.
Note that the ratios D r=o and A o=r would approach infinity when the divisor concentrations (½D ox and ½A red ) approach zero (refer to Equation 27). Thus to avoid numerical errors in the simulation, a small amount of ½D ox and ½A red are assumed to have been created during the transient period before reaching the steady-state.

Comparison with Existing ETC Energy Transduction Equations
In Table 3 the new chemiosmotic rate law is compared to the ETC rate equations from both the Korzeniewski and the Beard OXPHOS models using complex I as the example. Each of the three rate equations is derived using a thermodynamics approach with the same formulation for the thermodynamic force (Equation 51) associated with the overall reaction of a given complex (Equation 1). However, in spite of their common origins, the three approaches differ in the level of detail at which the complex is modeled, which in turn results in functional differences.
Korzeniewski's rate equation follows Onsager's linear force-toflux relationship, based on the local equilibrium approximation   (2) and (5) is introduced in the enzyme transition mechanism that short-circuits the normal cycle of the enzyme and uncouples the free energy transduction to drive the proton translocation. doi:10.1371/journal.pone.0014820.g008 where thermodynamic forces vary slowly [48]. In contrast, Beard's rate equation and the chemiosmotic rate law work with non-linear steady-states that are not bound to the same restrictions. In fact, the derivation of Beard's equation follows the same thermodynamic force-to-flux relationships up to Equations 7 and 8, but diverges in the representation of the forward and reverse fluxes in Equation 4, where it simply assumes a mass-action reaction. By modeling the overall reaction as an elementary mass-action reaction, Beard's equation does not account for the kinetics within the complex. Its single kinetic parameter, denoted the activity parameter x, is a scaling parameter that serves to adjust the magnitude of the net flux (similar to ½DE epp A in Equation 32), but does not capture the maximum turnover rate of the enzyme or its affinity towards substrates. In contrast, the V f ' max and V r' max of the chemiosmotic rate law represent the apparent maximum velocities of the reaction, which are modulated by both thermodynamic and kinetic factors through the three bounded functions in Equation 46 ( Figure 9).
Least-square fitting of the three rate equations over the entire range of the experimental time-series data from an isolated complex I kinetics assay allow a comparison of how well the equations can reproduce the original time-series through their respective parameters. Figure 12A shows that the chemiosmotic rate law gives a better approximation of the experimental data compared to Korzeniewski's linear function, or Beard's massaction based function. However, a more accurate comparison of the principles guiding the three rate equations is given by fitting them to the initial velocity of the time-series data, where the three rate equations should theoretically converge. Figure 12B shows that the output of the chemiosmotic rate equation, as well as the original experimental data, fall between the other two approaches. When the comparison is extended to all possible values of K R,D and K R,A in the chemiosmotic equation (Figure 12C), one can note that the output of the chemiosmotic equation becomes equivalent to the output of Beard's equation when both K R,D and K R,A are equal to one, while it approaches the output of Korzeniewski's equation when both K R,D and K R,A approach zero. This can be explained by noting that when both K R,D and K R,A are equal to one in Equation 25, there is a direct correlation between the external concentrations and the internal probabilities, and Equation 19 becomes essentially the same as the Beard's firstorder mass-action reaction representation of complex I. On the other hand, when both K R,D and K R,A approach zero, all probabilities in Equation 26 approach one. This results in a zeroorder reaction that is insensitive to changes in the external concentrations, which is similar to the simulation output of Korzeniewski's linear flux equation.
An important result emerging from these comparisons is that the new chemiosmotic rate law represents a more general formulation of the kinetic and thermodynamic behaviors for a chemiosmotic ETC complex; one which encompasses both Korzeniewski's and Beard's formulations, and can smoothly interpolate between them to cover a spectrum of biochemical behaviors. Therefore, for a given overall reaction and a thermodynamic force definition, the chemiosmotic rate law can replace the existing ETC rate equations in the OXPHOS models of Korzeniewski and Beard to incorporate biochemically relevant kinetic parameters which allows a more accurate specification of an ETC complex.

Kinetic Parameter Values and Sensitivity Analysis
This section covers the analysis of the kinetic parameter values and their sensitivity in the chemiosmotic rate law applied to each of the ETC complexes in isolation. A summary of the experimental values obtained for all the kinetic parameters of complex I, III, and IV are given in Table 4, and the corresponding sensitivity in Table 5. In addition, the sensitivity of the slippage parameter c with equal perturbations on the forward and reverse rate constants (f s~0 :5) is also included for comparison. The sensitivity of all parameters are determined by calculating the magnitude of variations needed to change the flux of an isolated complex v by 1%. Thus, the lower the values in Table 5, the more sensitive the flux is to changes in the corresponding parameter.
In summary, the sensitivity of V f ' max provides the least amount of information about a specific complex since it is constant across different complexes regardless of its parameter value or concentrations of reaction species. This is due to the fact that V f ' max is independent from all other parameter and variable values in Equation 46. Variations dK S,D and dK S,A depend only on their respective parameters K S,D and K S,A , and total substrate concentrations D T and A T . The smaller the value of K S,D or  K S,A is compared to its respective total substrate concentration, the larger the parameter variation is required to change the flux, and therefore the less sensitive the parameter is. Compared to the other parameters in Table 5, where only small variations dK S,D , dK S,A , and dc are needed to affect the flux, multiple-fold changes of K R,D and K R,A are necessary in order to affect the flux. These large values of dK R,D and dK R,A suggest that the flux is not very sensitive to changes in K R,D or K R,A . However, as shown previously in Figure 12C, small changes in the values of K R,D and K R,A can give rise to significant changes in the curvature of the simulated reaction time-series. This is due to the fact that D r=o and A o=r are time-dependent ratios of the species concentrations, which tend to infinity when there are only substrate species present, or become zero when there are only product species left. Thus, in the beginning of a time-series, when the values of D r=o and A o=r are large, K R,D and K R,A are small in comparison, and their variations have negligible effect on the flux. However, as D r=o and A o=r get smaller and closer to the values of K R,D and K R,A over time, their effect on the flux would become much more significant. In addition, the larger the value of the varying parameter, the earlier the reaction flux is affected by the concentration of the substrate (Figure 13). On the other hand, the smaller the value of the varying parameter, the longer the reaction maintains the same flux, resulting in a sharper change in the flux at the end of the time-series, when the reaction runs out of its reactant(s). Whereas a K R,A or K R,D value of zero makes the reaction independent from the respective reactant concentration and a value of one makes the reaction dependent on the reactant concentration in a first-order mass-action fashion, a value larger than one suggests an even higher dependence on the reactant concentration than first-order. In addition, from Equation 24, when either K R,D or K R,A is larger than one, the redox potential of the respective boundary electron transfer reaction is necessarily negative, indicating an energetically unfavorable reaction. These derived relationships are important for the interpretation of the experimental values of K R,A and K R,D in Table 4. For complex I, both K R,D and K R,A are lower than one. In contrast, for both complex III and complex IV, the cyctochrome c associated K R,D or K R,A value is larger than one, suggesting a higher-order   dependence on the concentration of cytochrome c consistently with the fact that two cyctochrome c are required for the corresponding reaction. Thus, even though complex III and IV are not very sensitive to the saturation concentration of cytochrome c, the redox state ratio of cytochrome c has a profound effect on the flux through the redox state function R. Lastly, from Equation 55, it is clear that the effect of dc depends largely on the redox potential of the system ( Table 2).

Network Sensitivity Analysis
In this section, we conduct a more global sensitivity and perturbation analysis with the same experimentally observable in vitro homogenate environment, but with all the ETC complexes functioning in tandem. The simple pathway model of the ETC used for this analysis consists of interactions between a driving dehydrogenase reaction, the three main electron transfer complexes (complex I, III, and IV), and the three electron carrier  redox pairs ( Figure 1A). Each complex is modeled using the chemiosmotic rate law applied with the parameter values derived from our experiments (Table 4), together with thermodynamic constants from the literature and saturation reactant concentrations used in our experiments (Table 2).
To analyze the sensitivity of a network, the most commonly used tool is the flux control coefficient (FCC) [49], which represents the relative change in the global steady state flux J resulting from an infinitesimal change in a property of an individual complex i, divided by the relative change of that complex's activity v i from the same infinitesimal change, and normalized by the corresponding steady state flux and complex activity. The complex with the largest FCC exerts the largest control on the flux at a particular steady state, as an increase in the activity of this complex would result in the largest overall flux increase. For the present study, an expanded definition of FCC is used such that the properties of interest are the individual parameters j i,j of the ith complex: This expanded definition expresses quantitatively the effect that small variations in the parameter j i,j have on the flux of the system J, if the effect of j i,j on the local complex activity v i is known. The FCCs are calculated for the simple ETC pathway by taking the same parameter variation as in Table 5 to find the corresponding changes in J numerically, and then dividing the result by v i . The results are summarized in Table 6, which shows that regardless of   Figure 12C close to the zero-order linear approximation. B. The curvature ( dJ dt ) time-series obtained from the second-derivative of the concentration time-series in Panel A show that the closer the concentration time-series is to the linear approximation, the larger the changes to the curvature time-series. C. The concentration time-series with respect to changes in K R,D while K R,A is held at 0.3 show that variations in K R,D produce a subset of the time-series close to the first-order mass action approximation. D. The curvature timeseries obtained from the second-derivative of the concentration time-series in Panel C show that the closer the concentration time-series is to the mass-action approximation, the smaller the changes to the curvature time-series, and also the lower the peak time of the curvature. doi:10.1371/journal.pone.0014820.g013 which parameter is varied, the control coefficients remain consistent for the same complex. This indicates that for the same small variation in v i , the parameters have the same effect on J without significantly affecting the steady state concentrations.
Although FCCs can show the relative control of the individual components in a pathway network, they have little predictive value as each FCC is computed for a single steady-state. As the state of the system changes, the values of the FCCs, may change as well. To overcome the single steady-state limitation of the FCC, we next apply threshold curve analysis [13] to study the effects of different parameter perturbations on our simplified ETC network across a continuum of steady-states. Compared to the FCC, a threshold curve is not limited to just infinitesimal changes in a parameter j i,j , but gives a complete range of variations that result in 0% to 100% decrease in both v i and J. The percent decreases in v i and J are plotted versus each other, and the relationship curve between them gives not only a measure of the global flux control by an individual enzyme complex, but possibly also a threshold value for v i beyond which the level of J would be significantly reduced.
An example is shown in Figure 14A, where variations in the V f ' max parameter of complex I have a direct one-to-one effect on the complex I activity v 1 , but a delayed effect on J. The two curves are combined to form the complex I threshold curve in Figure 14B, where it is paired with the threshold curves of complex III and IV to produce the specific threshold profile of the system. The threshold effects shown in these curves can be attributed to a combination of: (1) an excess of complex activity due to an excess of complex available in the system, which tends to produce a curve with a plateau phase, followed by a sharp decline in J; and (2) the buffering of individual complex activity perturbations by the metabolic network (kinetic properties of the enzymes, structure of the pathway network, concentrations of substrates, etc.), which is responsible for the smoothness of the curves [50]. As such, the clear threshold value of complex I, and the smooth threshold curves of complex III and IV indicate that complex I is in excess relative to complex III and IV, which affect J in a more gradual and controlling manner. Thus, the threshold profile gives the same conclusion as the FCC analysis, but offers a more complete view of a network's interactivity across the entire normalized range of J steady-states.   The same threshold profile can be observed irrespective of which kinetic parameter is used to establish the relationship between that of v i and J, so long as the two are always compared at the same steady-state. Figure 15 shows that this is true for the complex I threshold curve with respect to modulations in the V f ' max parameter, the saturation parameters (K S,D or K S,A ), and the redox state parameters (K R,D or K R,A ). However, in the case of the c parameter, the threshold is markedly changed because the Figure 16. Parameter Perturbations and the Threshold Curve. Each parameter in the rate law perturbs the complex I threshold curve and threshold profile in a different way. A. The saturation parameters (represented by K S,D ) do not alter the metabolic network significantly, which allows the threshold curve to maintain its shape and a threshold value. B. At K S,D~1 00, there are no perceivable effect on the threshold curves of the other two complexes. C. The redox state parameters (represented by K R,D ) have a significant effect on the metabolic network, and thus the shape of the threshold curve as it changes v 1 's order of dependency on the reactant concentrations. D. At K R,D~1 :0, the entire threshold profile is greatly affected by the changes in the metabolic network. E. The c parameter has no effect on the threshold curve until its value drops below 0.1. After which, the threshold curve quickly approaches that of a straight line. F. At c~0:01, v 1 becomes the sole rate-limiting reaction. doi:10.1371/journal.pone.0014820.g016 thermodynamic property of the complex, and the structure of the network are changed with the introduction of an alternative pathway for the electron in the complex.
The threshold profile of the system is consistent for all kinetic parameters because the threshold curve relates the relative changes instead of the absolute changes in v 1 and J. Thus, threshold curves based on different kinetic parameters cannot show how each parameter may alter the network. Instead, we look at how the threshold profile of the system is changed by a single value perturbation in a specific parameter. This is demonstrated in Figure 16, which compares how perturbations in the saturation parameters, redox state parameters, and the c parameter affect the complex I threshold curve and the threshold profile of the system. A visual metrics of comparing the threshold between the threshold curves is presented as the area between two consecutive curves, which represents the difference in the area under the two curves (color coded also by the difference in the color index of the two curves). The quantitative measure of each threshold curve in threshold profile is presented in Table 7 as the fraction of the total square plot area covered under the curve. From Figure 16A, increases in the value of a saturation parameter (represented by K S,D ) shifts the threshold value of the curve to the left but the plateau phase is maintained for a large range of parameter values. The perturbed threshold profile in Figure 16B shows that even at a value of 100mM, the effect of a saturation parameter is localized to the threshold of its own complex and has no significant effect on the threshold of the other two complexes. In contrast, Figure 16C shows that increases in the value of a redox state parameter (represented by K R,D ) has a significant effect on the shape and smoothness of the threshold curve as it changes v 1 's order of dependency on the reactant concentrations, indicative of a large perturbation effect on the metabolic network. At K R,D~1 :0, the perturbed threshold profile in Figure 16D shows that this perturbation in the metabolic network is propagated to complex III and IV, and affects their threshold curves greatly. The c parameter has essentially no effect on the threshold curve in Figure 16E until its value drops below 0.1. This suggests that the thermodynamic force of complex I, in the absence of the pmf, has an excess in the redox potential, and that only ten percent of the redox potential is necessary to drive the redox reaction of complex I. After dropping below 0.1, the threshold curve quickly approaches that of a straight line, showing direct effect on the flux. The corresponding threshold profile in Figure 16D shows that, at c~0:01, v 1 becomes so rate-limiting that v 3 and v 4 of complex III and IV are always in excess in relation to v 1 .

Conclusion
In short, we have developed a framework for modeling mitochondria bioenergetics using a new chemosmotic rate law to represent each ETC complex in the OXPHOS pathway. This modular framework subsumes and generalizes several previous approaches and relies on a small set of biochemically relevant parameters. We have conducted enzymatic assays to derive those kinetic parameters for complex I, III, and IV, and validated the predictions of the model against experimental concentration timeseries. These results, together with detailed sensitivity analyses, show that the parameters, originally derived from a simple reference system, provide a good breadth of feasible physiological responses of the complexes. In particular, threshold curves relating the flux of one complex to the global flux of the system show how each parameter in each complex can have a differential effect on the threshold curve of the corresponding complex, as well as on the overall threshold profile of the system ( Figure 14B and 16). The flexibility and accuracy of the model, coupled with the diverse range of behaviors it is capable of generating through transformation of the system, suggest that in the future this approach could enable the comparative modeling and analysis of mitochondria from different systems and pathological states. Table 7. Measure of Threshold in Threshold Profile of Figure 16 (Fraction of the Total Area Under the Curve).