Complex I is bypassed during high intensity exercise

Human muscles are tailored towards ATP synthesis. When exercising at high work rates muscles convert glucose to lactate, which is less nutrient efficient than respiration. There is hence a trade-off between endurance and power. Metabolic models have been developed to study how limited catalytic capacity of enzymes affects ATP synthesis. Here we integrate an enzyme-constrained metabolic model with proteomics data from muscle fibers. We find that ATP synthesis is constrained by several enzymes. A metabolic bypass of mitochondrial complex I is found to increase the ATP synthesis rate per gram of protein compared to full respiration. To test if this metabolic mode occurs in vivo, we conduct a high resolved incremental exercise tests for five subjects. Their gas exchange at different work rates is accurately reproduced by a whole-body metabolic model incorporating complex I bypass. The study therefore shows how proteome allocation influences metabolism during high intensity exercise.

Substrate efficiency with and without complex I bypass enzyme from yeast (NDI). A) Addition of the yeast NDI reaction to the model increased the attainable area in pareto plot of catalytic capacity vs substrate efficiency. The NDI enzyme was assumed to have a specific activity of 500 mmol mg -1 min -1 . B) The predicted flux through the GLY-PHOS shuttle was abolished in the presence of the NDI enzyme. C) the NDI reaction carried similar flux as GLY-PHOS reaction did in the absence of NDI.

Figure 3 Fitted slopes below and above the inflection point for 5 subjects.
A) Non-linear least squares fit of the inflection point in vO2, and its 95% confidence interval was calculated using the lsqcurvefit function in matlab. The slope before and after the inflection point (± s.e.) was calculated using a linear regression model and the fitlm function in matlab. The significance of the difference in slope was given as the tstatistic of the additional slope coefficient (H) for a linear model vO2=W:(1+H) for work rates (W) centered to the inflection point. For subjects 2-5 the exercise was divided in two sessions with a short break inbetween to accommodate requirements for comfort. We are aware that the slope has the unit ml[min] -1 J -1 s -1 , which could be further simplified, however, these are preferred units by convention. B) Median RER at low and high work rates. Statistics calculated with a paired two tailed t-test, tested for normality using the Lilliefors test.

Figure 4 Fitted models for multiple subjects.
A) best fit model (orange lines) for 4 subjects optimizing the parameters maintenance, scaling factor and fat oxidizing capacity (HMR_6911), the value for vO2 max was taken directly from the VO2 max test. The curve in absence of complex I bypass (grey line) shown for reference. B) A comparably poor fit was attained when the complex I capacity parameter (HMR_6921) was left unconstrained.

Figure 5
Sensitivity analysis of the model. All parameters of the multi-tissue model were perturbed by a random value uniformly sampled on the interval +-20% of the initial value. The result of 100 simulations shown above. The black curve shows the reference condition, the gray curves are individual perturbed simulations, and the red curve shows the average of all simulations. Histograms of dvO2 per dW at low (0-55% of wmax) and high (55-100% wmax) work rates.  at high work rates, with and without prior exercise. An ordinary differential equation system was set up with a basal lactate flux of 0.5 mol per h and an additional lactate flux of 1.8 mmol per h from the exercising muscle. The lactate filled a reservoir with an estimated volume of 37 L, corresponding to blood and muscle tissue. The lactate efflux from the reservoir was model as a Michaelis-Menten equation with a kM for lactate of 10.73 mM 1 and a manually fitted vmax of 6.6 mol per h. The oxygen flux in absence of lactate was estimated by an exponential equation from literature 2 as dynamic models have previously shown that the oxygen dynamics during the first few minutes of exercise can be derived from transient changes in phosphocreatine 3 . The oxygen flux in the presence of lactate was calculated by the fraction of oxygen flux derived from lactate and an assumed increase in oxygen expenditure when using lactate as carbon source of 16%, based on the increased oxygen expenditure at rest in the presence of lactate, 1.57 vs 1.35 L per min 2 (pseudo code in Supplementary Methods). A) Simulations (blue lines) show that lactate accumulates in blood and tissues before reaching steady state (blue dashed line) for subjects without prior exercise, in agreement with measurements from literature from n=8 subjects (orange dots) 2 , error bars show s.e.m. The baseline (orange dashed line) is taken a few minutes before the exercises. B) After the initial adjustment given by an empirical formula (black dashed line), there was a slow upward drift in oxygen flux driven by the increased usage of lactate as carbon source, which has lower oxygen efficiency than glycogen. C) Lactate concentrations decreased during the exercise when initial levels were higher than the steady state values due to prior exercise. D) High lactate concentrations increased the oxygen consumption rate, but there was no upward drift since lactate uptake was at full capacity under these conditions.    For values with reference to Brenda 11 the highest value for human was selected. If no value existed for human, values from other mammals were used as indicated. If there were no value available in Brenda a literature search was performed.    In literature there are also turnover rates that have been estimated in sub mitochondrial particles from assumed protein concentrations. These rates can be recalculated to specific activities using a molecular mass of 980 kDa 9 , but they were not included in the analysis since complex I had not been isolated. The rates are 250 s -1 (15.3 U per mg) 17 , 211.54 s -1 (12.9 U per mg) 18 , 166 s -1 (10.2 U per mg) 19 and 50 s -1 (3.1 U per mg) 20 .

Supplementary Discussion
Differential expression of metabolic enzymes is prevalent under paraphysiological conditions. There is often an inverse relationship between enzymes from OXPHOS and glycolysis in muscle in paraphysiological conditions 22 , consistent with a reallocation of the enzyme pool. Studies have shown that chronic hypoxia is accompanied by increased levels of OXPHOS proteins including complex I at the expense of glycolytic enzymes and lactate dehydrogenase, whilst the opposite occurs as the acute adaptation to hypoxia 22 . A different pattern of protein reallocation is seen in studies of hypoxia in skeletal muscle of mouse 23 ; there the activity of Complex I is doubled after 2 days of exposure but returns to baseline after 10 days, potentially due to adaptations that restore oxygen availability on the cellular level. Also other forms of metabolic adaptations have been proposed, e.g. rewiring of the TCA cycle using glutamine as carbon source 23 . Patients suffering from complex I deficiency have a lower VO2max compared with healthy control 24 , consistent with the assumption that complex I is required to support complex I bypass. For these patients, a metabolic bypass of complex I, using a membrane-permeable form of succinate, has been proposed as treatment, and shown effective in complex I deficient blood cells 25 . In muscular dystrophy the levels of OXPHOS proteins are increased 22 , but the oxidative capacity per gram of muscle is decreased and the emission of H2O2 by complex I is elevated 26 . A reduced ability to bypass complex I would be consistent with the impaired oxidative capacity despite increased levels of OXPHOS proteins, as well as with increased reliance on flux through the H2O2 producing complex I. It could also result in elevated membrane potential, which has been liked to increased production of H2O2 27 .

Manual curation of the muscle GEM
The genome scale model HMR 2.0 28 contains 3765 genes associated with over 8000 reactions and over 3000 unique metabolites and was used as a scaffold. Reactions that had previously 29 been identified as non-expressed in proteomics and transcriptomics data of human myocytes were removed from the model. A number of different manual curation steps were also performed as specified below resulting in a manually curated muscle GEM with 2951 genes, barely 7000 reactions and 2760 unique metabolites.
ATP synthesis is at the core of muscle metabolism. The stoichiometry of ATP synthase (HMR_6916) was adjusted from 4 H + to 3 H + per ATP in accordance with literature estimates, which gives a P/O ratio close to 2.5 when oxidizing NADH 30 . A stoichiometry of 4 H + is commonly used to account for the proton expenditure of pumping phosphate back in to the mitochondria, but this reaction is already present in the genome-scale model (HMR_5043). To be able to synthesize ATP the uptake of different carbon sources was enabled. A glycogen consumption reaction was added (EC 2.4.1.1). The reaction described the consumption of one glycosyl subunit and one phosphate to produce one glucose-1-phosphate. Uptake reactions of free fatty acids from the extracellular space to the cytosol were also added.
The ability to handle lactate appropriately is important for the model's performance. The idea of an intercellular lactate shuttle between the cytosol and mitochondria is controversial as the experimental data is conflicting and the thermodynamic feasibility has been questioned 31 . To prevent this opportunity in the model, the mitochondrial lactate dehydrogenase reaction (HMR_4280) was blocked. A hypothetical lactate dehydrogenase reaction using ferricytochrome as electron acceptor (HMR_8514) was also blocked due to lack of experimental evidence. It is listed as a probable D-lactate dehydrogenase based on similarities with yeast in Uniprot 32 . Finally a proton-independent lactate transporter was added. The MCT transporter is electron-neutral, one monocarboxylate anion is transported together with a proton 33 , since HMR 2.0 does not model charge (e.g. for lactate), the default formulation (HMR_5998) that relied on cotransport of protons imposed an artificial proton cost on lactate secretion.
To allow the transition between a reduced GEM and the muscle GEM, reactions involved in metabolism were curated. The glycerol phosphate shuttle was curated, the reaction from sn-glycerol-3-phosphate to DHAP (HMR_0483) was allowed to proceed in the canonical direction 34 . Alternative formulations using FAD (HMR_0482) were blocked since FAD is an enzyme-bound cofactor in the reaction. The same applied for the succinate to fumarate reaction (HMR_8743). For several fatty acid reactions (e.g. HMR_3211) the FAD was misannotated as ubiquinone, which bypassed the requirement of the ETF reaction (HMR_6911). These reactions were blocked together with a similar reaction involving proline (HMR_3838). To prevent other oxygen consuming reactions than complex IV (HMR_6914), reactions consuming oxygen in the cytosol were removed, excluding reactions involved in transport and fatty acid oxidation. Additionally, the reactions involved in NADPH and NADH were curated so that NADPH reactions were only involved in anabolism and NADH reactions in catabolism.

Reduction of the muscle GEM to a small-scale model
The muscle GEM was reduced to a small-scale model by removing reactions that did not carry flux when maximizing ATP synthesis under various conditions. Some additional restrictions were imposed to comply with the canonical reactions involved in ATP synthesis. All subcellular compartments were removed apart from the extra cellular compartment, the cytosol and the mitochondria. All transport between mitochondria and cytosol was blocked apart from the malate-aspartate shunt (HMR_3949, HMR_3825 and HMR_4852), uptake of pyruvate (HMR_4926) exchange of CO2, O2 and H2O (HMR_4922, HMR_4898 and HMR_4888) as well as exchange of phosphate, ATP and protons (HMR_5043, HMR_6328 and HMR_7638).
Water was allowed to freely exchange over the boundary (HMR_9047) but the uptake was blocked for all molecules apart from octanoic acid (HMR_9813), glycogen (HMR_9729) and O2 (HMR_9048) and the production was blocked for all molecules apart from CO2 (HMR_9058) and lactate (HMR_9135). The following simulations were performed in the constrained GEM, and reactions that did not carry flux were removed. Maximizing ATP in the presence of fatty acids and glycogen and after blocking complex I (HMR_6921), glycerol-phosphate shunt (HMR_0483), ATP synthase (HMR_6916) and oxygen uptake (HMR_9048) respectively. The resulting model had 95 reactions (out of which 66 were associated with EC codes) and 76 unique metabolites.

Construction of the multi tissue model
The stochiometric matrixes of three identical muscle models (m1, m2 and m3) were concatenated to a joint matrix by adding a compartment for exchange of metabolites between the models and exchange over the boundary (Figure 8). A few auxiliary reactions were added to the joint model, an ATP hydrolysis reaction in each compartment corresponding to the work performed by this sub model, for m3 this was coined maintenance. A reaction for the joint ATP synthesis by m1 and m2, a reaction corresponding to ventilation, exchanging O2 for CO2, and finally a reaction for the buffering of lactate by bicarbonate system (lactate => CO2 + H2O).

Pseudo code, Pareto plot (Figure1B)
Data: S = stoichiometric matrix (mets × rxns) for a reduced model of intermediary metabolism SA = specific activity vector (rxns) with data for each enzyme (µmol mg -1 min -1 ) Parameters: P = protein constraint g per gdw, an arbitrary number, 0<P≤ 1 Functions: [v f]=linprog(problem), solves a linear programming problem and outputs a vector v with the solution and a scalar f with the value of the objective function.

Output:
The estimated mass required for different flux distributions.

Algorithm:
1. Calculate a vector of weights w in the unit g mmol -1 h -1 from the specific activity values and assuming half saturation. (1) 2. Calculate the maximum ATP synthesis capacity under the specific activity constraint using linprog, where is an all zero vector apart from a 1 at the position of the ATP hydrolysis reaction.
3. For multiple ATP hydrolysis rates (r) on the interval [0 ATPmax], calculate the minimal substrate uptake, where is an all zero vector apart from a 6 at the position of the glycogen uptake reaction and 16 at the position of palmitate uptake.
4. Plot X vs Y, where = ATP hydrolysis rates , = ATP hydrolysis rates (4) For the extreme states, identify their flux vector v, and estimate the protein requirement as Pseudo code, fiber simulation ( Figure 2B) Data: S = stoichiometric matrix (mets × rxns) for a reduced model of intermediary metabolism SA = specific activity vector (rxns) with data for each enzyme (µmol mg -1 min -1 ) P = proteomics concentrations vector (rxns) for each enzyme (g per gdw) Output: The optimal flux distributions for different ATP hydrolysis rates assuming substrate minimization. Functions: [v f]=linprog(problem), solves a linear programming problem and outputs a vector v with the solution and a value f with the value of the objective function. Algorithm: 1. Calculate the vmaxes (mmol per h) of the reactions assuming half saturation 2. Calculate the maximum ATP production capacity under the vmax constraint using linprog, where is an all zero vector apart from a 1 at the position of the ATP hydrolysis reaction.
3. For multiple ATP hydrolysis rates (r) on the interval [0 ATPmax], calculate the minimal substrate uptake, where is an all zero vector apart from a 6 at the position of the glycogen uptake reaction and 8 at the position of octanoic acid uptake.
4. Plot fluxes of interest, e.g. X vs Y where = ATP hydrolysis rates , For the comparison of flux and flux capacity ( Figure 2C), for all i let Yi be the maximum value of vi amongst the iterations in step 3. X = vmax. Plot log of X vs log of Y.
Pseudo code, maximum oxygen uptake simulation ( Figure 2D) Data: S = stoichiometric matrix (mets × rxns) for a reduced model of intermediary metabolism SA = specific activity vector (rxns) with data for each enzyme (µmol mg -1 min -1 ) P = proteomics concentrations vector (rxns) for each enzyme (g per gdw) Output: The oxygen uptake rate at the maximum attainable ATP synthesis rate.
Functions: [ for the following conditions: • OC: unconstrained uptake of octanoic acid.
• GM: unconstrained flux through the added reaction aspartate -> glutamate (vasp-glu), that allows uptake of glutamate in exchange for aspartate. • GMS: as GM with the addition of unconstrained flux through the added reaction fumarate -> succinate (vfum-suc) that allows uptake of succinate in exchange for fumarate. • PGMSOC: as GMS with addition of unconstrained uptake of octanoic acid. S = stoichiometric matrix (mets × rxns) for a connected model with 3 muscle submodels. Estimated Parameters: dwMuscle = muscle mass (kg dry weight) muscleRatio = fraction of type 1 muscle type m2Efficency = enzyme capacity ratio (m2/m1). InternalWork = ATP expenditure required to overcome internal resistance at 0 W. type1tresh = ATP synthesis rate before m2 activation. type2target = target for m2 share of total ATP. HMR_6914 = complex IV max, calculated from proteomics data (mmol gdw -1 h -1 ) HMR_6921 = complex I max, calculated from proteomics data (mmol gdw -1 h -1 ) lactateBuffering = maximum rate of lactate buffering in blood (mol per h) Fitted Parameters: HMR_6911= ETF max (mmol gdw -1 h -1 ) scalingFactor = a scaling factor for the differences in enzyme concentrations amongst subjects. maintainance = basal ATP expenditure for maintenance (mol per h) Measured parameters: VO2max = maximum oxidative capacity (mol per h)

Output:
Flux distributions including exchange fluxes. Functions: [v f]=linprog(problem), solves a linear programming problem and outputs a vector v with the solution and a value f with the value of the objective function.

Notation:
We introduce the following notation for convenience, , is an element of the vector v corresponding to the reaction Y in the sub-model X. The m1 sub-model corresponds to oxidative type I muscle fibers, m2 to glycolytic type II fibers, b corresponds to the blood compartment, and s corresponds to the boundary.
2. Set lower bounds on ATP synthesis for maintenance in m3 and internal work in m1 in the lower bounds vector (lb).
7. Construct an objective function (the vector ) for minimization of oxygen and substrate utilization and to disfavor lactate synthesis and ATP synthesis by m2. [ 9. Plot relevant fluxes, e.g oxygen expenditure against watt (using conversion-factor from ATP).
Pseudo code, steady-state watt max simulation (Figure 5a) Data: See connected model.

Fitted Parameters:
See connected model.

Output:
Maximum ATP synthesis capacity under different conditions.

Functions:
See connected model.
2. Constrain the uptake to the specified substrates (glycogen, glucose, fat) and set arbitrary large parameter values to FAFactor and vO2max to study these conditions.
3. Calculate the maximum ATP synthesis rate using linprog, where is an all zero vector apart from a 1 at the position of the m1 ATPtoWork reaction and the joint ATP reaction.
4. Compare the ATP synthesis capacity for the different substrates.
Pseudo code, watt max simulation with substrate depletion (Figure 5b) Data: See connected model. GlycogenPoolSize = the pool sizes of muscle glycogen (mol). GlucosePoolSize = the pool sizes of liver glucose (mol). LactatePoolSize = the pool of lactate based on the maximum lactate concentration in blood (mol).

Fitted Parameters:
See connected model.

Output:
Maximum ATP synthesis capacity for different exercise durations.

Functions:
See connected model.

2.
To avoid ambiguity of the solutions set up a secondary objective function to disfavor lactate production and ATP synthesis by m2. Output: The predicted time dynamics of lactate and oxygen at different work rates. Algorithm: 1. Calculate basal lactate flux = (basalConcentration, K) 2. Solve an ODE problem for the accumulation of lactate under this condition.
3. Calculate the expected oxygen consumption from consumption of lactate as 4. Calculate oxygen consumption based on empirical formula from literature 2 : 6. Calculate the modified oxygen consumption rate as: