Initially, in order to calibrate the developed model, computer simulated dependencies of the rates of ROS production on the AA5, inhibitor of the Q-binding site, and succinate concentration were fitted to experimental data on SMP prepared from bovine heart mitochondria (Siebels and Dröse, 2013). Computer simulated results and experimental data presented in Figure 2 show a good fit for both dependencies of the H₂O₂ production rate on the AA5 concentration at the fixed succinate concentration (100 μM; Figure 2A) and on the varied succinate concentration at the fixed AA5 concentration (0.05 μM; Figure 2B). Fitting resulted in changes in values of some model parameters compare to the initial values including literature data. The new adjustable values of some parameters are presented in Supplementary Table S3. The most significant deviations of the adjusted values from the initial values are observed for the catalytic constants of the rate of ROS generation by each redox center and the catalytic constants of Q binding to the Q-binding center. Computer simulated time course of the total H₂O₂ production rate at changes in the AA5 concentration from 0 up to 0.05, 0.1, and 0.15 μM at fixed succinate concentration (100 μM; Figure 2C) confirms stationary fit computer simulated and experimental data presented in Figure 2A. Computer simulated dependencies of the stationary rate of QH2 production (succinate oxidation) on the succinate concentration at different fixed AA5 concentrations presented in Figure 2D show that the maximum reaction rate and the Michaelis constant decrease proportionally with increasing AA5 concentration which is typical for a competitive inhibitor.

Computer simulated dependence of stationary rates of ROS production by different redox centers of CII, namely: the flavin site, [3Fe-4S] cluster, and semiquinone, CII.Q⁻, at the Q-binding site on the succinate concentration with varying degrees of inhibition of the Q-binding site by AA5 are presented in Figure 3. These computer simulation results predict that only reduced FADH₂ in the unoccupied dicarboxylate state (Figures 3A,B) and FADH^• (Figure 3C) have the experimentally observed [5] bell-shaped dependence of the rate of ROS production on the succinate concentration upon inhibition of the Q-binding site by AA5. Succinate-dependence of the rate of ROS formation at these sites (Figures 3A–C) as well as the total rate of ROS generation by CII (Figure 3F) changes from a very small amplitude in the basal state (AA5 = 0) to the intermediate‐ and high-amplitude bell-shaped kinetics with a shift to small succinate concentration upon an intermediate (AA5 = 0.05 μM) and strong (AA5 = 0.15 μM) inhibition of the Q-binding site, respectively.

It is proposed that FADH₂ can generate both hydrogen peroxide, H₂O₂, and superoxide, O₂⁻, with the rates v₂₂ and v₂₃, respectively (Figure 1A; Supplementary Table S1). The ratio of the catalytic constants of H₂O₂, and O₂⁻ formation k₂₂ and k₂₃ are unknown, one authors (Siebels and Dröse, 2013) found 75% H₂O₂ and 25% O₂⁻ in rat heart mitochondria while another (Grivennikova et al., 2017) found that bovine heart mitochondrial respiratory CII generates ROS, mostly as superoxide. Our fitting results predict close values for these constants: k₂₂ = 0.027 μM⁻¹⋅s⁻¹ and k₂₃ = 0.019 μM⁻¹⋅s⁻¹, so v₂₂ (Figure 3A) and v₂₃ (Figure 3B) look similar although we have to take into account that v₂₃ is the rate of O₂⁻ production, so the contribution of v₂₃ in the total rate of H₂O₂ production two times less than v₂₂ because two molecules O₂⁻ give one molecule H₂O₂ in the process of subsequent dismutation (reaction 28) in Supplementary Table S1. The stationary total rate of H₂O₂ production by CII (vH₂O_2tot) was computed as the rate of H₂O₂ release from the mitochondrial matrix to cytosol that equal to the summary rate of H₂O₂ production by FADH₂ in assembled, v₂₂, and disintegrated, v_22d, states and dismutation of O₂⁻, v₂₈, in the matrix at the steady state (see Explicit functions in Mathematical model in Supplementary Data).

The computer simulated dependence of stationary rates of ROS (O₂⁻) production by [3Fe-4S] cluster, v₂₅, (Figure 3D) and semiquinone at the Q-binding site, v₂₆, (Figure 3E) on the succinate concentration show monotonic hyperbolic kinetics in the basal state (AA5 = 0) as well as upon an intermediate (AA5 = 0.05 μM) and strong (AA5 = 0.25 μM) inhibition of the Q-binding site by AA5. It should be pointed that AA5 increases the sensitivity of O₂⁻ production by [3Fe-4S] cluster, v₂₅, to succinate, i.e., decreases Michaelis constant, K_m25, and does not affect on the maximal rate v₂₅. On the contrary, AA5 decreases proportionally the maximal rate of O₂⁻ production by semiquinone at the Q-binding site, v₂₆, and K_m26.

Just like in the work (Markevich et al., 2019), inhibition of CIII was simulated by decreasing the catalytic constant k₂₉ of QH₂ oxidation in the mitochondrial inner membrane [reaction (29) in Supplementary Tables S1, S2]. Figure 4 shows that model simulation of the stationary rate of H₂O₂ production at k₂₉ = 0.005 s⁻¹ (k₂₉ = 1 s⁻¹in the uninhibited state) good fit to experimental data on SMP from rat heart mitochondria upon inhibition of Complex III by myxothiazol (1.6 ìM; Grivennikova et al., 2017). It should be emphasized that values of all parameters including adjustable are the same for both modeling simulation presented in Figures 2A,B, 4 except the total concentration of CII, CIIt, that is equal to 235 and 97 μM for Figures 2, 4, respectively. This fact strongly suggests that the calibration of the developed mathematical model is correct.

Computer simulation modeling results upon inhibition of CIII confirm suggested earlier hypothesis that the non-monotonic succinate-dependence of the rate of ROS generation in CII results from ROS formation at the flavin in the unoccupied dicarboxylate binding site (Quinlan et al., 2012). Computational modeling analysis shows that only FADH₂ and FADH^• have the experimentally observed bell-shaped dependence of the rate of ROS production on the succinate concentration upon inhibition of CIII and predict that inhibition of CIII has the same effects on the kinetics of ROS production by different CII redox centers (Supplementary Figure S1) as inhibition of the Q-binding site considered above (Figure 3). Such similar effect of inhibition of CIII and Q-binding site on the total ROS production by CII was pointed earlier as a result from suppression of SQR activity of CII (Siebels and Dröse, 2013). The stronger inhibition of CIII (a more decrease in k₂₉) results in the more amplitude of the bell-shaped succinate-dependence of the rate of ROS production by the flavin site (Supplementary Figures S1A–C), and the total rate of ROS production by CII (Supplementary Figure S1F) and more shift of the maximal rate of ROS production to the small succinate concentration. The dependence of the stationary rate of ROS formation by [3Fe-4S] cluster, v₂₅, (Supplementary Figure S1D) and semiquinone at the Q-binding site, v₂₆, (Supplementary Figure S1E) on the succinate concentration shows hyperbolic kinetics in the basal state (k₂₉ = 1 s⁻¹) as well as upon an intermediate (k₂₉ = 0.1 s⁻¹) and strong (k₂₉ = 0.01 s⁻¹) inhibition of CIII with very small amplitude. Thus, these computer simulation data predict that the effect of inhibition of CIII and the Q-binding site on ROS production by different sites of CII is very similar.

The hypothesis that a decrease in the rate of ROS production by flavin in the unoccupied dicarboxylate binding site at the high succinate concentration may be a reason for bell-shaped dependence of the total rate of ROS production by CII upon inhibition of CIII was first proposed in (Quinlan et al., 2012). It is really easy to understand a decrease in the concentration of ROS-producing site unoccupied FADH₂ due to succinate binding to FADH₂ in the reaction (5; Figure 1A; Supplementary Table S1). More difficult to understand hyperbolic dependence of the rate of ROS formation on the succinate concentration by [3Fe-4S] cluster and semiquinone at the Q-binding site as well as the concentration of these sites and other Fe-S clusters (Supplementary Figure S2) located downstream flavin. Moreover, experimentally observed catalytic activity of CII has also hyperbolic dependence on the succinate concentration (Quinlan et al., 2012), although at the first glance, bypass ROS production rate by these sites as well as the mainstream electron flow in CII, i.e., the rate of QH₂ production, v₂₀, should have non-monotonic dependence on the succinate concentration because of bell-shaped dependence of the concentration of FADH₂ and FADH^• (Supplementary Figure S2) that are substrate and product, respectively, in the mainstream and bypass reactions (6) and (23) (Figure 1A).

Computer simulation results presented in Figure 5 explain this situation. First, the computer simulated stationary rate v₂₀ really has the hyperbolic dependence on the succinate concentration (Figure 5A) and decreases at decreasing k₂₉, i.e., decreasing QH₂ oxidation that results in a suppression of SQR activity. Second, Figures 5B,C show why the rates v₆ (that equal v₂₀ in the steady state) and v₂₃ (Figure 1A; Supplementary Table S1) that have one the same substrate and one the same product have qualitatively different kinetics, hyperbolic (v₆), and bell-shaped (v₂₃). Figure 5B shows that unidirectional mainstream electron flows in forward (from FADH₂ to oxidized [2Fe-2S] cluster), v_{6_forward}, and reverse (from reduced [2Fe-2S]⁻ cluster to FADH^•), v_{6_reverse}, direction (expression for them in the section Supplementary Data Mathematical model) really follow nonmonotonic FADH₂ and FADH^• concentration and have bell-shaped kinetics. However, the netto rate: v_{6_netto} = v_{6_forward −} v_{6_reverse} has hyperbolic kinetics due to common-mode inphase changes with a close amplitude in v_{6_forward} and v_{6_reverse} (Figure 5B). Values of these rates v_{6_forward} and v_{6_reverse} are so high compare to v_{6_netto} and close each other that curve for them are indistinguishable in Figure 3B. On the contrary, v_{23_reverse} is very small compare to v_{23_forward} (Figure 5C) due to a small concentration of the one product of this reaction, superoxide O₂⁻. That is why the netto rate of superoxide O₂⁻ production by FADH₂, v_{23_netto} = v_{23_forward −} v_{23_reverse} follows the forward rate v_{23_forward} only and has bell-shaped kinetics.

It is worth noting that the succinate concentration required for optimal (peak) ROS generation by CII upon inhibition of CIII or Q-binding site is different in different experiments and changes from approximately 50–500 μM in the experiments with both SMP and intact mitochondria (Quinlan et al., 2012; Siebels and Dröse, 2013; Grivennikova et al., 2017). It was proposed earlier (Grivennikova et al., 2017) that a high optimal succinate concentration up to 400–500 μM in intact mitochondria may be related to limitations of succinate permeability into the mitochondrial matrix, where succinate dehydrogenase active site is located. These limitations were simulated in the model by changes in the rate constant of succinate binding to FAD, k₁, because diffusion is one of the steps of succinate binding. Computational modeling results on changes in the rate constant k₁ presented in Supplementary Figure S3 support, in part, this hypothesis. A decrease in k₁ results in a shift of the optimal succinate concentration to the high values of succinate with simultaneous decrease in the maximal rate of ROS production. These changes resemble the effect of inhibition of CIII and Q-binding site. Therefore, it seems most likely that different succinate optimal values of ROS production observed in different experiments are related to both different power of inhibition of CIII or Q-binding site in different experiments and transport limitations of succinate.

It should point out that our preliminary computational modeling analysis of the simplified model of assembled CII without heme b as an electron carrier in the electron transfer pathway from succinate to Q at the Q-binding site (Markevich et al., 2019) predicts similar effects of inhibition of CIII on ROS production by CII, that is the high-amplitude bell-shaped dependence of the rate of ROS production by CII on the succinate concentration.

Computational analysis of the present extended model with heme b, i.e., included the thermodynamic cycle [3Fe-4S]↔heme b↔Q↔[3Fe-4S] [reactions (9, 11, and 13) for the first electron and reactions (17, 18, and 21) for the second electron transfer] that was experimentally studied in detail in (Anderson et al., 2014) shows that the value of midpoint potential of the heme b, E_m(b), does not affect on the rates of QH₂ and ROS production by CII in the steady state, although affects on the rates of electron transfer in the thermodynamic cycle including heme b (Figure 6). The effect of changes in E_m(b) values from −185 mV used in the present model for bovine heart CII (Supplementary Table S2) up to +36 mV for Escherichia coli (Anderson et al., 2014) on the stationary rates of electron transfer in CII and reduction of heme b was studied computationally.

Figure 6A shows that the stationary concentration of reduced heme b, b⁻, increases slowly up to the small concentration, 12 μM, that is about 5% of the total heme b concentration that equal to 235 μM, with an increase in the succinate concentration at E_m(b) = −185 mV while increasing E_m(b) up to +36 mV results in almost complete reduction of heme b, an increase in b⁻ concentration up to 230 μM, or 98% of the total heme b concentration, at 40 μM succinate concentration. However, despite a strong difference in the reduced steady-state heme b concentration at different E_m(b) values the rates of QH₂ production, v₂₀, i.e., the SQR activity, ROS production, and vH₂O_2tot, remain unchanged (Figure 6B).

As to the rates of electron transfer related to the thermodynamic cycle [3Fe-4S]↔heme b↔Q↔[3Fe-4S], that is reactions (9, 11, and 13) for the first electron and reactions (17, 18, and 21) for the second electron transfer, the computer simulation results presented in Figures 6C,D for E_m(b) = −185 mV and Figures 6E,F for E_m(b) = +36 mV show that these rates change very much when changing E_m(b). First of all, it should be pointed that input electron flow entering the thermodynamic cycle from [4Fe-4S] cluster, v₈ for the first electron and v₁₆ for the second electron transfer, remain unchanged at different E_m(b) values and equal to the output flow, i.e., the rate of QH₂ production, v₂₀. Other curves in Figures 6C–F describe the rates of electron flows inside the thermodynamic cycle. The most interesting of them is v₁₁ at E_m(b) = +36 mV (Figure 6E) that describes the first electron transfer between heme b and Q. Computer simulation results predict that the first electron transfer occurs in the reverse direction from semiquinone to heme b at the high value of E_m(b) equal to +36 mV. However, this effect is compensated by increasing the first electron flow in forward direction from [3Fe-4S] cluster to Q, v₁₃.

It is necessary to emphasize that it was taken into account “the principle of detail balancing” (Hearon, 1953) for thermodynamic cycles at all computations that requires the product of equilibrium constants along a cycle to be equal to 1. For the thermodynamic cycle [3Fe-4S]↔heme b↔Q↔[3Fe-4S] this means the following relations: K_eq13 = K_eq9·K_eq11 and K_eq21 = K_eq17·K_eq18. The values of K_eq9, K_eq11, and K_eq17, K_eq18 at E_m(b) = −185 mV are presented in Supplementary Table S2. Increasing E_m(b) from −185 up to +36 mV results in the following changes in these equilibrium constants: K_eq9 = K_eq17 = 5.55∙10⁻⁵ exp. [(185 + 36) F/RT] = 5.55∙10⁻⁵∙6,905 = 0.38; K_eq11 = 2.72/6905 = 3.94∙10⁻⁴; and K_eq18 = 3.269∙10⁶/6905 = 473.4.

That is, the values of equilibrium constants K_eq13 and K_eq21 kept unchanged while K_eq9, K_eq11, and K_eq17, K_eq18 are varied with increasing E_m(b) value.

These results imply that a high level of ROS production by CII induced by mutations in SDH cyt b (Ishii et al., 1998) is not related with changes in E_m(b). It is most likely as was pointed earlier (Tran et al., 2007; Anderson et al., 2014) that heme b plays more the structural role stabilizing CII as a heterotetramer than for catalysis in CII as an electron carrier. That is, the heme b plays an important role in assembly CII because as was shown in (Yu et al., 1987), in vivo, SDH is anchored to the inner membrane with the cytochrome b₅₆₀. Thus, mutations in the cytochrome b large subunit (SDHC), of CII, that induce oxidative stress and lead to apoptosis (see for review Ishii et al., 2007) results in a suppression of the SQR activity of assembled CII due to its disintegration and it is very likely that in this case the dependence of stationary rates of ROS production on the succinate concentration has the same features as upon inhibition of CIII or Q-binding site considered above.