Mathematical Modeling of ROS Production and Diode-like Behavior in the SDHA/SDHB Subcomplex of Succinate Dehydrogenases in Reverse Quinol-Fumarate Reductase Direction

Succinate dehydrogenase (SDH) plays an important role in reverse electron transfer during hypoxia/anoxia, in particular, in ischemia, when blood supply to an organ is disrupted, and oxygen is not available. It was detected in the voltammetry studies about three decades ago that the SDHA/SDHB subcomplex of SDH can have such a strong nonlinear property as a “tunnel-diode” behavior in reverse quinol-fumarate reductase direction. The molecular and kinetic mechanisms of this phenomenon, that is, a strong drop in the rate of fumarate reduction as the driving force is increased, are still unclear. In order to account for this property of SDH, we developed and analyzed a mechanistic computational model of reverse electron transfer in the SDHA/SDHB subcomplex of SDH. It was shown that a decrease in the rate of succinate release from the active center during fumarate reduction quantitatively explains the experimentally observed tunnel-diode behavior in SDH and threshold values of the electrode potential of about −80 mV. Computational analysis of ROS production in the SDHA/SDHB subcomplex of SDH during reverse electron transfer predicts that the rate of ROS production decreases when the tunnel-diode behavior appears. These results predict a low rate of ROS production by the SDHA/SDHB subcomplex of SDH during ischemia.


Introduction
Succinate dehydrogenase (SDH) is one of the key enzymes of cell energy metabolism, linking the Krebs cycle and the electron transport chain (ETC). SDH reversibly oxidizes succinate to fumarate and transfers the electrons produced by this reaction to the membrane quinone pool, providing ubiquinol QH 2 for oxidative phosphorylation in the cell. SDH can have strong nonlinear properties in both directions, such as hysteresis and bistability in forward succinate-quinone oxidoreductase activity [1] and a "tunnel-diode" behavior in reverse quinol-fumarate reductase direction [2]. Although the "tunnel-diode" behavior, that is, a strong drop in the rate of fumarate reduction as the driving force (over potential) is increased, has been detected in the voltammetry studies about three decades ago [2], however, molecular and kinetic mechanisms of this phenomenon remain unclear and studied in this work.
It is important to note that the "tunnel-diode" behavior in the reverse electron transfer in complex II is observed at the level of the hydrophilic SDHA/SDHB subcomplex. Figure 1 schematically shows the reverse electron transfer in complex II. In this case, the electron donor is quinol QH 2 in the CoQ (quinone-binding) site localized in the hydrophobic SDHC/SDHD subcomplex in the inner mitochondrial membrane. Two electrons are sequentially transferred, one after the other, to the soluble SDHA/SDHB subcomplex The authors who first discovered the tunnel-diode behavior [2] suggested that this phenomenon may be based on a strong decrease in the binding of fumarate or release of succinate when electron transfer is accelerated in the opposite direction, that is, fumarate reduction. We tested this hypothesis by analyzing computer-simulated steady-state rates of reverse electron transfer in the SDHA/SDHB subunits of SDH with a change in the binding/release constants of fumarate/succinate with the active center of SDH during the fumarate reductase reaction and analyzed a dependence of threshold values of the electrode potential when diode-like behavior is observed on different kinetic parameters and concentration of reduced/oxidized redox centers of CII.  Reverse electron transfer in SDH occurs in conditions of oxygen deficiency (hypoxia/anoxia). It was found [3] that complex I (CI) and dihydroorotate dehydrogenase (DHODH) can still deposit electrons into ETC when oxygen is not available as the terminal electron acceptor of the respiratory chain. In this case, cells lacking oxygen reduction accumulate ubiquinol, driving SDH in the reverse direction and depositing electrons onto fumarate, resulting in succinate generation [3] which may be very important during ischemia.
It was proposed that in ischemia when blood supply to an organ is disrupted and oxygen is not available, reversal of mitochondrial complex II (CII) results in ischemic succinate accumulation, that is, ubiquinone Q reduction by CI driving reversal of CII, with fumarate as an electron acceptor yielding succinate oxidation that drives injury at reperfusion [4]. However, the mechanism of experimentally observed ischemic succinate accumulation is controversial [5]. It was shown recently [6] that CII reversal is possible in hypoxic mitochondria but is not the primary succinate source in hypoxic cardiomyocytes or ischemic hearts. That is, ischemic succinate is generated by canonical Krebs cycle activity mostly from α-ketoglutarate and upstream metabolites rather than by mitochondrial CII reversal [6]. A detailed analysis of different models of ischemic succinate accumulation was performed by Chinopoulos in the review [5]. In particular, it was noted that in ischemia, SDH reverses, forming succinate only minorly due to a diode-like property. And besides, mammalian mitochondria lack suitable quinones that could support SDH reversal.
Thus, the main goal of this work was not only to analyze changes in the kinetics of electron transfer responsible for the diode-like behavior of SDH in reverse quinol-fumarate reductase direction but also to study conditions when diode-like behavior is observed.
The authors who first discovered the tunnel-diode behavior [2] suggested that this phenomenon may be based on a strong decrease in the binding of fumarate or release of succinate when electron transfer is accelerated in the opposite direction, that is, fumarate reduction. We tested this hypothesis by analyzing computer-simulated steady-state rates of reverse electron transfer in the SDHA/SDHB subunits of SDH with a change in the binding/release constants of fumarate/succinate with the active center of SDH during the fumarate reductase reaction and analyzed a dependence of threshold values of the electrode potential when diode-like behavior is observed on different kinetic parameters and concentration of reduced/oxidized redox centers of CII.
Besides, we analyzed ROS production in the SDHA/SDHB subcomplex of SDH during the reverse quinol-fumarate reductase activity and predicted that the rate of ROS production decreases at a decrease in the rate of succinate release from the active site of SDH when FADH • and FADH 2 are proposed as the major ROS producing redox centers.
These results predict a low rate of ROS production by the SDHA/SDHB subcomplex of SDH during hypoxia/anoxia, that is, during ischemia.
(D) Simultaneous changes in equilibrium constants K eq13 , K eq14 , and K eq16 due to an increase in the rate constants k −14 and k −16 . The values of the rate constants, k 14 and k 16 , are assumed unchanged and equal to 10 s −1 . Black solid curve (1) corresponds The total rate of reverse electron transfer (V rev_tot ) in SDH is computed as the sum of rates of oxidation of the cluster [2Fe-2S] in reactions 5,6,9,11,13,15 and equal to V rev_tot = V 5 + V 6 + V 9 + V 11 + V 13 + V 15 . The concentration of fumarate and succinate is equal to 1000 and 50 µM, respectively. Figure 2A shows that when a decrease in constants K eq10 and K eq12 occurs due to a decrease in k on rate constants, k 10 and k 12 , a tunnel-diode behavior is observed, but only at very high negative threshold values of the electrode potential E out (about −200 mV), which is significantly higher than the experimentally observed values (about −80 mV [2]). If a decrease in constants K eq10 and K eq12 occurs due to only an increase in k off rate constants, k -10 and k -12 ( Figure 2B), then tunnel-diode behavior is completely absent.
On the contrary, Figure 2C,D shows a decrease in the rate of dissociation of succinate from the active center at reduced FAD forms, that is, a decrease in equilibrium constants K eq14 and K eq16 due to a decrease in the rate constants k 14 and k 16 ( Figure 2C) or an increase in the rate constants k −14 and k −16 ( Figure 2D) upon FAD reduction results in a strong tunnel-diode behavior. Moreover, the threshold potential may be close to the experimentally observed values of about −60 [7] or −80 [2] mV at k 14 = k 16 = 10 −2 s −1 , K eq13 = 0.24 µM −1 , K eq14 = K eq16 = 0.25 µM (curve 3 in Figure 2C).
It is important to note here that a decrease in K eq14 and K eq16 is accompanied by a simultaneous increase in K eq13 according to the "detailed balance" relations [8], i.e., K eq13 = K eq8 · K eq9 /K eq14 . This is surprisingly consistent with the experimentally observed [9] strong increase in the midpoint redox potential of the FAD/FADH· pair when succinate binds to SDH. The midpoint potential of the cluster [2Fe-2S] also increases with succinate binding, but to a much lesser extent [9], so that leads to an increase in the constant equilibrium K eq13 .

E out Dependency of the Rate of Succinate Production in SDH
In order to better understand the mechanism of the tunnel behavior in SDH in the reverse direction, it is convenient to analyze the steady-state E out dependence of the rates V 8 , V 14 , and V 16 at different values of the binding/release constants of fumarate/succinate with the active center of SDH. These dependencies are presented in  case, is about −200 mV and is very far from the experimentally observed values of about −60 [7] or −80 [2] mV.
If the equilibrium constants Keq10 and Keq12 decrease only by increasing the reverse fumarate release constants k-10 and k-12, then nonmonotonicity in total succinate production, i.e., the tunnel effect, is completely absent ( Figure 6).

Figure 3.
Computer simulated dependence of the steady-state rate of succinate production in SDH on the electrode potential Eout at different values of the rate constants of succinate release from the dicarboxylate binding site. The rates of succinate production in reactions 8, 14, and 16 were calculated under simultaneous changes in the values of equilibrium constant Keq13 due to proposed succinate-dependent changes in midpoint potential Em(FAD.suc/FADH • .suc) compared to Em(FAD/FADH • ), as well as equilibrium constants Keq14 and Keq16 of the release of succinate from the active site upon FAD reduction due to a decrease in the rate constants, k14 and k16. (A-D) Black solid curve (1) corresponds to the total rate of succinate production Vsuc_tot = V8 + V14 + V16; blue dashed curve (2)-V8; red dash-dot curve (3)-V16; green dash-dot-dot curve (4)-V14. All the parameter values for results presented in (A-D) are the same as for curves (1-4) in Figure 3C. That is, values of k14, k16, Keq13, Keq14, and Keq16 for the results in (A-D) are the same as for curve (1), (2), (3), (4) in Figure 2C, respectively, and presented in the Figure 2C caption. Computer simulated dependence of the steady-state rate of succinate production in SDH on the electrode potential E out at different values of the rate constants of succinate release from the dicarboxylate binding site. The rates of succinate production in reactions 8, 14, and 16 were calculated under simultaneous changes in the values of equilibrium constant K eq13 due to proposed succinatedependent changes in midpoint potential E m (FAD.suc/FADH • .suc) compared to E m (FAD/FADH • ), as well as equilibrium constants K eq14 and K eq16 of the release of succinate from the active site upon FAD reduction due to a decrease in the rate constants, k 14 and k 16 . (A-D) Black solid curve (1) corresponds to the total rate of succinate production V suc_tot = V 8 + V 14 + V 16 ; blue dashed curve (2)-V 8 ; red dash-dot curve (3)-V 16 ; green dash-dot-dot curve (4)-V 14 . All the parameter values for results presented in (A-D) are the same as for curves (1-4) in Figure 3C. That is, values of k 14 , k 16 , K eq13 , K eq14 , and K eq16 for the results in (A-D) are the same as for curve (1), (2), (3), (4) in Figure 2C, respectively, and presented in the Figure 2C    Computer simulated dependence of the steady-state rate of succinate production in SDH on the electrode potential E out at different values of the rate constants of succinate binding to the dicarboxylate binding site. The rates of succinate production in reactions 8, 14, and 16 were calculated under a simultaneous increase in the values of equilibrium constant K eq13 , as well as a decrease in equilibrium constants K eq14 and K eq16 of the release of succinate from the active site upon FAD reduction due to an increase in the rate constants k −14 and k −16 . (A-D) Black solid curve (1) corresponds to the total rate of succinate production V suc_tot = V 8 + V 14 + V 16 ; blue dashed curve (2)-V 8 ; red dash-dot curve (3)-V 16 ; green dash-dot-dot curve (4)-V 14 . All the parameter values for results presented in (A-D) are the same as for curves (1-4) in Figure 2D. That is, values of k 14 , k 16 , K eq13 , K eq14 , and K eq16 for the results in (A-D) are the same as for curve (1), (2), (3), (4) in Figure 2D, respectively, and presented in the Figure 2D   Computer simulated dependence of the stationary rate of succinate production in SDH on the electrode potential E out at different values of the rate constants of fumarate binding to the dicarboxylate binding site. The rates of succinate production in reactions 8, 14, and 16 were calculated under simultaneous decrease in the values of equilibrium constants K eq10 and K eq12 of binding of fumarate to the active site upon FAD reduction due to a decrease in k on rate constants, k 10 and k 12 . (A-D) Black solid curve (1) corresponds to the total rate of succinate production V suc = V 8 + V 14 + V 16 ; blue dashed curve (2)-V 8 ; red dash-dot curve (3)-V 16 ; green dash-dot-dot curve (4)-V 14 . All the parameter values for results presented in (A-D) are the same as for curves (1, [3][4][5] in Figure 2A. That is, values of k 10 , k 12 , K eq10 , and K eq12 for the results in (A-D) are the same as for curves (1), (3), (4), (5) in Figure 2A, respectively, and presented in the Figure 2A caption. Figure 6. Computer simulated dependence of the stationary rate of succinate production in SDH on the electrode potential Eout at different values of the rate constants of fumarate release from the dicarboxylate binding site. The rates of succinate production in reactions 8, 14, and 16 were calculated under simultaneous decrease in the values of equilibrium constants Keq10 and Keq12 due to an increase in koff rate constants, k-10 and k-12. The values of kon rate constants, k10 and k12, are assumed unchanged and equal to 1 µM −1 s −1 . (A-D) Black solid curve (1) corresponds to the total rate of succinate production Vsuc = V8 + V14 + V16; blue dashed curve (2)-V8; red dash-dot curve (3)-V16; green dash-dot-dot curve (4)-V14. All the parameter values for results presented in (A-D) are the same as for curves (1, [3][4][5] in Figure 2B in the main text. That is, values of k10, k12, Keq10, and Keq12 for the results in (A-D) are the same as for curves (1), (3), (4), (5) in Figure 2B, respectively, and presented in the Figure 2B caption in the main text.

Computational Analysis of the ROS Production Rate in the SDHA/SDHB Subcomplex during Reverse Electron Transfer
Initially, it was analyzed dependence of the concentration of different redox centers of the SDHA/SDHB subcomplex in the reduced state on Eout. This is important for at least two reasons. First, we have to understand how the concentration of ROS-producing redox centers in the reduced state changes with a change in Eout. And secondly, it is necessary to know the interdependence of Eout and the concentration of the cluster [3Fe-4S] in the reduced state in order to extrapolate the Eout dependence of various rates shown in Figures 2-6 on the dependence of the same rates on the concentration of the reduced cluster [3Fe- Figure 6. Computer simulated dependence of the stationary rate of succinate production in SDH on the electrode potential E out at different values of the rate constants of fumarate release from the dicarboxylate binding site. The rates of succinate production in reactions 8, 14, and 16 were calculated under simultaneous decrease in the values of equilibrium constants K eq10 and K eq12 due to an increase in k off rate constants, k -10 and k -12 . The values of k on rate constants, k 10 and k 12 , are assumed unchanged and equal to 1 µM −1 s −1 . (A-D) Black solid curve (1) corresponds to the total rate of succinate production V suc = V 8 + V 14 + V 16 ; blue dashed curve (2)-V 8 ; red dash-dot curve (3)-V 16 ; green dash-dot-dot curve (4)-V 14 . All the parameter values for results presented in (A-D) are the same as for curves (1, [3][4][5] in Figure 2B in the main text. That is, values of k 10 , k 12 , K eq10 , and K eq12 for the results in (A-D) are the same as for curves (1), (3), (4), (5) in Figure 2B, respectively, and presented in the Figure 2B caption in the main text.
It should be noted that the total current in the SDH in the reverse direction is exactly equal to twice the total rate of free succinate production because two electrons are transferred during the formation of one succinate molecule, which is V rev_tot = 2·V suc_tot , where V suc_tot = V 8 + V 14 + V 16 .
The steady-state E out dependencies of the rates V 8 , V 14 , and V 16 at different values of the succinate release constants from the active center of SDH are presented in Figure 3.
First of all, it should be emphasized that the electron transfer rate constants from the cluster [2Fe-2S] to FAD are several orders of magnitude higher than the dissociation constants of succinate. Therefore, electron transfer rates, V 13 and V 15 , significantly prevail over the V 8 and V 14 rates, respectively, at the branching points with succinate dissociation, such as FAD.suc and FADH·.suc. This means that an increase in the rates of reactions 13 and 15 inhibits the succinate dissociation in reactions 8 and 14, respectively. Figure 3A shows that the E out dependence of the succinate dissociation rate V 8 (curve 2) has a nonmonotonic behavior at the basal values of all equilibrium constants. That is, it decreases when the electrode potential E out decreases below a certain threshold value. This drop is due to the strong predominance of the electron transfer rate from the cluster [2Fe-2S] to FAD, V 13 , over the rate V 8 at very low E out values. However, the total succinate production rate V suc_tot (curve 1) increases monotonically with a decrease in E out due to a strong increase in the dissociation rate V 16 (curve 3). The rate V 14 (curve 4) always remains low due to the high rate V 15 .
The 10-fold decrease in succinate dissociation constants k 14 and k 16 from 10 to 1 s −1 ( Figure 3B) results in a strong nonmonotonic (a tunnel-diode behavior) in the E out dependence of the rate V suc_tot (curve 1). This happens because the succinate dissociation rate V 16 becomes very low (curve 3) and cannot compensate for the drop in the rate V 8 (curve 2) when E out decreases below the threshold value. The rate V 14 (curve 4) remains very low. With a further decrease in the constants k 14 and k 16 from 10 to 0.1 ( Figure 3C) and 0.01 s −1 (Figure 3D), the qualitative behavior of the rates V suc_tot , V 8 , V 14 , and V 16 is the same as in the previous case, only a drop in V suc_tot begins at more positive E out values up to −60 mV observed in experiments [7]. Moreover, the total succinate production rate V suc_tot is completely determined only by the rate V 8 (curves 1 and 2).
Thus, the results of computer simulation of a decrease in the rate constants of succinate dissociation k 14 and k 16 and, correspondingly, in K eq14 and K eq16 , taking into account simultaneous increase in K eq13 , can well quantitatively explain the experimentally observed tunnel-diode behaviour in the reverse transfer of electrons to SDH.
Quite different properties of the appearance of nonmonotonicity in the rate V suc_tot are observed when the equilibrium constants K eq14 and K eq16 decrease due to an increase in the constants k −14 and k −16 of the succinate binding to the active center of SDH ( Figure 4). Figure 4A is identical to Figure 3A. First of all, k −14 and k −16 must be increased by 1000 times compared to the basal value of 0.04 µM −1 s −1 , i.e., up to 40 µM −1 s −1 in order for nonmonotonicity to appear in the E out dependence of the V suc_tot ( Figure 4B). Just as in the previous case, the drop in the rate V suc_tot (curve 1) when E out decreases below the threshold value occur due to a decrease in the rate V 16 (curve 3), which is no longer able to compensate for the drop in the rate V 8 (curve 2). There is a further decrease in the rates V 8 (curve 2) and V 16 (curve 3) to almost 0 with E out less than the threshold value when k −14 and k −16 increase up to 400 ( Figure 4C) and 4000 µM −1 s −1 ( Figure 4D). In this case, the total rate of succinate release, V suc_tot , is determined only by the rate V 14 (curves 1 and 4).
As can be seen from Figure 4D (see also curve 4 in Figure 2D), the E out threshold values can approach the experimentally observed values of about −60 mV. However, the values of the equilibrium constant K eq13 satisfying the "detailed balance" relations have to be very high K eq13 = 24 µM −1 at very high rate constants k −14 and k −16 and do not correspond to the experimentally observed values of midpoint potentials E m (FAD.suc/FADH • .suc) and E m ([2Fe-2S]).
As for the change in the succinate generation rates with a decrease in the binding of fumarate to the active center of SDH during FAD reduction, which corresponds to the data on reverse electron transfer in Figure 2A,B, these data are presented in Figures 5 and 6. Computer simulation analysis shows that if a decrease in the equilibrium constants K eq10 and K eq12 occurs only due to a decrease in the fumarate binding constants k 10 and k 12 , then the nonmonotonicity in the E out dependence of V suc_tot is explained as well as in previous cases by a strong drop in the rate V 16 ( Figure 5). However, the threshold potential, in this case, is about −200 mV and is very far from the experimentally observed values of about −60 [7] or −80 [2] mV.
If the equilibrium constants K eq10 and K eq12 decrease only by increasing the reverse fumarate release constants k -10 and k -12 , then nonmonotonicity in total succinate production, i.e., the tunnel effect, is completely absent (Figure 6).

Computational Analysis of the ROS Production Rate in the SDHA/SDHB Subcomplex during Reverse Electron Transfer
Initially, it was analyzed dependence of the concentration of different redox centers of the SDHA/SDHB subcomplex in the reduced state on E out . This is important for at least two reasons. First, we have to understand how the concentration of ROS-producing redox centers in the reduced state changes with a change in E out . And secondly, it is necessary to know the interdependence of E out and the concentration of the cluster [3Fe-4S] in the reduced state in order to extrapolate the E out dependence of various rates shown in Figures 2-6 on the dependence of the same rates on the concentration of the reduced cluster [3Fe-4S]. Then we will be able to predict the possibility of achieving a tunnel effect in real physiological conditions.
The E out dependence of the concentration of FADH • and FADH 2 , as well as the cluster [3Fe-4S] in the reduced state is shown in Figure 7. These redox centers were previously proposed [10][11][12] as the main ROS generators in the subcomplex SDHA/SDHB of SDH. As can be seen from the results presented in Figure 7, only changes in the concentration of FADH • on E out have nonmonotonic dependence ( Figure 7A,B). Moreover, it is interesting that these nonmonotonic changes are only at a simultaneous decrease in the rate constants, k 14 and k 16 , of succinate dissociation from reduced FAD forms ( Figure 7B) can have the experimentally observed values of the threshold potential of about −60 [7] or −80 [2] mV. It is important to note that the concentration of FADH • and FADH 2 decreases with a decrease in the rate of release of succinate from the active center, that is, with a decrease in the equilibrium constants K eq14 and K eq16 ( Figure 7B,D). It should also be noted an important property of the dependence of the concentration of the reduced cluster [3Fe-4S], [3Fe-4S] − , on E out ( Figure 7E). The concentration of [3Fe-4S] − reaches values close to the maximal value of 235 µM already at E out about −60 mV. This means that the "tunneldiode" behavior in reverse quinol-fumarate reductase direction observed experimentally at −60 mV is possible only with maximal [3Fe-4S] cluster reduction, which is apparently difficult to achieve under physiological conditions. At the same time, the dependence of the concentration of the reduced cluster [2Fe-2S] − on E out ( Figure 7F) is much weaker and reaches a maximum value of 235 µM at values of E out about −200 mV. Net reversal of the mammalian SDH complex has been considered thermodynamically unfavorable because the standard reduction potential of ubiquinone (UQ) is slightly greater than that of fumarate [7]. However, because the reduction potential of UQ and fumarate are very close to each other (~10 mV apart), it was shown that UQH 2 accumulation drives the net reversal of the SDH complex in mammalian cells upon suppression of O 2 reduction [3].
The E out dependence of the stationary rates of ROS production is presented in Figure 8. The total stationary rate of H 2 O 2 production, VH 2 O 2 , was computed as the rate of H 2 O 2 release from the mitochondrial matrix to the cytosol, V 22 , that equal to the summary rate of H 2 O 2 production by FADH 2 , V 17 , and dismutation of O 2 · − , V 21 , in the matrix at the steady state (see Explicit functions in Mathematical model in Supplementary Materials). Values of the catalytic constants of ROS formation by different redox centers in the subcomplex SDHA/SDHB were taken from our previous computational model of CII in the forward succinate-quinone oxidoreductase (SQR) direction [13]. That model was calibrated by fitting the computer-simulated results to experimental data obtained on submitochondrial particles prepared from bovine [14] and rat heart [11] mitochondria upon inhibition of the Q-binding site by atpenin A5 and Complex III (CIII) by myxothiazol, respectively. It was shown [13] that only reduced flavin adenine dinucleotide (FADH 2 ) in the unoccupied dicarboxylate state and flavin semiquinone radical (FADH • ) feature the experimentally observed bell-shaped dependence of the rate of ROS production on the succinate concentration upon inhibition of CIII or Q-binding site of CII, i.e., suppression of SQR activity. At the same time, the maximal rate of ROS production was about 1000 pmol/min/mg mitochondrial protein.  and FADH 2 (C) concentration on E out at a simultaneous decrease in k on rate constants, k 10 and k 12 , binding of fumarate to the dicarboxylate binding site upon FAD reduction. In this case, the values of the reverse k off rate constants k -10 and k -12 are assumed unchanged and equal to 50 s −1 . Black solid curve (1) corresponds to k 10 = k 12 = 1 µM −1 s −1 , K eq10 = K eq12 = 2·10 −2 µM −1 ; blue dashed curve (2)-k 10 = k 12 = 10 −1 µM −1 s −1 , K eq10 = K eq12 = 2·10 −3 µM −1 ; red dashdot curve (3)-k 10 = k 12 = 10 −2 µM −1 s −1 , K eq10 = K eq12 = 2·10 −4 µM −1 ; green dash-dot-dot curve (4)-k 10 = k 12 = 10 −3 µM −1 s −1 , K eq10 = K eq12 = 2·10 −5 µM −1 ; pink short-long curve (5)k 10   Simultaneous decrease in equilibrium constants Keq10 and Keq12 of binding of fumarate to the dicarboxylate binding site upon FAD reduction due to a decrease in kon rate constants, k10 and k12. In this case, the values of the reverse koff rate constants k-10 and k-12 are assumed unchanged and equal to 50 s −1 . Parameter values for curves (1-5) are the same as for curves (1-5) in Figure 2A. (B) Simultaneous decrease in equilibrium constants Keq10 and Keq12 upon FAD reduction due to an increase in koff rate constants, k-10 and k-12. In this case, the values of kon rate constants, k10 and k12, are assumed unchanged and equal to 1 µM −1 s −1 . Parameter values for curves (1-4) are the same as for curves (1-4) in Figure 2B. (C) Simultaneous changes in equilibrium constants Keq13, Keq14, and Keq16 describe a decrease in the rate of dissociation of succinate from reduced FAD forms due to a decrease in the rate constants, k14 and k16. The values of the rate constants, k-14 and k-16 are assumed unchanged and equal to 0.04 µM −1 s −1 . Parameter values for curves (1-4) are the same as for curves (1-4) in Figure  2C. (D) Simultaneous changes in equilibrium constants Keq13, Keq14, and Keq16 due to an increase in the rate constants k-14 and k-16. The values of the rate constants, k14 and k16, are assumed unchanged and equal to 10 s −1 . Parameter values for curves (1-5) are the same as for curves (1-5) in Figure 2D. The total rate of H2O2 production, VH2O2, was computed as the rate of H2O2 release from the mitochondrial matrix to the cytosol, V22, that equal to the summary rate of H2O2 production by FADH2, V17, and dismutation of O2• − , V21, in the matrix. The concentration of fumarate and succinate is equal to 1000 and 50 µM, respectively.

Dependence of the Rate of Succinate Production on the Fumarate Concentration in SDH
Computer-simulated dependencies of the stationary rates of succinate production in the reverse direction of SDH on the fumarate concentration at different values of the electrode potential Eout are presented in Figure 9A  (A) Simultaneous decrease in equilibrium constants K eq10 and K eq12 of binding of fumarate to the dicarboxylate binding site upon FAD reduction due to a decrease in k on rate constants, k 10 and k 12 . In this case, the values of the reverse k off rate constants k -10 and k -12 are assumed unchanged and equal to 50 s −1 . Parameter values for curves (1-5) are the same as for curves (1-5) in Figure 2A. (B) Simultaneous decrease in equilibrium constants K eq10 and K eq12 upon FAD reduction due to an increase in k off rate constants, k -10 and k -12 . In this case, the values of k on rate constants, k 10 and k 12 , are assumed unchanged and equal to 1 µM −1 s −1 . Parameter values for curves (1-4) are the same as for curves (1-4) in Figure 2B. (C) Simultaneous changes in equilibrium constants K eq13 , K eq14 , and K eq16 describe a decrease in the rate of dissociation of succinate from reduced FAD forms due to a decrease in the rate constants, k 14 and k 16 . The values of the rate constants, k −14 and k −16 are assumed unchanged and equal to 0.04 µM −1 s −1 . Parameter values for curves (1-4) are the same as for curves (1-4) in Figure 2C. (D) Simultaneous changes in equilibrium constants K eq13 , K eq14 , and K eq16 due to an increase in the rate constants k −14 and k −16 . The values of the rate constants, k 14 and k 16 , are assumed unchanged and equal to 10 s −1 . Parameter values for curves (1-5) are the same as for curves (1-5) in Figure 2D. The total rate of H 2 O 2 production, VH 2 O 2 , was computed as the rate of H 2 O 2 release from the mitochondrial matrix to the cytosol, V 22 , that equal to the summary rate of H 2 O 2 production by FADH 2 , V 17 , and dismutation of O 2 · − , V 21 , in the matrix. The concentration of fumarate and succinate is equal to 1000 and 50 µM, respectively.
As we noted in Supplementary Materials in the section "Mathematical model. Dimension of local and whole mitochondrial concentration and rates." of this work, in order to compare computer-simulated rates of ROS production presented in Figure 8 in µM/s with experimentally observed rates expressed in pmol/min/mg protein the computer-simulated rates should be multiplied by a factor of 220, i.e., 1 µM/s = 220 pmol/min/mg mitochondrial protein. As the simulation results show, only a simultaneous decrease in equilibrium constants K eq10 and K eq12 of binding of fumarate to the dicarboxylate binding site upon FAD reduction due to a decrease in k on rate constants, k 10 and k 12 , ( Figure 8A) can result in very high rates of ROS production more than 1000 pmol/min/mg mitochondrial protein.
That is clear because the decrease in the rate constants k 10 and k 12 results in an increase in the concentration of the major ROS-producing redox centers FADH • and FADH 2 . However, our preliminary analysis shows that a decrease in the rate constants of fumarate binding to the dicarboxylate binding site upon FAD reduction, k 10 and k 12 , is unlikely to account for experimentally observed values of the threshold potential during the "tunnel-diode" behavior and should not be considered as a possible high rate of ROS production at reverse SDH activity. All other changes in the rate constants of binding or release of succinate/fumarate have a small effect on the rate of ROS production by the subcomplex SDHA/SDHB in the reverse direction of SDH ( Figure 8B-D). Thus, a decrease in the rate of succinate release from the active center during the reduction in FAD to FADH 2 that quantitatively explains the experimentally observed tunnel-diode behavior in SDH from the beef heart [2] and Escherichia coli [7] mitochondria results in a very low rate of ROS production in the reverse direction of SDH ( Figure 8C,D) when FADH • and FADH 2 are proposed as the major ROS producing redox centers.
However, it should be pointed out the recent work [12] in which it was shown that in the absence of respiratory chain inhibitors, the model analysis revealed the [3Fe-4S] iron-sulfur cluster as the primary O 2 · − source. In this case, taking into account the very high concentration of [3Fe-4S] − at relatively small values of E out as shown in Figure 8E, it is very likely to expect a high rate of ROS production by this cluster [3Fe-4S] − of SDH during a "tunnel-diode" behavior in reverse quinol-fumarate reductase direction in the absence of respiratory chain inhibitors. Thus, real changes in the rate of ROS production by the subcomplex SDHA/SDHB in the reverse direction of SDH during hypoxia/anoxia, that is, during ischemia, depends on the real experimental condition.

Dependence of the Rate of Succinate Production on the Fumarate Concentration in SDH
Computer-simulated dependencies of the stationary rates of succinate production in the reverse direction of SDH on the fumarate concentration at different values of the electrode potential E out are presented in Figure 9A-D. Figures A-D differ in the values of the equilibrium constants K eq13 , K eq14 , and K eq16 , which change due to a decrease in the rate constants of succinate release from the active center k 14 and k 16 from high control values ( Figure 9A) to very small values ( Figure 9D) when electron transfer is accelerated in the opposite direction. As can be seen from Figure 9, the dependence of the total succinate generation rate on the fumarate concentration has the usual hyperbolic character for any values of E out and equilibrium constants K eq13 , K eq14 , and K eq16 . At the same time, with a decrease in the values of the rate constants k 14 and k 16 , the maximum rate of succinate production begins to decrease when the negative values of E out exceed some threshold values (about −200 mV in Figure 9B), which the less in absolute value (about −150 and −100 mV), the less the rate constants of succinate release from the active center in Figure 9C,D, respectively. generation rate on the fumarate concentration has the usual hyperbolic character for any values of Eout and equilibrium constants Keq13, Keq14, and Keq16. At the same time, with a decrease in the values of the rate constants k14 and k16, the maximum rate of succinate production begins to decrease when the negative values of Eout exceed some threshold values (about −200 mV in Figure 9B), which the less in absolute value (about −150 and −100 mV), the less the rate constants of succinate release from the active center in Figure  9C,D, respectively.  . Computer simulated dependence of the stationary rate of succinate production in SDH on the fumarate concentration at different values of the electrode potential E out and rate constants of succinate release from the dicarboxylate binding site upon FAD reduction. (A-D) Simultaneous changes in equilibrium constants K eq13 , K eq14 , and K eq16 describe a decrease in the rate of succinate dissociation from reduced FAD forms due to a decrease in the rate constants k 14 and k 16 . The values of the rate constants, k −14 and k −16 are assumed unchanged and equal to 0.04 µM −1 s −1 . E out values for black solid curve curves (1) in all figures correspond to 0 mV; blue dashed curves (2) to −50 mV; red dash-dot curves (3) to −100 mV; green dash-dot-dot curves (4) to −150 mV; pink short-long curves (5) to −200 mV. Values of equilibrium constants K eq13 , K eq14 and K eq16 and the rate constants, k 14

Kinetic Model of Reverse Electron Transfer in SDHA/SDHB Subunits of SDH
A kinetic scheme of reverse electron transfer and O 2 · − /H 2 O 2 production underlying a mechanistic computational model of SDH activity in the reverse direction (fumarate reduction) in subunits SDHA and SDHB is presented in Figure 10. This kinetic scheme includes the following redox centers: (a) three iron-sulfur clusters: [2Fe-2S], [4Fe-4S], and [3Fe-4S]) located in SDHB subunit ( Figure 10A) that transfer electrons one at a time to produce succinate in the SDHA subunit; (b) flavin adenine dinucleotide, FAD, located in SDHA subunit ( Figure 10B).
Reverse electron transfer reactions in the SDHA/SDHB subunits of SDH include both the mainstream electron pathway from [3Fe-4S] cluster to fumarate and bypass reactions resulting in O 2 · − /H 2 O 2 formation. These bypass reactions are marked in red in the kinetic scheme ( Figure 10C). This scheme is supported by the numerous literature data on electron transfer pathways between different redox centers of SDH [10,11,15,16].

Kinetic Model of Reverse Electron Transfer in SDHA/SDHB Subunits of SDH
A kinetic scheme of reverse electron transfer and O2• − /H2O2 production underlyin mechanistic computational model of SDH activity in the reverse direction (fumarate duction) in subunits SDHA and SDHB is presented in Figure 10. This kinetic scheme cludes the following redox centers: (a) three iron-sulfur clusters: [2Fe-2S], [4Fe-4S], a [3Fe-4S]) located in SDHB subunit ( Figure 10A) that transfer electrons one at a time produce succinate in the SDHA subunit; (b) flavin adenine dinucleotide, FAD, located SDHA subunit ( Figure 10B).  Reverse electron transfer reactions in the SDHA/SDHB subunits of SDH include both the mainstream electron pathway from [3Fe-4S] cluster to fumarate and bypass reactions resulting in O2• − /H2O2 formation. These bypass reactions are marked in red in the kinetic scheme ( Figure 10C). This scheme is supported by the numerous literature data on electron transfer pathways between different redox centers of SDH [10,11,15,16].
The entire reaction network of electron transfer and O2• −/ H2O2 production corresponding to this kinetic scheme in Figure 10 consists of 22 reactions that are described in detail in Table 1. Midpoint redox potentials, rate constants, and concentrations are taken from the experimental data (see Table 2 and references therein).
The reverse electron transfer in the soluble subcomplex SDHA/SDHB begins with the reduction in the cluster [3Fe-4S] (reaction 1) in the SDHB subunit with subsequent single electron transfer from the [3Fe-4S] cluster through [4Fe-4S] iron-sulfur center (reaction 2) The entire reaction network of electron transfer and O 2 · − /H 2 O 2 production corresponding to this kinetic scheme in Figure 10 consists of 22 reactions that are described in detail in Table 1. Midpoint redox potentials, rate constants, and concentrations are taken from the experimental data (see Table 2 and references therein).
The reverse electron transfer in the soluble subcomplex SDHA/SDHB begins with the reduction in the cluster [3Fe-4S] (reaction 1) in the SDHB subunit with subsequent single electron transfer from the [3Fe-4S] cluster through [4Fe-4S] iron-sulfur center (reaction 2) to the cluster [2Fe-2S] (reaction 3) as shown in Figure 10A. The initial reduction in the cluster [3Fe-4S] (reaction 1) can be carried out by various electron donors, the nature of which depends on experimental conditions. In physiological conditions, electrons in reaction 1 are exchanged with the quinone pool. A detailed kinetic scheme of the exchange of the first and second electrons between the [3Fe-4S] cluster and the quinone pool in CII is presented in work [1]. In voltammetry experiments [2,7], electrons are exchanged directly between SDH and a graphite electrode. The rate of this exchange is controlled by the electrode potential (E out ) and is considered in detail in Tables 1 and 2. Figure 10B presents a kinetic scheme of chemical reactions of the reversible oxidoreduction of succinate, fumarate, and FAD, catalyzed by the SDH flavoprotein subunit A (SDHA). These reactions involve binding/dissociation of fumarate/succinate to/from the dicarboxylate binding site of the SDHA subunit as well as a single electron transfer from the iron-sulfur cluster [2Fe-2S] of the hydrophilic SDHB subunit to FAD in the SDHA subunit. Reactions of the reactive oxygen species formation in the soluble subcomplex SDHA/SDHB of SDH are presented in Figure 10C.

Release of hydrogen peroxide (H 2 O 2 ) from the mitochondrial matrix to cytosol
22  The reference for the equilibrium constant K eq . c The reference for the rate constant of direct reaction k forward . d The used value of K eq7 = 2778 is calculated from the relation K eq7 · K eq8 · K eq12 = exp (2 · F · (E (fum/suc) − E (FAD/FADH2))/R · T) = 555.6, where midpoint redox potentials E (FAD/FADH2) = −79 mV (pH 7.0) and E (fum/suc) = 0 mV (pH 7.0) [22], respectively, and F, R, and T have the usual meaning. In the general case, when the values of the K eq8 and K eq12 can vary, K eq7 = 555.6/K eq8 · K eq12 . e Relations between equilibrium constants according to 4 thermodynamic cycles in the kinetic scheme: (1) K eq5 = 1/K eq4 · K eq6 ·K eq7 · K eq8 ; (2) K eq6 = K eq11 · K eq12 /K eq10 ; (3) K eq13 = K eq8 · K eq9 /K eq14 ; (4) K eq15 = 1/K eq7 · K eq12 · K eq13 · K eq16 ; f The used value was taken from [21], which was calculated from experimental data on Mn-SOD activity in mitochondria of cardiac cells [23]. Conserved moieties (in µM). The pool of electron carriers. In this work, we assumed that the concentration of the SDHA/SDHB subunits of SDH is equal to the concentration of Complex II. According to [24], the content of Complex II in cardiac mitochondria is 0.209 nmol Complex II/mg of mitochondrial protein. Translation of whole membrane concentration expressed in nmol/mg mit.prot. to local protein concentration expressed in µM presented in [21]. We have shown earlier [21] that 1 nmol/mg of protein corresponds to 273 µM when normalized to the mitochondrial volume (V mit ). If the concentration is normalized to the inner mitochondrial membrane volume (V imb ) it should be additionally taken into account that the ratio W imb = V imb /V mit = 0.24 [21]. Therefore 0. The reduction in fumarate to succinate in the flavoprotein subunit SDHA shown in Figure 10B, can occur via various alternative pathways. Fumarate can first bind to the dicarboxylate binding site of SDH when FAD is oxidized (reaction 4), followed by the transfer of the first (reaction 5) and second (reaction 6) electrons from the cluster [2Fe-2S] to FAD with the formation of the complexes FADH • .fum and FADH2.fum, respectively. In reaction 7, fumarate is reduced to succinate with simultaneous oxidation of FADH 2 to FAD. Then, in reaction 8, the FAD. suc complex dissociates with the release of succinate.
In another pathway of reducing fumarate to succinate, FAD in the unoccupied dicarboxylate state is sequentially reduced first to FADH • in reaction 9 and then to FADH 2 in reaction 11 receiving the first and second electrons, respectively, from the cluster [2Fe-2S]. In this case, fumarate can bind to the active site when FAD is in FADH • (reaction 10) or FADH 2 (reaction 12) state, respectively. After that, fumarate reduction to succinate and succinate release occur as well as in the previous pathway that is in reactions 6-8 or 7, 8.
In addition, reactions 13-16 represent the third pathway of two electrons transfer from the cluster [2Fe-2S] to the FAD and the release of succinate when succinate is initially bound to the active center of SDH.
These three pathways of formation of succinate from fumarate presented in Figure 2B include 4 thermodynamic cycles in which the initial and final states are identical (reactions 4-8; 6, 10-12; 8, 9, 13, 14; 7, 12, 13, 15, 16). Therefore, the equilibrium constants of the reactions along any cycle must satisfy so-called "detailed balance" relationships [8]. These detailed balance relations require the product of the equilibrium constants along a cycle to be equal to 1, as at equilibrium, the net flux through any cycle vanishes. Therefore, such relations decrease the number of independent rate constants in a kinetic model and imply that when any one of the equilibrium constants in any cycle changes, other constants in this cycle should automatically change.

Computational Model of Reverse Electron Transfer in SDHA/SDHB Subunits of SDH. Mathematical Model
A computational model consisting of 13 ordinary differential equations (ODEs) for the kinetic scheme presented in Figure 2, plus 4 moiety conservation equations, was derived from the reaction networks using the law of mass action and Michaelis (more exactly: Henri-Michaelis-Menten [25]) equation for all 22 kinetic processes. The models were implemented in DBSolve Optimum software available at the website http://insysbio.ru, accessed on 6 December 2022. The details of the mathematical model describing oxidized and reduced states of different carriers and electron flow through SDHA/SDHB subunits of SDH are presented in Supplementary Materials. Values of model parameters, rate constants, and concentration of different electron carriers were taken from experimental data in the literature on thermodynamics and kinetics of electron transfer in CII of the respiratory chain (Tables 1 and 2). Additionally, the model is presented in SBML format in a separate file: Markevich_Final Reverse Scheme SDHA-B.xml (Supporting information).

Conclusions
A mechanistic computational model of reverse electron transfer in the subcomplex SDHA/SDHB of SDH was developed to account for a "tunnel-diode" behavior in reverse quinol-fumarate reductase direction in the SDHA/SDHB subcomplex detected in the voltammetry studies about three decades ago [2]. The model consists of a system of 13 ordinary differential equations that describe the dependence of the rates of electron flows and succinate production in SDH on the electrode potential E out . It was found that the "tunnel-diode" behavior in the SDHA/SDHB subcomplex of SDH, that is, a strong drop in the rate of fumarate reduction as the driving force is increased, is accounted for by a decrease in succinate release from the active site of SDH during reduction in flavin adenine dinucleotide, FAD, to FADH 2 . In particular, the model simulation predicts that experimentally observed threshold values of the electrode potential of about −60 [7] or −80 [2] mV can be explained quantitatively by a decrease in the rate constants of succinate release, k 14 , and k 16 , from the control value of 10 s −1 (curve 1 in Figure 2C) to k 14 = k 16 = 10 −2 s −1 (curve 4 in Figure 2C).
It is important to emphasize that a decrease in dissociation constants of succinate from the active site is accompanied by a simultaneous shift of equilibrium from FAD to FADH· according to the "detailed balance" relations [8] in thermodynamic cycles of succinate formation and release as we noted in Section 3.1. "E out dependency of the total rate of reverse electron transfer in SDH." This is surprisingly consistent with the experimentally observed [9] strong increase in the midpoint redox potential of the FAD/FADH· pair when succinate binds to SDH.
Computational analysis of ROS production in the SDHA/SDHB subcomplex of SDH during the reverse quinol-fumarate reductase activity predicts that the rate of ROS production decreases at a decrease in the rate of succinate release from the active site of SDH when FADH • and FADH 2 are proposed as the major ROS producing redox centers. The concentration of FADH • and FADH 2 decreases with a decrease in the equilibrium constants K eq14 and K eq16 as shown in Figure 7B,D. However, the concentration of the reduced cluster [3Fe-4S] − of SDH is maximal during a "tunnel-diode" behavior in reverse quinol-fumarate reductase direction, so it is very likely to expect a high rate of ROS production by this cluster if [3Fe-4S] is the primary O 2 · − source redox center in the absence of respiratory chain inhibitors as proposed in [12].
As to the participation of reverse electron transfer in SDH in ischemic succinate accumulation, it seems that it depends on real physiological conditions since SDH can be a major succinate source in ischemic in some cases [4]. However, sometimes ischemic succinate is generated by canonical Krebs cycle activity mostly from α-ketoglutarate and upstream metabolites rather than by mitochondrial CII reversal [6]. Our computational analysis predicts that a "tunnel-diode" behavior in the reverse quinol-fumarate reductase direction is observed only at the maximal reduction in the [3Fe-4S] cluster. So, it is unlikely that a tunnel effect can block ischemic succinate accumulation, as proposed in the review [5]. It seems more likely, that complex I and dihydroorotate dehydrogenase can support reverse electron transfer in SDH when oxygen is not available as the terminal electron acceptor of the respiratory chain, as shown in [3].