📒 Heiske 2017

Comprehensive mathematical model of oxidative phosphorylation valid for physiological and pathological conditions1



  • Various models have hence been established to address different problems with respect to the mitochondrial energy metabolism.
    • Magnus and Keizer have used two complementary rate equations to describe the respiratory chain activity, one for calculating the oxygen consumption and the other one for assessing the flux of the translocated protons
    • All equations included a dependency on the membrane potential ΔΨ, and Cortassa et al. modified these equations in order to take into account also ΔμH
    • Yugi and Tomita: impact of ΔμH was not considered
    • Korzeniewski and Zoladz and Beard: rate equation based on equilibrium thermodynamics or on mass action law. Do not exhibit saturation kinetics
  • develop an OXPHOS model, where each OXPHOS complex is described by a separate rate equation, which reproduces the enzyme kinetics over a wide range of substrate and products concentrations, as well as ΔμH
  • Michaelis–Menten type rate equation based on a sequential random binding mechanism, respecting thermodynamics



  • derived from the previous models of Beard and Korzeniewski and Zoladz
  • modified the rate equations for some reactions in order to take into account classical kinetic parameters and to respect thermodynamics

Rate equations for the OXPHOS complexes and influence of the electrochemical potential difference

  • the kinetics of complex I in the absence of the electrochemical potential difference ΔμH can be well described by a Henri–Michaelis–Menten type rate equation based on random binding

  • general reaction scheme:

  • In the absence of ΔμH:

  • Haldane relationship:

  • equilibrium constant is related to the Gibbs energy of reaction

  • Under the influence of the electrochemical potential difference

  • The resulting Haldane equation

  • modulate the maximum rates k and the Km by an exponential expression

  • distribution of ΔGH

  • Satisfying the thermodynamic relationship

  • General forms of reaction rates

Rate equations of ETC complexes

Rate equation for the ADP/ATP carrier (AAC, ANT)

Rate equation for the proton leak

  • Korzeniewski and Zoladz model (exponential relationship)

Rate equation for the Pi/OH exchange

Rate equation for the adenylate kinase reaction

Parameters of the OXPHOS complexes

Model parameters

Mitochondrial geometry and contents

Kinetic parameters of the respiratory chain (class A)

  • initial rate measurements on complex III and complex IV on bovine heart mitochondria
Complex III

Complex IV

Parameters values taken from literature data (class B)

  • kinetic constants for the ATP synthase and the AAC, dissociation constants of mADP, mATP, and phosphate

Empirical parameters (class C)

Model validation

  • model can reproduce the Bose dataset, and is able to describe the steady‐state fluxes through the OXPHOS complexes under phosphorylating and nonphosphorylating conditions

Mitochondrial ATP production

  • how tight the mitochondrial ATP production is coupled to the respiratory chain activity. This efficiency can be characterized by the P/O ratio
  • The P/O ratio obtained at steady state was 2.64 which is in the range of experimental determined values 47 and which is very close to 2.67, the theoretical maximum
  • only an exponential‐shaped proton leak flux in function of ΔΨ can correctly reproduce a near zero leak flux at state 3 as well as the given leak flux at state 4

Flux–force relationship

  • uncouplers
  • respiration substrate
  • ADP (from 0 to 1.3 mM)
  • ATP synthase inhibitor
  • inhibition of the ETC (complex I)
  • in good agreement with the experimental flux–force curves of Amo and Brand 49 (Fig. 6A) and those of Hafner et al. 50 (Fig. 6B)

Simulation of time courses
Extension of the differential equation system of the OXPHOS model

Flux control coefficients

  • Flux control coefficients (FCCs):

approximated by

  • The FCCs of C1, C3, and C4, the ATP synthase, and the AAC are in good agreement with the experimental data. The FCC for the PiOH is lower in our model (due to lower Km)

Application to threshold curves
Antimycin A effects
Antimycin A effects
oligomycin effects
oligomycin effects
threshold curves at state 3
threshold curves at state 4

  • This corresponds to the tendencies observed experimentally for complex I and III in rat liver mitochondria 60 and for complex IV in human hepatic cancer cell line (HepC2)


New OXPHOS rate equations

  • generic OXPHOS rate equation based on a Henri–Michaelis–Menten‐like rate equation, with influence of ΔμH as well as the substrate and product affinities (Km)
  • new generic rate equation could be adapted to the number of substrate and products and to the stoichiometric factors of each respiratory chain complex
  • For complex III, we did not take into account both ubiquinone sites (ubiquinone_o and ubiquinone_i). Using one apparent Km for ubiquinone and one for ubiquinol was thus more convenient and lead to satisfactory results for the description of the initial rate measurements
  • these new rate equations are able to reproduce the steady‐state kinetics experimentally found for respiratory chain complexes I, III, and IV in the absence of ΔμH
  • The approach respects thermodynamics and allows us to reproduce realistic flux–force relationships

Model parameterization

  • a large part of the parameters (class A and B) correspond to classical kinetic parameters which can be determined experimentally
  • used the dataset of Bose in order to estimate the class C parameters
  • the number of ζ and γ (free energy distribution in polarized mitochondria) to adjust can be importantly reduced (six parameters less to adjust), as the Km values of ubiquinone and ubiquinol of both, complexes I and III, and those of ADP and ATP of the ATP synthase were found to be insensitive to the membrane potential

Activation of complex III by phosphate

  • activation of the respiratory chain upon addition of Pi (experimental data of Bose)
  • for simplicity we modeled only complex III to be influenced by Pi, similar as it has been done by Beard

Model validation

  • the model was able to describe well the large dataset of Bose et al. 17, comprising respiration rates, ΔΨ, matricial proton concentrations, and the redox states of NADH and cytochrome c
  • correct P/O ratio of 2.64
  • the simulated time courses of ΔΨ, O2, and the redox state of NADH showed a realistic behavior and were in good agreement with experimental data
  • successfully confronted our model with different experimental datasets

Application to threshold curves and future applications

  • simulate threshold curves and predicted thereby the influence of an inhibition of a given OXPHOS system component on the respiratory flux
  • For state 4 respiration, the simulated thresholds were higher than for state 3
  • These good agreements could not be obtained when using OXPHOS rate equation based on mass action as described by Beard

Materials and methods

Conversion of catalytic activities (nmol·min−1·mg mito. prot.−1) into model rate constants (mol·s−1·(L mito. water space)−1

  • 1 g mitochondrial protein corresponds to 2.28 mL mitochondrial water diffusion space

f_prot = 7.3E-6

Conversion of respiration rates in nmol O2·min−1·nmol cytochrome a−1 into model rate constants (mol·(L mito. water space)−1·s−1)

  • fcyta: the cytochrome a content in mammalian heart mitochondria. Its value was estimated within the above given range together with other adjustable model parameters

Comparison of Km* values with literature Km values

  • we make the approximations that k is the same for both kinetic models, and that for [S] = Km, we have v = k/2


  1. Heiske M, Letellier T, Klipp E. Comprehensive mathematical model of oxidative phosphorylation valid for physiological and pathological conditions. FEBS J. 2017;284(17):2802-2828. doi:10.1111/febs.14151. FEBS↩︎