Kinetic Mathematical Modeling of Oxidative Phosphorylation in Cardiomyocyte Mitochondria

Oxidative phosphorylation (OXPHOS) is an oxygen-dependent process that consumes catabolized nutrients to produce adenosine triphosphate (ATP) to drive energy-dependent biological processes such as excitation-contraction coupling in cardiomyocytes. In addition to in vivo and in vitro experiments, in silico models are valuable for investigating the underlying mechanisms of OXPHOS and predicting its consequences in both physiological and pathological conditions. Here, we compare several prominent kinetic models of OXPHOS in cardiomyocytes. We examine how their mathematical expressions were derived, how their parameters were obtained, the conditions of their experimental counterparts, and the predictions they generated. We aim to explore the general landscape of energy production mechanisms in cardiomyocytes for future in silico models.


Introduction
The metabolic process of oxidative phosphorylation (OXPHOS) consumes chemical energy from catabolism to generate adenosine triphosphate (ATP), the energy currency of the cell. OXPHOS complexes sit on the inner mitochondrial membrane (IMM), consisting of electron transport chain (ETC) complexes I, II, III, and IV, as well as ATP synthase (complex V) [1] (Figure 1). The ETC complexes transfer electrons from reducing equivalents (NADH, FADH2) to oxygen molecules in controlled steps and extract the released chemical energy to pump protons out of the mitochondrial matrix, generating a proton electrochemical gradient called the proton motive force (pmf) that drives ATP synthase to produce ATP [2][3][4]. Cells can couple ATP hydrolysis to energy-consuming biological processes, e.g., muscle contraction, ion homeostasis, and macromolecule synthesis. OXPHOS energy generation is more efficient than that of anaerobic glycolysis. Mitochondrial OX-PHOS produces approximately 30 ATP molecules for one glucose molecule, compared to two ATP molecules in glycolysis [5]. Thus, mitochondria are famous as "the powerhouse of the cell [6]".
Oxygen is an excellent electron acceptor, with a redox potential of~800 mV per electron under physiological conditions [7]. Nonetheless, some oxygen molecules are not reduced completely to water in the ETC. Instead, they become reactive oxygen species (ROS): superoxide, hydrogen peroxide, and hydroxyl radicals [8]. These highly reactive molecules can serve as a redox signaling mechanism from the mitochondria under physiological conditions [9][10][11]. However, excessive ROS can damage biomolecules, such as mitochondrial DNA (mtDNA), mitochondrial proteins for metabolism (e.g., aconitase and complex I), and the peroxidation of membrane lipids crucial for OXPHOS function [12,13]. The ROS scavenging systems include superoxide dismutase (SOD), cytochrome c, catalase, glutathione (GSH), and thioredoxin (Trx) systems [14,15]. heart failure, ischemic-reperfusion injury, diabetic cardiomyopathy, and chemotherapyinduced cardiomyopathy [38]. Investigating the complex interactions and inner workings of OXPHOS and related metabolic systems using systems biology approaches and computer simulation models can complement the knowledge gathered from in vitro and in vivo studies, identify candidate mechanisms, and corroborate/propose novel hypotheses [39][40][41][42]. For these reasons, we review the mathematical modeling of OXPHOS complexes, focusing on cardiomyocyte bioenergetics models.

Law of Mass Action
The simplest way to describe a biochemical reaction rate is to assume a sufficiently large number of particles spreading evenly in a solution. Since the reaction rate is proportional to the number of particles colliding with each other, it is also proportional to the product of substrate concentrations to the power of their stoichiometric coefficients. For an elemental reversible reaction aA + bB = cP + dQ, the reaction rate in the PQ-forming direction can be described as: where the forward rate constant ( ) is often an adjustable parameter determined by the literature and model fitting; the backward rate constant is constrained by the equilibrium Mitochondrial dysfunction is involved in a variety of pathological conditions, such as type 2 diabetes [16][17][18][19][20][21], liver toxicity [22,23], neurodegenerative diseases [24][25][26][27][28], arrhythmia, and heart failure [29][30][31][32][33][34]. In this review, we focus on the kinetic modeling of ATP generation in cardiomyocyte mitochondria because (i) many computational models and experimental data on cardiomyocyte mitochondria are available; (ii) cardiomyocytes rely on mitochondria-derived ATP for action potential firing, muscle contraction, and ion homeostasis in each heartbeat [35]. Mitochondria can take up to 30% of the cardiomyocyte volume [36]. However, cardiomyocytes maintain metabolite stability under varying workloads [37]; (iii) the OXPHOS system is interconnected with other metabolic pathways, such as the citric acid cycle, mitochondrial calcium handling, metabolite transfer, and ROS generation and scavenging. The interactions between these systems are related to cardiomyocyte pathologies such as altered calcium dynamics and electrophysiology, heart failure, ischemic-reperfusion injury, diabetic cardiomyopathy, and chemotherapy-induced cardiomyopathy [38]. Investigating the complex interactions and inner workings of OX-PHOS and related metabolic systems using systems biology approaches and computer simulation models can complement the knowledge gathered from in vitro and in vivo studies, identify candidate mechanisms, and corroborate/propose novel hypotheses [39][40][41][42]. For these reasons, we review the mathematical modeling of OXPHOS complexes, focusing on cardiomyocyte bioenergetics models.

Law of Mass Action
The simplest way to describe a biochemical reaction rate is to assume a sufficiently large number of particles spreading evenly in a solution. Since the reaction rate is proportional to the number of particles colliding with each other, it is also proportional to the product of substrate concentrations to the power of their stoichiometric coefficients. For an elemental reversible reaction aA + bB = cP + dQ, the reaction rate in the PQ-forming direction can be described as: where the forward rate constant (k f ) is often an adjustable parameter determined by the literature and model fitting; the backward rate constant is constrained by the equilibrium constant (K eq ). The reactants (A, B, P, and Q) can be substrates or different stages of the enzyme catalyzing the reaction.
The law of mass action is the basis of the famous Michaelis-Menten enzyme kinetics and more complex schemes.

Quasi-Steady State Enzyme Kinetics
Michaelis-Menten kinetic expressions are derived from the law of mass action in enzyme catalysis reactions. Take an irreversible reaction S → P catalyzed by enzyme E as an example, assuming that the first substrate binding reaction is reversible and the second product formation reaction is irreversible.
Suppose we treat the intermediate enzyme-substrate complex (ES) as the fast species and use the conservation relationship (Et = ES + E = constant) to constrain the free enzyme (E). In that case, the quasi-steady-state assumption reduces the two reactions into one by time-scale separation. The product formation rate is then described as a Michaelis-Menten hyperbolic function of the concentration of substrate S [43][44][45][46][47], where V max is the maximal reaction rate when the enzyme is saturated, and K m is the Michaelis constant, the concentration of S when the reaction rate reaches half of V max .
For enzyme reactions involving more than one substrate/product, more combinations, such as a random binding order, a compulsory binding order, ping-pong mechanisms, and dead-end intermediates, are possible [47,48]. The number of parameters and the complexity of the mathematical description increase as more binding and conversion scenarios are considered [43,49]. For simplicity, the enzyme and enzyme-substrate complexes are often treated as fast species, and their transitions are assumed to be in a quasi-steady-state (the rates of change of the transition states are set to zero). The state occupancies and the enzyme turnover rate can be solved systematically by linear algebra or by the Hill-King-Altman diagram method [50,51].
For example, for a reversible enzyme-catalyzed reaction, S + E ES EP E + P the reaction scheme can be drawn as in Figure 2. kij are rate constants for the enzyme state i transitioning to state j. After applying quasi-steady state approximation for the enzyme stages, the state occupancies of each stage can be solved by linear algebra.
Alternatively, the state occupancies can be solved by the King-Altman-Hill diagram method [50]. This graph-based method first reduces the edges of the diagram to be a set of simple connected graphs without cycles. For example, the reaction cycle of Figure 2a can be broken down into three graphs: Λ, <, and > [50]. The weight of each stage (w i ) is the sum of the products of the transition rates indicated in Figure 2b. Alternatively, the state occupancies can be solved by the King-Altman-Hill diagram method [50]. This graph-based method first reduces the edges of the diagram to be a set of simple connected graphs without cycles. For example, the reaction cycle of Figure 2a can be broken down into three graphs: Λ, <, and > [50]. The weight of each stage ( ) is the sum of the products of the transition rates indicated in Figure 2b.
The quasi-steady state reaction flux is then the product of the turnover rate (v) and the enzyme activity. For enzymes with more states and state transitions, a computer program called KAPattern [51] can automatically derive the weights of each state and the overall reaction rate.

Lumped Models of the Electron Transport Chain
A simplistic approach to describe reaction rates in the ETC is to treat the ETC complexes (I, II, III, and IV) as a single-unit proton pump.
Korzeniewski's OXPHOS models [52,53] assumed that complex I and III reactions The quasi-steady state reaction flux is then the product of the turnover rate (v) and the enzyme activity. For enzymes with more states and state transitions, a computer program called KAPattern [51] can automatically derive the weights of each state and the overall reaction rate.

Lumped Models of the Electron Transport Chain
A simplistic approach to describe reaction rates in the ETC is to treat the ETC complexes (I, II, III, and IV) as a single-unit proton pump.
Korzeniewski's OXPHOS models [52,53] assumed that complex I and III reactions were close to equilibrium. The overall reaction rate of NADH reducing cytochrome c linearly depended on the thermodynamic span (∆G 13 ), the deviation between the current condition and the equilibrium. ∆G 13 was determined by the ratios of NAD/NADH and oxidized/reduced cytochrome c couples, the mitochondrial membrane potential, and the proton motive force.
Cells 2022, 11, 4020 [54] built a lumped complex I-III-IV ETC model ( Figure 3) based on the six-state proton pump model by Pietrobon and Caplan [55] to study calcium dynamics and bioenergetics in mouse pancreatic beta-cell mitochondria. This model showed that the rates of proton translocation and oxygen consumption depended on the ratio of NADH to NAD and the mitochondrial membrane potential (∆Ψ m ).
where J HR is the proton translocation rate, J O is the oxygen consumption rate, V T is the thermal voltage (~26.7 millivolts at 37 • C), ρ res is the concentration of the respiratory complex, ∆Ψ B is the bulk phase boundary potential, g is the correction factor for ∆Ψ m , and r a , r b , r c1 , r c2 , r 1 , r 2 , and r 3 are constant parameters for the ETC. This lumped model did not include the ROS generation mechanism or the redox states of ubiquinone and cytochrome c pools.  [55] to study calcium dynamics and bioenergetics in mouse pancreatic beta-cell mitochondria. This model showed that the rates of proton translocation and oxygen consumption depended on the ratio of NADH to NAD and the mitochondrial membrane potential ( ). where is the proton translocation rate, is the oxygen consumption rate, is the thermal voltage (~26.7 millivolts at 37 °C), is the concentration of the respiratory complex, is the bulk phase boundary potential, g is the correction factor for , and , , , , , , and are constant parameters for the ETC. This lumped model did not include the ROS generation mechanism or the redox states of ubiquinone and cytochrome c pools.  [54]. The six-state proton pump consists of proton translocation (5 → 6 → 1 → 2), NADH oxidation (2 → 3 → 4), and a slip reaction (5 → 2). Each cycle represents one electron transfer from NADH to oxygen, coupled with the translocation of 6 protons from the matrix to the intermembrane space (IMS). The arrows represent the forward reaction direction, although all reactions are reversible. Reprinted/adapted with permission from Ref. [54].  [54]. The six-state proton pump consists of proton translocation (5 → 6 → 1 → 2), NADH oxidation (2 → 3 → 4), and a slip reaction (5 → 2). Each cycle represents one electron transfer from NADH to oxygen, coupled with the translocation of 6 protons from the matrix to the intermembrane space (IMS). The arrows represent the forward reaction direction, although all reactions are reversible. Reprinted/adapted with permission from Ref. [54]. The mitochondrial energetics model by Cortassa et al. [56] adapted Magnus and Keizer's model [54] and replaced the dependence on the mitochondrial membrane potential (∆Ψ m ) with the proton motive force (pmf, ∆p). The model added succinate oxidation flux via complexes II, III, and IV, depending on the pmf and the FADH2 to FAD ratio.
where V He is the proton translocation rate from complex I-III-IV, V He(F) is the proton translocation rate from complex II-III-IV, V O 2 is the oxygen consumption rate, ∆p is the pmf, ρ res(F) is the concentration of complex II-III-IV, V T is the thermal voltage, ρ res is the concentration of the respiratory complex, ∆Ψ B is the bulk phase boundary potential, g is the correction factor for ∆Ψ m , and r a , r b , r c1 , r c2 , r 1 , r 2 , and r 3 are constant parameters for the ETC.

Complex I
Mitochondrial complex I, also called NADH: ubiquinone oxidoreductase or type I NADH dehydrogenase, is the largest complex (with 46 subunits in mammalian mitochondria) in the ETC [57][58][59]. Its primary function is to transfer two electrons from one reduced nicotinamide adenine dinucleotide (NADH) molecule to one ubiquinone (coenzyme Q, denoted as Q) and pump four protons from the mitochondrial matrix to the intermembrane space. The chemical reaction is summarized as follows: The core molecular structure of complex I, sitting in the IMM, consists of one hydrophilic peripheral arm, protruding toward the mitochondrial matrix, and one hydrophobic membrane arm [58]. The peripheral arm conducts electron transfer from NADH to ubiquinone and consists of the NADH dehydrogenase module, the chain of iron-sulfur (FeS) clusters, and the ubiquinone reductase module. The membrane arm contains four parallel proton pumps; each translocates one proton once the bound ubiquinone is fully reduced and triggers a conformational change in the membrane arm. The redox process in complex I is reversible. Under a high pmf and a highly reduced Q pool with a large QH2 to Q ratio, the direction of the reaction favors QH2 oxidation, NAD+ reduction, and proton translocation from the IMS to the matrix. This is known as reverse electron transport (RET) [10]. Conversely, forward electron transport (FET) means that electrons flow 'normally' from complex I to ubiquinone.
The forward electron flow in complex I is NADH → FMN → FeS cluster chain [57][58][59]. The NAD+/NADH couple can transfer only two electrons simultaneously, while the FeS cluster chain can handle only Cells 2022, 11, 4020 7 of 37 one electron at a time. Thus, FMN acts as a buffer between the NAD+/NADH couple and the FeS cluster chain because FMN can switch between three redox states: fully oxidized FMN, FMNH radical, and fully reduced FMNH2 [60]. Likewise, the ubiquinone molecule bound to the ubiquinone reductase receives electrons from the FeS clusters to form the ubisemiquinone radical (Q), and finally to the ubiquinol molecule (QH2). The rate of electron transport from NADH to FeS cluster N2 is much faster than the overall turnover rate of complex I (a time scale of 100 µs vs. 5 ms), with the FeS clusters mostly in the reduced state under physiological conditions, indicating that ubiquinone binding, reduction, and release, are rate-limiting steps [58]. The molecular structure of complex I also implies that the conformational changes induced by ubiquinone reduction couple the translocation of four protons from the matrix to the IMS [57].

Reactive Oxygen Species (ROS) Generation Mechanisms
Complex I is considered one of the major ROS producers in the mitochondria [12]. Several possible ROS-producing mechanisms are supported by experimental evidence in complex I [58,[61][62][63]. Most studies recognize that the flavin site (IF) in the NADH dehydrogenase module and the quinone site (IQ) in the quinone reductase module are two major sources of ROS generation in complex I [64]. These two sites have different characteristics. ROS generation in the IF site is proportional to fully reduced flavin (FADH2) not bound by NAD/NADH [63]. The FeS cluster N1a can reduce the lifetime of superoxideforming flavosemiquinone radical (FMNH•) [65]. ROS generation at the IF site is enhanced by blocking electron transport at the IQ site, for instance, by rotenone. ROS generation in the IQ site, on the other hand, is dependent on the fraction of bound semiquinone radicals, which is small during FET, but significant during RET. The inhibition of quinone reductase by rotenone decreases ROS generation during RET. Thus, whether rotenone enhances or inhibits ROS generation in complex I depends on several factors, including different substrates (glutamate/malate vs. succinate), the pmf, and the redox state of the ubiquinone pool, all of which determine the direction of complex I reaction. Kinetic modeling can help us elucidate the complicated dependence of the reaction rate on respiring substrates and the pmf in complex I.

Kinetic Models of Complex I
In Korzeniewski's rate skeletal muscle OXPHOS models [53,66], the complex I reaction rate (J C1 ) depended hyperbolically on its thermodynamic span (∆E C1 ), the deviation of current condition from thermodynamic equilibrium.
Later, in a rat skeletal muscle OXPHOS model [67] and a rat cardiomyocyte OXPHOS model [68] showed a linear relationship between the complex I reaction rate (J C1 ) and its thermodynamic span (∆E C1 ).
The thermodynamic span (∆E C1 ) can be derived from substrate (NAD/NADH and ubiquinone/ubiquinol couples) concentrations and the pmf (∆p). Beard's OXPHOS model [69] also assumed that complex I operated near-equilibrium under physiological conditions and used a simple law of mass action descriptor for complex I reaction flux (J C 1 ).
where the enzyme activity (X C 1 ) is an adjustable parameter fitted by experimental data of porcine cardiomyocyte mitochondria [70], and the equilibrium constant (K C 1 eq ) depends on the proton concentration, mitochondrial membrane potential (∆Ψ m ), and Gibbs free energy change (∆ r G 0 C 1 ) of the reaction. The descriptors for ROS generation were not included. Heiske et al. [71] built an OXPHOS model, which described the complex I reaction rate (v C 1 ) based on Henry-Michaelis-Menten kinetics.
where k f is the forward reaction constant, k b is the backward reaction constant, and K x m is the apparent Michaelis constant for species x (NADH, NAD, QH2, or Q). The reaction constants and apparent Michaelis constants are thermodynamically constrained by the change in the Gibbs free energy in the reaction. Markevich and Hoek's mitochondrial ROS model [63] also included a complex I model using the law of mass action for its 17 reactions. Reactions of NAD/NADH and Q/QH2 couples, the IF (flavin) [60] and IQ (Q-reduction) site, and the FeS clusters of complex I were included to identify the source of ROS formation. They concluded that the majority of ROS production in complex I came from the IF site during FET and the IQ site during RET.
Gauthier et al. [62,72] built a seven-state complex I model ( Figure 4) based on the Magnus-Keizer [54,55] six state proton pump model. The enzyme went through a cycle of NADH oxidation, ubiquinone reduction, and proton translocation (state 1→ 2 → 3 → 4 → 7 → 5 → 6 →1). A superoxide-generating pathway (state 4 → 2) was included. where the enzyme activity ( ) is an adjustable parameter fitted by experimental data of porcine cardiomyocyte mitochondria [70], and the equilibrium constant ( ) depends on the proton concentration, mitochondrial membrane potential ( ), and Gibbs free energy change ( ) of the reaction. The descriptors for ROS generation were not included. Heiske et al. [71] built an OXPHOS model, which described the complex I reaction rate ( ) based on Henry-Michaelis-Menten kinetics.
where is the forward reaction constant, is the backward reaction constant, and is the apparent Michaelis constant for species x (NADH, NAD, QH2, or Q). The reaction constants and apparent Michaelis constants are thermodynamically constrained by the change in the Gibbs free energy in the reaction.
Markevich and Hoek's mitochondrial ROS model [63] also included a complex I model using the law of mass action for its 17 reactions. Reactions of NAD/NADH and Q/QH2 couples, the IF (flavin) [60] and IQ (Q-reduction) site, and the FeS clusters of complex I were included to identify the source of ROS formation. They concluded that the majority of ROS production in complex I came from the IF site during FET and the IQ site during RET.
Gauthier et al. [62,72] built a seven-state complex I model ( Figure 4) based on the Magnus-Keizer [54,55] six state proton pump model. The enzyme went through a cycle of NADH oxidation, ubiquinone reduction, and proton translocation (state 1→ 2 → 3 → 4 → 7 → 5 → 6 →1). A superoxide-generating pathway (state 4 → 2) was included.  [72]. The seven-state proton pump is adapted from the six-state proton pump model of Magnus and Keizer [54,55]. The model includes proton pumping (5 → 6 → 1→ 2), NADH oxidation (2 → 3 → 4), ubiquinone reduction (4 → 7 → 5), and superoxide generation (4 → 2). Each cycle represents one electron transfer from NADH to ubiquinone or oxygen. The arrows represent the forward reaction direction, although all reactions are reversible. Reprinted/adapted with permission from Ref. [62].  [72]. The seven-state proton pump is adapted from the six-state proton pump model of Magnus and Keizer [54,55]. The model includes proton pumping (5 → 6 → 1→ 2), NADH oxidation (2 → 3 → 4), ubiquinone reduction (4 → 7 → 5), and superoxide generation (4 → 2). Each cycle represents one electron transfer from NADH to ubiquinone or oxygen. The arrows represent the forward reaction direction, although all reactions are reversible. Reprinted/adapted with permission from Ref. [62]. The complex I model by Bazil et al. [61] includes five states by the number of electron(s) the complex holds ( Figure 5). They assumed that the electron transfer inside complex I was much faster than reactions with substrates; electrons were distributed in steady-state across the redox centers: FMN, FeS cluster N2, and the Q-site. To eliminate internal state variables, the state transitions were also assumed to be in a quasi-steady state. The enzyme turnover and ROS generation rates were functions of the pmf and substrate concentrations. The kinetic parameters were fitted by the kinetic data and the ROS data from bovine heart mitochondria and constrained by thermodynamics in redox chemistry. They reproduced the reaction rates between NAD/NADH and Q/QH2 couples and the ROS production rates of superoxide and hydrogen peroxide. The complex I model by Bazil et al. [61] includes five states by the number of electron(s) the complex holds ( Figure 5). They assumed that the electron transfer inside complex I was much faster than reactions with substrates; electrons were distributed in steadystate across the redox centers: FMN, FeS cluster N2, and the Q-site. To eliminate internal state variables, the state transitions were also assumed to be in a quasi-steady state. The enzyme turnover and ROS generation rates were functions of the pmf and substrate concentrations. The kinetic parameters were fitted by the kinetic data and the ROS data from bovine heart mitochondria and constrained by thermodynamics in redox chemistry. They reproduced the reaction rates between NAD/NADH and Q/QH2 couples and the ROS production rates of superoxide and hydrogen peroxide.  [61]. The 5 states represent the number of electrons held by complex I. The model includes peroxide generation reactions (green arrows) taking two electrons from complex I to oxygen, superoxide generation reactions (orange arrows) taking one electron from complex I to oxygen, NADH oxidation reactions (blue arrows) taking two electrons from NADH to complex I, and finally ubiquinone reduction reactions (black arrows) taking two electrons from complex I to ubiquinone, coupled with 4-proton translocation. The figure was reprinted/adapted with permission from Ref. [61].

Molecular Structure and Reaction Mechanism
Mitochondrial complex II is also called succinate dehydrogenase (SDH), or succinatecoenzyme Q reductase (SQR). Complex II couples the citric acid cycle (CAC) with the ETC by transferring electrons from succinate to ubiquinone (Q), forming fumarate and ubiquinol (QH2) [73]: The molecular structure of complex II consists of four subunits: SDHA (succinate dehydrogenase site, with FAD as the prosthetic group), SDHB (FeS clusters), SDHC and SDHD (ubiquinone reductase site) [74]. The electron transfer inside complex II is summarized as follows: In the forward direction, succinate is oxidized to fumarate, and Q is reduced to QH2. Nevertheless, the reaction is reversible, depending on the concentrations of succinate, fumarate, Q, and QH2. Significant succinate accumulation is found in ischemic cells and the strong forward reaction upon reoxygenation is responsible for RET-related ROS generation and ischemic-reperfusion injury [9,10,75]. Complex II activity and ROS generation from RET are inhibited by oxaloacetate and malonate [76]. Figure 5. Schematic of the complex I model by Bazil et al. [61]. The 5 states represent the number of electrons held by complex I. The model includes peroxide generation reactions (green arrows) taking two electrons from complex I to oxygen, superoxide generation reactions (orange arrows) taking one electron from complex I to oxygen, NADH oxidation reactions (blue arrows) taking two electrons from NADH to complex I, and finally ubiquinone reduction reactions (black arrows) taking two electrons from complex I to ubiquinone, coupled with 4-proton translocation. The figure was reprinted/adapted with permission from Ref. [61].

Molecular Structure and Reaction Mechanism
Mitochondrial complex II is also called succinate dehydrogenase (SDH), or succinatecoenzyme Q reductase (SQR). Complex II couples the citric acid cycle (CAC) with the ETC by transferring electrons from succinate to ubiquinone (Q), forming fumarate and ubiquinol (QH 2 ) [73]: The molecular structure of complex II consists of four subunits: SDHA (succinate dehydrogenase site, with FAD as the prosthetic group), SDHB (FeS clusters), SDHC and SDHD (ubiquinone reductase site) [74]. The electron transfer inside complex II is summarized as follows: In the forward direction, succinate is oxidized to fumarate, and Q is reduced to QH 2 . Nevertheless, the reaction is reversible, depending on the concentrations of succinate, fumarate, Q, and QH 2 . Significant succinate accumulation is found in ischemic cells and the strong forward reaction upon reoxygenation is responsible for RET-related ROS generation and ischemic-reperfusion injury [9,10,75]. Complex II activity and ROS generation from RET are inhibited by oxaloacetate and malonate [76].

Kinetic Models of Complex II
In the guinea pig cardiomyocyte mitochondrial model by Cortassa et al. [56], SDH irreversibly oxidized succinate (SUC) into fumarate, and was inhibited by oxaloacetate (OAA) and fumarate (FUM). However, in this model, succinate oxidation was not coupled to proton pumping in the ETC.
where k cat is the catalytic constant, E T is the SDH concentration, K m is the Michaelis constant for succinate, K OAA i is the inhibition constant for OAA, and K FUM i is the inhibition constant for fumarate.
Wei's acid-base equilibria model [77] (also used in [14]) coupled ETC and CAC reactions by assuming that the FADH2 to FAD ratio in the Cortassa et al. model [56] was determined by the succinate to fumarate ratio. The reaction rate was also inhibited by OAA. However, because complexes II, III, and IV were lumped, there were no reactions of the ubiquinone/ubiquinol couple.
where V He(F) is the proton translocation rate, V O 2 SDH is the oxygen consumption rate, ∆p is the pmf, ρ res(F) is the concentration of complex II-III-IV, K SDH res is the equilibrium constant of succinate oxidation, P SUC is the binding polynomial for unprotonated succinate, K i is the inhibition constant for OAA, ∆Ψ B is the bulk phase boundary potential, g is the correction factor for ∆Ψ m , and r a , r b , r c1 , r c2 , r 1 , r 2 , and r 3 are constant parameters for the ETC (complexes II, III, and IV).
Demin's OXPHOS models [62,78,79] assumed succinate excess. The reaction rate of SDH depended hyperbolically on the fraction of oxidized ubiquinone in the Q pool.
The complex II formulation in the hepatocyte mitochondria model by Mogilev skaya et al. [80] used random order rapid equilibrium Bi-Bi enzyme kinetics. The model was used to study salicylate toxicity to hepatocytes by inhibiting CAC activity.
Manhas et al. [84] devised a five-state complex II model to study how substrates and inhibitors affect its reaction rate and ROS generation ( Figure 6). The authors assumed that electron transport within complex II is much faster than substrate binding (succinate, fumarate, Q/QH2 couple, and oxygen). The five states represent how many electrons are held by complex II, ranging from zero to five; the state transitions represent reactions between complex II and its substrates. The thermodynamically constrained parameters were fitted against various datasets from bovine, pig, and guinea pig heart mitochondria experiments. The authors reproduced the result of succinate oxidation kinetics with and without inhibitors (atpenin and malonate) and ROS production rates from complex II. The authors concluded that the primary source of ROS is the [3Fe-4S] iron-sulfur complex under most physiological conditions. Manhas et al. [84] devised a five-state complex II model to study how substrates and inhibitors affect its reaction rate and ROS generation ( Figure 6). The authors assumed that electron transport within complex II is much faster than substrate binding (succinate, fumarate, Q/QH2 couple, and oxygen). The five states represent how many electrons are held by complex II, ranging from zero to five; the state transitions represent reactions between complex II and its substrates. The thermodynamically constrained parameters were fitted against various datasets from bovine, pig, and guinea pig heart mitochondria experiments. The authors reproduced the result of succinate oxidation kinetics with and without inhibitors (atpenin and malonate) and ROS production rates from complex II. The authors concluded that the primary source of ROS is the [3Fe-4S] iron-sulfur complex under most physiological conditions.  [84]. The 5 states represent the number of electrons in complex II. The model includes peroxide generation reactions (green arrows) taking two electrons from complex II to oxygen, superoxide generation reactions (orange arrows) taking one electron from complex II to oxygen, succinate oxidation reactions (blue arrows) taking two electrons from NADH to complex II, and finally ubiquinone reduction reactions (black arrows) taking one electron from complex II to ubiquinone. The figure was reprinted/adapted from Ref. [84].
Markevich et al. [85] developed a complex II model using mass action kinetics, which included substrate binding and electron transport processes with a total of 31 reactions. The parameters were fitted against experimental data from bovine and rat heart mitochondria. The authors also investigated the difference between intact and disintegrated complex II, where SDHA/SDHB subunits were separated from SDHC/SDHD subunits.  [84]. The 5 states represent the number of electrons in complex II. The model includes peroxide generation reactions (green arrows) taking two electrons from complex II to oxygen, superoxide generation reactions (orange arrows) taking one electron from complex II to oxygen, succinate oxidation reactions (blue arrows) taking two electrons from NADH to complex II, and finally ubiquinone reduction reactions (black arrows) taking one electron from complex II to ubiquinone. The figure was reprinted/adapted from Ref. [84].
Markevich et al. [85] developed a complex II model using mass action kinetics, which included substrate binding and electron transport processes with a total of 31 reactions. The parameters were fitted against experimental data from bovine and rat heart mitochondria. The authors also investigated the difference between intact and disintegrated complex II, where SDHA/SDHB subunits were separated from SDHC/SDHD subunits. The authors concluded that in disintegrated complex II, Q-binding site inhibition, and complex III inhibition had a high-amplitude bell-shaped dependence on succinate in ROS generation. In contrast, intact complex II under physiological conditions had a hyperbolic dependence on succinate in ROS generation.

Complex III
Respiratory complex III is also called the cytochrome bc1 complex, coenzyme Q: cytochrome c-oxidoreductase. It transfers electrons from ubiquinol (QH 2 ) to cytochrome c (cytc) by the following reaction:

. Molecular Structure and Reaction Mechanism
Complex III translocates protons outward via Mitchell's Q-cycle [86]. The first and rate-limiting reaction is ubiquinol oxidation to ubiquinone at the Qo (outer) site. The two electrons from ubiquinol bifurcate toward two paths. One electron goes through the high potential pathway (Qo → iron-sulfur protein (ISP) → cyt c1 → cyt c) to reduce cytochrome c; the other goes through the low potential pathway (Qo → cyt bL → cyt bH → Qi) to reduce ubiquinone at the Qi (inner) site. Two ubiquinol molecules oxidized at the Qo site complete a full cycle and release four protons to the intermembrane space. Meanwhile, one ubiquinone molecule at the Qi site is reduced to ubiquinol, absorbing two protons from the mitochondrial matrix. Overall, the Q-cycle translocates two charges for every ubiquinol oxidized at the Qo site with two electrons transferred [87].
On the inner mitochondrial membrane, complex III monomers form dimers, allowing electron transfer between their two cyt bL centers. However, electron transfer between monomers has been considered nonsignificant; the dimer acts as two independent monomers under physiological conditions [88][89][90].

Reactive Oxygen Species Generation Mechanisms
In addition to complex I, complex III is another major source of mitochondrial-derived reactive oxygen species. (mtROS) Semiquinone radicals (SQ) at the Qo site can reduce oxygen to form superoxide. There are two mechanisms for SQ formation, the semiforward and the semireverse mechanisms. The semiforward mechanism states that after ubiquinol (QH2) at the Qo site donates one electron to the iron-sulfur protein (ISP), the remaining SQ intermediate reduces oxygen to form superoxide. [89] In contrast, the semireverse mechanism argues that the SQ intermediate is too unstable to stay and play a significant role in superoxide generation. Instead, QH2 donates all its electrons downstream, one to ISP and the other to bL. [89] Superoxide forms when bL is highly reduced due to high pmf or the downstream pathway is blocked by inhibitors such as antimycin A. An electron can jump from reduced bL to a ubiquinone molecule at the Qo site, forming a transient SQ intermediate, which reduces oxygen to superoxide [87,93]. Both mechanisms of superoxide generation are strongly dependent on the pmf, but the semireverse mechanism implies that the rate of superoxide generation is maximized when the ubiquinone pool is moderately reduced, which is more consistent with experimental findings on antimycin-induced ROS generation from complex III [87]. Complex III releases superoxide to both sides of the IMM [94].

Kinetic Models of Complex III
In Korzeniewski's rate skeletal muscle OXPHOS models [53,66], the complex III reaction rate (J C3 ) depended hyperbolically on its thermodynamic span (∆E C3 ), the deviation of current condition from thermodynamic equilibrium.
A rat skeletal muscle OXPHOS model [67] and a rat cardiomyocyte OXPHOS model [68] by Korzeniewski et al. showed a linear relationship between the complex I reaction rate (J C3 ) and its thermodynamic span (∆E C3 ).
The thermodynamic span (∆E C3 ) can be derived from substrate (ubiquinone/ubiquinol and oxidized/reduced cytochrome c couples) concentrations, the mitochondrial membrane potential, and the pH difference across the IMM (embedded in the pmf, ∆p).
Beard's OXPHOS model [69] assumed that the complex III reaction operated near equilibrium and used the law of mass action for the reaction rate. Activation from inorganic phosphate ([Pi] m ) was added to explain the data in a mitochondrial experiment [70]. However, a later study [95] found that complex III activity was not enhanced directly by inorganic phosphate. The apparent activation by inorganic phosphate was likely due to associated substrate transport and nonsteady-state experimental conditions.
Heiske et al. [71] built an OXPHOS model, where they described the complex III reaction rate (J c3 ) based on Michaelis-Menten kinetics. The apparent activation by inorganic phosphate was inherited from Beard's OXPHOS model [69].
where k f is the forward reaction constant, k b is the backward reaction constant, and K x m is the apparent Michaelis constant for species x (ubiquinone, ubiquinol, reduced and oxidized cytochrome c). The reaction constants and apparent Michaelis constants are thermodynamically constrained by the change in Gibbs free energy in the reaction. Demin et al. [79] modeled Mitchell's Q-cycle using the law of mass action. The system included ubiquinone diffusion on the IMM and electron transfer between ubiquinone, complex III redox centers (ISP, cyt c1, bL, and bH), and cyt c. The model used the semiforward mechanism for ROS generation; in the model, superoxide was generated from the oxygen oxidizing Qo site ubisemiquinone radical. The model consisted of two variants: the minimal and channeled models. The minimal model treated the Qo site ubisemiquinone radical as part of the Q pool, not binding to complex III. In the channeled model, on the other hand, Qo site ubisemiquinone radicals were bound to complex III, next to either ISP or heme bL. The superoxide generation rate showed exponential dependence on the mitochondrial membrane potential in the minimal model but sigmoid dependence in the channeled model. The minimal model was later integrated into the redox balance model by Gauthier et al. [62].
Demin et al. [78] also built a newer model by coupling the redox center states of ISP, cyt bL, and cyt bH to ensure electron bifurcation. The semiquinone radical at the Qo site remained bound until bL was oxidized by bH. The locking mechanism prevented both electrons from ubiquinol from going through ISP and bypassing the Q-cycle.
Selivanov's complex III model [96] included all possible combinations of binding and redox states, including ubiquinone (at both the p-side and n-side), cytochromes bL and bH, ISP, and cytochrome c1, creating a 400-state ODE model. They found that there are two steady-state branches in complex III in simulations and experiments. Increasing the succinate supply or reducing the oxygen concentration switched complex III from a low ROS-generating mode to a high ROS-generating mode, which could explain the increased ROS production from ischemia-reperfusion.
Guillaud et al. [97] constructed a 16-state complex III model in which superoxide was produced via the semireverse mechanism by reduced cytochrome bL and oxidized ubiquinone. Mass action kinetics described the transition between the redox states of ISP, cyt bL, cyt bH, and the Qi site. The model reproduced the bell-shaped dependence of ROS generation on the fraction of reduced quinone in the experiments.
Markevich and Hoek [63] developed a complex III model using mass action kinetics to test competing ROS generation mechanisms under a variety of proton motive forces, substrate levels, and inhibitors. They concluded that the late dissociation of ISP, that is, ISP leaving cyt bL after ubiquinone, is more consistent with experimental data.
Bazil et al. [98] developed an algebraic descriptor of complex III to reproduce experimental results of electron transport, ROS generation, and bistability. The model used a six-state scheme to represent the number of electrons on the Qo and Qi sites. Kinetic rate constants were determined by experimental data, pH, the pmf, and changes in the Gibbs free energy for state transitions. The net turnover flux was derived from state occupancies at steady state, obtained by solving a linear system of equations. Later, a complex III functional dimer model [99] was built upon a similar principle (Figure 7). The dimer model had six-states representing zero to five electrons shared by the redox centers. The dimer model not only reproduced bistability in electron transfer and ROS generation but also showed that low-dose antimycin partially blocking the Qi site stimulates cytochrome c reduction, which can only be explained via a dimer model without using a different parameter set. The dimer model also revealed that ROS were generated predominantly via the semireverse mechanism by reduced cytochrome bL under most conditions. ROS generation from the semiforward mechanism (by semiquinone radical) became significant only with a highly reduced ubiquinone pool under antimycin blocking the Qi site. rameter set. The dimer model also revealed that ROS were generated pre the semireverse mechanism by reduced cytochrome bL under most cond eration from the semiforward mechanism (by semiquinone radical) bec only with a highly reduced ubiquinone pool under antimycin blocking th

Complex IV
Respiratory complex IV, also called cytochrome c oxidase (CcO), is the terminal enzyme of the ETC, receiving electrons from cytochrome c and catalyzing the reduction of oxygen to water [100]. For every oxygen molecule reduced, four cytochrome c molecules are oxidized, and eight positive charges are translocated outward against the mitochondrial membrane potential.

Molecular Structure and Reaction Mechanism
In complex IV, electrons flow from cytochrome c → CuA → cytochrome a → cytochrome a3-CuB binuclear center (BNC) to oxygen. Because every cytochrome c molecule carries only one electron, it requires four reduced cytochrome c molecules to fully reduce one oxygen molecule to water.
The active site of complex IV to oxygen is the BNC and the catalytic cycle consists of two phases: the eu-oxidase phase and the peroxidase phase [100]. In the eu-oxidase phase, an oxygen molecule binds doubly reduced BNC and is reduced to peroxide. In the peroxidase phase, BNC receives two electrons, fully reducing peroxide to water. Under physiological conditions, complex IV operates continuously and each electron transfer from cytochrome c is coupled to the outward translocation of one proton [100][101][102]. However, slip reactions (electron transfer without proton pumping) have been reported in isolated complex IV [103], with two protons instead of four being pumped per oxygen molecule reduced. It is also hypothesized that under a high pmf, complex IV trades proton pumping efficiency for a higher enzyme turnover via the slip pathway.
The BNC also evens out the redox potential changes for each electron transfer, keeping the midpoint potential for oxygen reduction at approximately −800 mV per electron [102]. The tight interactions between BNC and oxygen ligands also minimize ROS generation during the reaction [12,104].
Complex IV activity is inhibited by nitric oxide (NO) [105,106], cyanide (CN − ) [107], carbon monoxide (CO) [108], hydrogen sulfide (H 2 S) [109], azide (N 3 − ) [110], and formate [111], the toxicity of which can block the electron transport chain and ATP synthesis. Among complex IV inhibitors, nitric oxide (NO) is the most physiologically relevant because under physiological conditions, a substantial fraction (up to half) of complex IV activity is inhibited by NO [105,106]. NO inhibits complex IV in two modes: competitive inhibition with oxygen binding and uncompetitive inhibition in which NO is reduced to nitrite. NO binding to complex IV is also sensitive to near infrared and red photons, causing the dissociation of NO, and freeing up available complex IV for OXPHOS [106,112].

Kinetic Models of Complex IV
The early OXPHOS model by Korzeniewski and Froncisz [52] included a three-state (A1, A2, and A3) complex IV reaction cycle [113] that depended on the phosphorylation potential.
The reaction rate depended on the ratio of oxidized and reduced cytochrome c couple, the oxygen concentration, and the phosphorylation potential term (AAP), which is the quotient of ATP by ADP and inorganic phosphate concentrations.
In Korzeniewski's later OXPHOS models for hepatocytes [114] and myocytes [53], the descriptors for complex IV were simplified and dependent on the pmf. The reaction rate depended on the oxygen concentration, the reduced cytochrome c concentration (c 2+ ), and the availability of reduced cytochrome a3 (a 2+ ). a 2+ was assumed in rapid equilibrium and was derived from the ratio of oxidized and reduced cytochrome c couple, the mitochondrial membrane potential (∆Ψ), the pH difference across the inner mitochondrial membrane (embedded in the pmf, ∆p), and the midpoint potentials of cytochrome c (E mc ) and cytochrome a3 (E ma ) in the binuclear center.
Beard's OXPHOS model [69] used mass action kinetics for cytochrome c and oxygen and included a hyperbolic dependence for the oxygen concentration.
The reaction rate of complex IV (J C 4 ) depended on the enzyme activity X C 4 , the oxygen concentration ([O 2 ]) and affinity (K O 2 ), the cytochrome c pool size ([∑ cytc]), the mitochondrial membrane potential (∆Ψ m ), and the apparent equilibrium constant of complex IV (K C 4 eq ), which depended on the proton concentration, the ∆Ψ m , and the changes in Gibbs free energy (∆ r G 0 C 4 ). Heiske et al. [71] built an OXPHOS model that described the complex IV reaction rate (v c4 ) based on Henry-Michaelis-Menten kinetics.
where k f is the forward reaction constant, k b is the backward reaction constant, and K x m is the apparent Michaelis constant for species x (oxygen, reduced and oxidized cytochrome c). The reaction constants and apparent Michaelis constants are thermodynamically constrained by the change in the Gibbs free energy.
Demin et al. [78] built a reaction cycle of four states (Y, Yr, YO, and YOH) for complex IV with mass action kinetics as part of their OXPHOS model (Figure 8) cally constrained by the change in the Gibbs free energy.
Demin et al. [78] built a reaction cycle of four states (Y, Yr, YO, and YOH) IV with mass action kinetics as part of their OXPHOS model (Figure 8). The mi redox balance model by Gauthier et al. [62] used the King-Altman-Hill diagr on this complex IV model to eliminate intermediate state variables of Y, Yr, YO Wilson et al. [116] modeled complex IV around the redox states of the bin copper center. Between state transitions, the equilibrium constants came fro point potentials of oxygen and binuclear center reduction reactions; the kinet were fitted from experiments on isolated mitochondria. The steady-state o sumption rate came from mass action kinetics. The authors showed that the re were highly dependent on the energy state Q, acidity, and redox state of the c pool. They also revealed that the apparent Michaelis constant for oxygen the energy state Q increases, similar to the results by Krab et al. [115].
Pannala et al. [117] improved the complex IV model by Wilson et al. [1 verting the arbitrary energy state Q to the proton motive force (1 Q = 1.5 p competitive and uncompetitive inhibition by nitric oxide (NO), adding a proto mechanism coupled with each cyt c oxidation, assuming all steps are reversible to strict thermodynamic constraints, and deriving steady-state reaction rat Wilson et al. [116] modeled complex IV around the redox states of the binuclear ironcopper center. Between state transitions, the equilibrium constants came from the midpoint potentials of oxygen and binuclear center reduction reactions; the kinetic constants were fitted from experiments on isolated mitochondria. The steady-state oxygen consumption rate came from mass action kinetics. The authors showed that the reaction rates were highly dependent on the energy state Q, acidity, and redox state of the cytochrome c pool. They also revealed that the apparent Michaelis constant for oxygen increases as the energy state Q increases, similar to the results by Krab et al. [115].
Pannala et al. [117] improved the complex IV model by Wilson et al. [116], by converting the arbitrary energy state Q to the proton motive force (1 Q = 1.5 pmf), adding competitive and uncompetitive inhibition by nitric oxide (NO), adding a proton pumping mechanism coupled with each cyt c oxidation, assuming all steps are reversible to conform to strict thermodynamic constraints, and deriving steady-state reaction rates from the King-Altman-Hill diagram method. This model reproduced turnover rates under a range of pmf and cyt c redox states, as well as NO inhibition under a range of oxygen levels and cyt c redox states.
where n is the number of protons translocated per ATP, also called the H+/ATP ratio, which is approximately three in yeast and mammalian mitochondria [121,122].

Molecular Structure and Reaction Mechanism
Mitochondrial complex V is a membrane-bound molecular motor on the inner mitochondrial membrane (IMM). Complex V can be separated by water solubility into a soluble globular F 1 catalytic sector in the matrix (subunit α, β, γ, δ, and ε) and a membranebound F 0 proton-translocating sector (subunit c-ring, a, b, d, e, f, g, A6L, and oligomycin sensitivity-conferring protein (OSCP)). Mechanically, complex V can be divided into the rotor part (subunit c-ring, γ, δ, and ε) and the stator part (the remainder) [120,123,124].
Complex V acts like a turbine generator. The proton gradient across the inner mitochondrial membrane rotates the turbine (oligomeric c ring) and the stalk (gamma subunit) connected to the turbine [124]. In the active site of alpha3/beta3, MgADP and Pi are collected and assembled into MgATP. The mechanical interactions between the rotating stalk and the stator subunits α3/β3 induce a conformational change, altering the binding affinity of MgADP/Pi/MgATP substrates. MgATP is then "pushed" out of complex V, completing the catalytic cycle [120].
Complex V proteins often form a row of V-shaped dimers on the apex of mitochondrial cristae to enhance the efficiency of ATP synthesis. Because proton-pumping ETC complexes concentrate on mitochondrial cristae, the proton gradient is the greatest across the cristae interior and the intermembrane space, allowing efficient ATP synthesis [124][125][126].
Complex V can be inhibited by a range of substances. Oligomycin blocks the proton channels in the c-ring and is commonly used in mitochondrial metabolism studies [123].
Recent studies by Juhaszova et al. [130,131] demonstrated that complex V can also harness the "potassium motive force" to generate ATP in experiments on proteoliposome (PL)-reconstituted purified ATP synthase and isolated mitochondria. The potassium ion influx is driven by the mitochondrial membrane potential and the potassium concentration gradient.
Potassium ions entering the mitochondrial matrix are translocated back out by the mitochondrial K+/H+ exchanger (mKHE), with a 1:1 stoichiometry of protons entering the mitochondrial matrix [132,133].
Therefore, these study results are fully compatible with Mitchell's chemiosmotic mechanism.

Kinetic Models of Complex V
In Korzeniewski's rate skeletal muscle OXPHOS models [53,66,67], the ATP synthesis rate (J C 5 ) was nonlinearly dependent on the thermodynamic force (γ) of complex V. The thermodynamic span was the displacement from the thermodynamic equilibrium and was determined by the pmf (∆p), the ratio of ADP and inorganic phosphate to ATP, and the free energy change of ATP synthesis (∆G p0 ). The H:ATP ratio (n A ) was 8/3.
Beard's OXPHOS model [69] assumed that complex V operated close to equilibrium and utilized mass action kinetics for its substrates: ATP, ADP, and inorganic phosphate.
where X F 1 is the activity of complex V, and K F 1 eq is the apparent equilibrium constant depending on the proton gradient, binding polynomials of substrates (P ATP , P ADP , P PI ), change in Gibb's free energy (∆ r G 0 F 1 ), and the mitochondrial membrane potential (∆Ψ m ). Heiske et al. [71] built an OXPHOS model that described the complex V reaction rate (J c5 ) based on Henry-Michaelis-Menten kinetics.
where k f is the forward reaction constant, k b is the backwards reaction constant, and K x m is the apparent Michaelis constant for species x (ATP, ADP, inorganic phosphate). The reaction constants and apparent Michaelis constants are thermodynamically constrained by the change in the Gibbs free energy.
The beta cell mitochondria model by Magnus and Keizer [54] assumed that complex V was a six-state proton pump operating in reverse, using the same framework as the ETC in the same model. The ATP synthesis rate and the proton consumption rate strongly depended on both the mitochondrial membrane potential and the phosphorylation potential. The H + /ATP stoichiometry was higher than 3:1 due to the presence of slip reactions.
where V ATPase is the ATP generation rate, V Hu is the proton flux through complex V, V T is the thermal voltage, ρ F1 is the concentration of complex V, K F1 eq is the apparent equilibrium constant, ∆Ψ B is the phase boundary potential, and p a , p b , p c1 , p c2 , p 1 , p 2 , and p 3 are the parameters of complex V. The mitochondrial energetics model for guinea pig ventricular cardiomyocytes by Cortassa et al. [56] replaced the dependence on the mitochondrial membrane potential with the dependence on the pmf in Magnus and Keizer's model (Figure 9). A factor of 100 was also added to address reversibility under low pmf. Later models [134][135][136][137][138][139] also utilized this schematic for mitochondrial energetics in guinea pig ventricular cardiomyocytes. REVIEW 20 was also added to address reversibility under low pmf. Later models [134][135][136][137][138][139] als lized this schematic for mitochondrial energetics in guinea pig ventricular cardio cytes.  [54] and Wei [77]. In addition to the proton flux ( ), the model also considered the potassium (   Figure 9. Schematic of the complex V model by Magnus and Keizer [54], assuming complex V is a 6-state proton pump operating in reverse. Not accounting for slip reactions (5 → 2), one ATP molecule is generated for every 3 protons translocated from the intermembrane space (IMS) to the matrix. The figure was reprinted/adapted from Magnus and Keizer [54]. The complex V model by Cortassa et al. in 2022 [140] incorporated ATP generation from the potassium gradient, based on the studies of Juhaszova et al. [130,131]. The mathematical descriptions were based on the model of Magnus and Keizer [54] and Wei et al. [77]. In addition to the proton flux (J H C 5 ), the model also considered the potassium (J K C 5 ) and the sodium fluxes (J Na C 5 ) through complex V. All three fluxes contributed to ATP synthesis (J ATP C 5 ).
The complex V descriptors in Nguyen et al.'s cardiomyocyte OXPHOS model [141] used Hill kinetics for the mitochondrial membrane potential and random-order Bi-Uni kinetics for ADP, ATP, and inorganic phosphate (Pi). The calcium activation of ATP synthesis was also considered. The complex V reaction was assumed to be fully coupled; thus, the H:ATP ratio was exactly 3:1.
where V C 5 is the ATP synthesis rate, V max is the activity of complex V, K ATP , K ADP , and K Pi are the dissociation constants for ATP, ADP, and inorganic phosphate, respectively; K v is the voltage activation factor, and K Ca is the calcium activation factor. The complex V model by Metelkin et al. [79,142] to study the interactions between complex V and the adenine nucleotide translocator (ANT) used a similar scheme and included the concentrations of protons on both sides of the IMM.

Summaries of the OXPHOS Models
In 1967 Chance proposed the first OXPHOS model using operational flux expressions [143], Bohnensack [144] and Holzhutter [145] later constructed the first kinetic OX-PHOS model and used thermodynamic force to drive the OXPHOS reactions [146]. These studies have established the foundation of kinetic modeling of OXPHOS, and more detailed and sophisticated models have been proposed since then [146]. In this section, we summarize the mitochondrial OXPHOS computational modeling mentioned in the previous sections in Table 1. Because many of the studies included two or more respiratory complexes, we provide overall remarks in this section rather than in the respective OXPHOS complex sections above.
The rat hepatocyte [52,114] and skeletal muscle cell OXPHOS models [53,67] by Korzeniewski et al. included OXPHOS complexes I, III, IV, and V, along with proton leak, substrate dehydrogenation, adenine nucleotide translocator (ANT), phosphate carrier (PiC), and adenylate kinase (AdK). Cytochrome c reduction (by complex I and complex III) and ATP synthesis (by complex V) rates were linked to their thermodynamic spans, which is the displacement from thermodynamic equilibria. Complex I, III, and V reactions were assumed to operate near equilibrium [147]. Cytochrome c oxidation (by complex IV) used mass action kinetics. The oxygen-binding reduced heme a3 was assumed to be in rapid equilibrium and derived by oxidized to reduced cytochrome c ratio and the pmf. Korzeniewski et al. used these models to study the metabolic flux control of ATP generation and consumption in different metabolic modes (glycolysis, fatty acid oxidation, and gluconeogenesis) in hepatocytes [114,148], metabolic flux control and the impact of respiratory complex deficiencies in skeletal muscle cells [53], comparing the negative feedback and parallel activation mechanisms of ATP synthesis and consumption in skeletal muscle cells [66,67,149], and the regulation of OXPHOS reactions in cardiomyocytes [68]. They demonstrated that parallel activation to NADH generation, proton pumping by respiratory complexes, ATP synthesis, and ATP consumption could explain why the ATP to ADP ratio barely changes when the work intensity increases. The OXPHOS models by Korzeniewski et al. are simple, thermodynamically consistent (for complexes I, III, V, and ANT), and in good agreement with experimental data of ATP synthesis, proton leak, and oxygen consumption. However, they did not include detailed descriptions of the citric acid cycle, calcium dynamics, and ROS generation/scavenging mechanisms.
Saito et al. [150] built a cardiomyocyte mitochondrial model, which integrated complex I, III, IV, and V models by Korzeniewski et al. [52,53,67,114], citric acid cycle reactions, ion (proton, sodium, calcium, potassium, and phosphate) and respiring substrate (malate and glutamate) dynamics, and high energy phosphate buffering by adenylate kinase and creatine kinase. The authors used this model to study why the metabolite levels of NADH and ATP remain relatively stable under various workloads. They found that inorganic phosphate (Pi) stimulated ATP synthesis more effectively than calcium.
Beard's OXPHOS model [69] included complexes I, III, IV, and V, as well as ANT, PiC, AdK, potassium leakage, and potassium efflux through mKHE. The reaction rates of the OXPHOS complexes were described by mass action kinetics, assuming they operated near equilibrium. The reaction rate constants were also constrained by the change of Gibbs free energy in each complex and fitted by the Bose et al. dataset [70] of porcine heart and skeletal muscle mitochondria. However, this mitochondrial model did not include explicit citric acid cycle processes, ROS generation/scavenging, or mitochondrial calcium dynamics. This model required complex III to be explicitly activated by inorganic phosphate to fit the experimental data. The simple mass action kinetics descriptors did not follow saturation kinetics in enzyme-catalyzed reactions.
Wu et al. [82] extended Beard's OXPHOS model [69] and integrated citric acid cycle reactions (including complex II, succinate dehydrogenase) and substrate exchange across the IMM. The parameters were estimated based on the Bose et al. dataset [70] of porcine heart and skeletal muscle mitochondria and the LaNoue et al. dataset [151] of rat heart mitochondria. The authors concluded that ADP and NAD regulated the citric acid cycle fluxes. In contrast, ATP synthesis was regulated by inorganic phosphate, which activated complexes III and V. However, ROS generation/scavenging and calcium dynamics were not included in this model. Bazil et al. [83] extended the mitochondrial model by Wu et al. [82] and integrated sodium, calcium, and potassium influx/efflux reactions as well as mitochondrial matrix volume dynamics due to osmolarity changes caused by the matrix potassium concentration.
Heiske et al. [71] built an OXPHOS model based on the OXPHOS models of Beard [69] and Korzeniewski and Zoladz [67]. The authors used Henry-Michaelis-Menten kinetics instead of mass action kinetics for the rate equations of the OXPHOS complexes (I, III, IV, V) to address saturation in enzyme kinetics. The parameters were constrained by the changes in the Gibbs free energy and fitted by their experimental data [71,152] and the Bose et al. dataset [70]. However, this model did not include citric acid cycle reactions, ROS generation and scavenging, and calcium dynamics.
The complex III monomer [98] and dimer [99] [84] have used quasisteady state assumptions to describe the overall reaction rates of substrate consumption and ROS formation. Each model had several redox states representing the number of electrons in the redox centers. Redox reactions between the complex and substrates were state transitions. The kinetic parameters were constrained by Gibbs free energy changes of the reactions and fitted by various datasets, primarily from bovine cardiomyocyte mitochondria experiments. By applying the quasi-steady state assumption and setting the rate of change of each state to zero, the overall reaction rates could be solved from a system of linear equations of states and their transition rates. These respiratory complex models have reproduced redox reaction rates under various conditions. They have been used to study the source of ROS generation in ischemia/reperfusion (I/R) injury [81,153] and to test the hypothesis of calcium phosphate inhibiting complex I during calcium overload in a computational-experimental approach [154].
Markevich et al. built complex I [63], III [63], and II [85] models using mass action kinetics. These models have been used to study ROS generation mechanisms and steady-state ROS production rates in various substrate inputs and a range of mitochondrial membrane potentials. Because their mass action kinetics models included detailed reaction steps of substrate binding, proton pumping, enzyme state transitions, and ROS generation, the authors can test different ROS generation mechanisms and simulate complex inhibition by changing the kinetic constants. The minimal mitochondrial model by Magnus and Keizer [54] included a lumped respiratory chain (complex I-III-IV) and ATP synthase (complex V) for the OXPHOS system as well as adenine ANT, calcium influx via the mitochondrial calcium uniporter (MCU), and calcium efflux via the mitochondrial sodium-calcium exchanger (mNCX). Matrix proton, inorganic phosphate, and sodium concentrations were assumed to be constant, so their dynamics were not included. ROS production and scavenging systems were also absent. The OXPHOS system used the six-state proton pump scheme from Pietrobon and Caplan [55]. The (quasi-)steady-state rates of NADH oxidation, proton translocation, and ATP synthesis were derived by the King-Altman-Hill diagram method [50]. Magnus and Keizer later incorporated this mitochondrial model with cellular glycolysis and ion channels to simulate the mouse beta-cell membrane potential, ATP, and calcium oscillations in the presence of glucose [155,156]. The mitochondrial energetics model by Cortassa et al. [56] integrated the calcium dynamics and OXPHOS from the Magnus-Keizer model [54] and the citric acid cycle (CAC) reactions from the Jafri-Dudycha model [157]. The CAC provided NADH for the ETC and was stimulated by mitochondrial calcium. The model was used to simulate the generation and consumption of NADH and ATP in cardiomyocyte mitochondria under various workloads, demonstrating that mitochondrial calcium can enhance energy supply. However, the succinate oxidation in the respiratory chain did not couple with the succinate dehydrogenase reaction in the CAC. Later the same group integrated the mitochondrial model with electrophysiology and muscle contraction systems and constructed the excitationcontraction coupling/mitochondrial energetics (ECME) model of guinea pig ventricular cardiomyocytes [136]. The simulations were in good agreement with the experimental results of force-frequency relationships, mitochondrial NADH, and mitochondrial calcium changes in response to varying workloads.
Bertram et al. [158] simplified Cortassa et al.'s model [56], known as the BLPS model, to simulate mitochondrial ATP synthesis. The electron transport chain and ATP synthase were described by Michaelis-Menten kinetics for the substrates (NADH, ADP) and sigmoid dependence on the mitochondrial membrane potential. Saa and Siqueira [159] improved and corrected the BLPS model to study the metabolic response under a range of glycolysis fluxes and oscillatory calcium levels.
Nguyen et al. [141] built a cardiomyocyte mitochondrial model that integrated the ETC model from Magnus and Keizer [54], the CAC reactions from the Jafri-Dudycha model [157], and the calcium uniporter and ATP synthase model from the authors' previous model [160]. The authors added NCX, NHX, and PiC to simulate inorganic phosphate, sodium, and proton dynamics. The models can reproduce mitochondrial calcium concentrations and influx under periodic stimulation, consistent with cardiomyocyte experiments. They also showed that increasing cytosolic sodium levels enhanced calcium extrusion from the mitochondria via the NCX, leading to decreased mitochondrial calcium concentrations and ATP synthesis.
Cortassa et al. [135] built a ROS-induced ROS release (RIRR) model based on their previous model [56]. The authors integrated CAC, OXPHOS, ROS production, ROSinducible anion channels for ROS efflux, and ROS scavenging by catalase, SOD, and the glutathione (GSH) system. Since there were no mechanistic descriptions for ROS production, it was modeled by taking an adjustable fraction of oxygen consumption to produce superoxide. They simulated the ROS-induced ROS release (RIRR) phenomenon and reproduced mitochondrial membrane potential depolarization and recovery events observed in guinea pig ventricular cardiomyocytes under oxidative stress [161].
Zhou et al. [138] incorporated the ECME model [136] and the RIRR model [135] to study action potential shortening and insensitivity to electrical stimuli due to oxidative stress-induced mitochondrial depolarization, ATP synthase reversal, decreased ATP to ADP ratio, and the opening of ATP-sensitive potassium channels. The ECME-RIRR model has been used to study potential oscillations in a mitochondrial network [137,162], arrhythmogenesis from blocked action potential propagation [139], and abnormal action potential patterns due to ROS disrupting calcium dynamics in the cytosol [163,164].
Wei et al. [77] extended the mitochondrial energetics model by Cortassa et al. [56] by integrating proton, phosphate, sodium, and calcium dynamics. The authors included mathematical descriptors for PiC, MCU, and mNCX. The authors also improved the complex II model by coupling succinate oxidation in the ETC and the CAC. The model results were in good agreement with mitochondrial membrane potential, mitochondrial NADH, and pH changes in guinea pig cardiomyocyte experiments. However, ROS production and scavenging mechanisms were not included in this model.
Kembro et al. [14] extended the mitochondrial model of Wei et al. [77] and the RIRR model of Cortassa et al. [135] to build a two-compartment (cytosol and mitochondria) ROS production, transport, and scavenging model. The mitochondrial energetic-redox (ME-R) model included mitochondrial NADPH generation from transhydrogenase (THD) and NADPH-producing isocitrate dehydrogenase (IDH2), ROS efflux via anion channels (for superoxide) and simple diffusion (for hydrogen peroxide), and ROS scavenging by catalase, SOD, and the GSH and Trx systems [15]. The fatty acid oxidation model by Cortassa et al. [165] added beta oxidation steps for palmitate to the two-compartment ROS generation, transfer, and scavenging model by Kembro et al. [14] The experimentalcomputational approach showed that ROS excess upon palmitate addition resulted from impaired GSH and Trx ROS scavenging systems and palmitate-induced uncoupling.
Cortassa et al. [140] extended the model of Kembro et al. [14] by including a complex V model utilizing the potassium gradient to generate ATP and an ETC model coupled to mitochondrial matrix volume dynamics, which was controlled by potassium content. The authors reproduced experimental data of changes in matrix volume, potassium uptake, and mitochondrial membrane potential from low to high energy demand.
Gauthier et al. [62] built a mitochondrial model, called the ETC-ROS model, to study the redox balance hypothesis [166], which stated that ROS production is minimal and energy generation is maximal when the cellular environment is moderately reduced. Thus, a mechanism of ROS generation from the respiratory complexes was needed. The complex I model was derived from the six-state proton pump scheme of Pietrobon and Caplan [55]. The complex I model included an oxygen-reducing pathway for superoxide production and an additional state for ubiquinone reduction. It reproduced the guinea pig cardiomyocyte mitochondrial data of Aon et al. [166]. The complex II, III, and IV models were modified from the hepatocyte mitochondrial model by Demin et al. [78,79] to reproduce the macrophage OXPHOS data from Kim et al. [167]. The complex V model was based on the mitochondrial model by Wei et al. [77]. The ROS scavenging systems included the SOD, GSH, and Trx systems. The authors showed that mitochondrial ROS production was minimized when the redox environment (NAD/NADH and NADP/NADPH couples) was moderately reduced, supporting the redox-optimized ROS balance (R-ORB) hypothesis [166].
Later Gauthier et al. reported a new mitochondria model [72] that incorporated the OXPHOS and ROS generation parts from the ETC-ROS model [62] and the ROS scavenging parts from the ME-R model [14] to study ROS production in heart failure. Their model can reproduce respiration rates, NADH to NAD ratios, and ROS production rates at state three (ADP present) and state four (ADP absent). The model also demonstrated that either blocking mitochondrial calcium uniporter or increasing cytosolic sodium concentrations (as in heart failure) suppressed mitochondrial calcium, NADH, NADPH, and GSH levels while ROS production increased. De Oliveira et al. [168] extended this model [72] and included doxorubicin-induced acute and chronic cardiac toxicity. Their model can simulate acute toxicity by inactivating OXPHOS complexes and enhancing ROS production. In contrast, chronic toxicity was simulated by ROS damaging mitochondrial DNA (mtDNA), inhibiting the long-term activity of OXPHOS complexes and leading to a vicious cycle as faulty OXPHOS complexes tend to produce more ROS and damage mtDNA even further. The authors predicted the dose response of both acute and chronic doxorubicin toxicity. They also showed that extending the duration of the iron chelator therapy partially mitigated chronic doxorubicin toxicity and protected cardiomyocytes against mitochondrial dysfunction.

Other OXPHOS Aspects and Concluding Remarks
In this review, we focus on the mathematical modeling of the OXPHOS complexes, how their reaction mechanisms have been described, and where and how ROS come from. Although we cover the descriptors from complex I to complex V (Table 1), to build an integrated mitochondrial bioenergetics model, additional aspects related to OXPHOS must be addressed.
For instance, ATP generation can be stimulated by mitochondrial calcium, ADP, and inorganic phosphate to meet varying workload demands. Mitochondrial calcium levels increase with cytosolic calcium levels upon excitation-contraction coupling and stimulate dehydrogenases in the CAC. ADP and inorganic phosphate are substrates of complex V to synthesize ATP. Takeuchi et al. [40] conducted a comprehensive review of modeling the tight regulation of ATP levels via calcium, ADP, and inorganic phosphate dynamics.
ROS generation and scavenging regulation is another important topic in mitochondrial OXPHOS modeling. ROS play a crucial role in physiological metabolic signaling and pathological oxidative stress. Pereira et al. [170] and Mazat et al. [171] discussed the mathematical ETC models of ROS production that could be further extended in studies of oxidative stress and signaling. For instance, Kembro et al. [14] built a two-compartment (cytosolic and mitochondrial) model integrating ROS generation from the ETC, ROS scavenging by SOD, catalase, the GSH and the Trx systems, NADPH production, and superoxide-induced mitochondrial depolarization. The authors found that mitochondrial membrane potential oscillation frequencies and amplitudes were dependent on the counteracting forces of ROS production and the antioxidant capacity [172].
Mitochondrial OXPHOS coordinates with other mitochondrial components to regulate ATP production, transportation, and utilization. ATP generated from complex V must be exported from the mitochondrial matrix to the cytosol via ANT, which catalyzes the 1:1 exchange of free ATP to free ADP across the IMM [142,173].
The creatine kinase (CK) system shuttles and buffers ATP synthesis from the OX-PHOS and consumption by biological processes, exchanging high energy phosphates between ATP/ADP and creatine/creatine phosphate couples [174]. Mitochondrial CK (mtCK) couples functionally and structurally with ANT to channel ADP recycling, enhance ATP synthesis, and reduce ROS generation [175][176][177]. Mathematical modeling has shown that mtCK can maintain cellular homeostasis and dampen mitochondrial metabolic oscillation [178].
The ETC is fueled by the CAC. CAC-generated NADH is oxidized in complex I and the CAC intermediate succinate is oxidized in complex II. CAC is coupled with carbohydrate (glycolysis), lipid (beta-oxidation) and amino acid metabolism. Some mod-els keep these inputs as constant parameters [56] for simplicity, while others study the impact of varying fuel inputs. For instance, in a glucose-sensing beta cell model by Fridlyand et al. [179], glycolysis flux can enhance CAC, ETC, and mitochondrial ATP synthesis fluxes and induce insulin secretion. A cardiomyocyte mitochondrial model by Cortassa et al. [165] showed that excessive fatty acids can lead to increased ROS production and mitochondrial uncoupling, both computationally and experimentally.
The ETC architecture in vivo could be more complex than the homogeneous assumption in the ODE-based mathematical models. Respiratory complexes I, III dimer, and IV could assemble into a supercomplex, channeling the electron carriers of ubiquinone and cytochrome c between the complexes to increase efficiency, which may deviate from the pool behavior assumed in ODE models [1,[180][181][182]. The ETC complexes have also been shown to concentrate on the cristae membrane, the infoldings of the inner mitochondrial membrane, with complex V (ATP synthase) dimers sitting on the apex of the cristae. The membrane potential was shown to be higher across the crista IMM than across the noncrista IMM. Each crista can hold different membrane potentials, acting as independent ATP-generating units [125].
To date, several integrative cellular models have been developed to study the interactions of mitochondrial energetics and the cellular environment. Computational models have been used to study the energy supply and demand in cardiomyocytes [37,136,177,183], glucose sensing and insulin secretion in pancreatic beta cells [156,179], metabolic interactions between neurons and astrocytes in the brain [184][185][186], hepatocellular respiration in the liver [187], and metabolic adaptation during exercise in the skeletal muscle [188,189]. Due to tissue differences in energetics demands, mitochondrial distributions and functions are very different in each tissue [190]. To consider mitochondrial energetics, the model parameters must be adjusted and validated by tissue-specific data. The mitochondria and OXPHOS model could also apply to heart failure studies or iPSC studies in patients. For example, induced pluripotent stem cell-derived cardiomyocytes (iPSC-CMs) provide a tool for disease modeling, including mitochondrial cardiomyopathy and metabolic disorders with cardiac phenotypes [191]. Cardiac computational modeling [192,193] with bioenergetics details can provide insight into the mechanisms in disease development.
In this review, we focus on deterministic kinetic models of OXPHOS, which emphasize the dynamic, continuous, and time-dependent interactions of each complex and substrate. Stochastic models of the ETC complex that focus on discrete events of electron transfer within a complex are also available [65,90,194,195]. While these models provide more details for the underlying molecular mechanism than deterministic kinetic models, they require more simulation time than deterministic models. They are more challenging to integrate with other modeling components. Another way to build OXPHOS ATP generation models is constraint-based modeling (CBM), which considers biochemical reactions as a whole and uses a stoichiometric matrix of metabolites and reactions to deduce the best plausible steady-state target flux (e.g., ATP synthesis rate) under some constraints in the reaction fluxes [196]. Although CBM requires far less parameters to fit and could handle genome-scale metabolic model (GEMMs) containing thousands of reactions and hundreds of metabolites, CBM cannot capture transient dynamic responses and complex interactions between components in the same way as kinetic models [197].
In conclusion, we review various mathematical kinetic modeling strategies for OX-PHOS complexes. The OXPHOS is a highly nonlinear system with complicated time-variant responses on a range of inputs, such as respiring substrates, ions (including protons, phosphate, magnesium, calcium, sodium, and potassium), and the mitochondrial membrane potential. Lumped models require the least number of parameters and intermediate state variables but omit some mechanistic details. Mass action kinetic models give the most detailed inner workings but include many parameters that must be fitted and intermediate state variables that must be computed. Quasi-steady state kinetics models for individual complexes sit in between the two. The OXPHOS system is also interconnected to other mitochondrial bioenergetics systems, as represented mathematically by their inputs: depen-dence on respiring substrate concentration and the pmf, and by their outputs: generating ATP, pmf, and ROS. We also briefly discuss other bioenergetics topics not covered in detail, such as the CAC, ATP export, the creatine shuttle, ROS production/scavenging, and tissue/species differences. These OXPHOS modeling strategies and other interconnected systems can help us to build an integrative bioenergetics model for cellular metabolism.

Conflicts of Interest:
The authors declare no conflict of interest.