BayFlux: A Bayesian method to quantify metabolic Fluxes and their uncertainty at the genome scale

Metabolic fluxes, the number of metabolites traversing each biochemical reaction in a cell per unit time, are crucial for assessing and understanding cell function. 13C Metabolic Flux Analysis (13C MFA) is considered to be the gold standard for measuring metabolic fluxes. 13C MFA typically works by leveraging extracellular exchange fluxes as well as data from 13C labeling experiments to calculate the flux profile which best fit the data for a small, central carbon, metabolic model. However, the nonlinear nature of the 13C MFA fitting procedure means that several flux profiles fit the experimental data within the experimental error, and traditional optimization methods offer only a partial or skewed picture, especially in “non-gaussian” situations where multiple very distinct flux regions fit the data equally well. Here, we present a method for flux space sampling through Bayesian inference (BayFlux), that identifies the full distribution of fluxes compatible with experimental data for a comprehensive genome-scale model. This Bayesian approach allows us to accurately quantify uncertainty in calculated fluxes. We also find that, surprisingly, the genome-scale model of metabolism produces narrower flux distributions (reduced uncertainty) than the small core metabolic models traditionally used in 13C MFA. The different results for some reactions when using genome-scale models vs core metabolic models advise caution in assuming strong inferences from 13C MFA since the results may depend significantly on the completeness of the model used. Based on BayFlux, we developed and evaluated novel methods (P-13C MOMA and P-13C ROOM) to predict the biological results of a gene knockout, that improve on the traditional MOMA and ROOM methods by quantifying prediction uncertainty.


Term Description Cyclic flux
Simultaneous forward and backwards flux within a single (reversible) metabolic reaction.Quantified as f orward f lux − backwards f lux.Sometimes referred to as an "exchange flux" in 13 C MFA literature [1], but we avoid using the word exchange to avoid confusion with extracellular exchange fluxes.(Extracellular) exchange flux Flux of reactions which exchange metabolites outside of the organism or cell, for example the uptake of an extracellular nutrient.Often, in a liquid cell culture, this can be measured directly by a change in extracellular concentration over time.

Edge points
Points on the edge of the feasible flux polytope, used by the ACHR (Artificial Centering Hit and Run algorithm) [2] as reference points to determine movement directions for the sampler.These are typically the maximum and minimum fluxes for each reaction in the metabolic model, as determined through linear optimization.Often called "warm-up points" in ACHR literature, we use the word 'edge' to avoid confusion with MCMC (Markov Chain Monte Carlo) warm-up samples.Warm-up samples Samples collected during the initial period of a MCMC (Markov Chain Monte Carlo) run, during which the sampler is unlikely to have yet reached a stationary distribution.These samples are typically discarded when analyzing the posterior probability distribution.For the CBMKr (Carbamate kinase) flux, the probability distributions from P-13 C ROOM is much closer to the experimentally obtained probability distribution for the knockout (computed from experimental data using BayFlux) indicating a more accurate prediction.The probability distribution from P-13 C MOMA is only marginally closer to the experimentally obtained probability distribution for the knockout since the peak of the distribution is a bit closer than the single value FBA MOMA prediction.Here, we show only one reaction as an illustration of the approach, whereas the true demonstration that the predictions are better must involve all fluxes and is shown in Fig. 7 in the main text.

Fig A .
Fig A.Fit labeling patterns are the same for both approaches, BayFlux and 13CFLUX2.Mass distributions vectors (MDVs) for experimental data (grey), the fits for the best sample for BayFlux (dark blue, ten million samples), and the fits for 13CFLUX2 (red) are closer together than the standard deviation of the experimental error.Units for horizontal axis are the number of extra neutrons in the measured metabolite[3].

Fig B.
Fig B.Flux results for BayFlux and 13CFLUX2 are essentially the same for fifteen different inputs.The best BayFlux sample (dark blue, ten million samples, y axis) and its credible interval overlap with the 13CFLUX2 best fit (in red, x axis) and its confidence interval for all fifteen different inputs and fluxes.The fifteen different inputs were obtained by randomly changing the exchange fluxes, but keeping the same labeling profiles for the metabolites.Both axes are in units of mmol/gDW/h.We set the 13CFLUX2 confidence bounds that could not be determined to the overall maximum BayFlux bounds for comparison purposes.

Fig C .
Fig C. 13CFLUX2 vs. BayFlux for 15 different inputs: Best Sample out of 10 million from BayFlux (in dark blue) vs. 13CFLUX2 (in red).We set the 13CFLUX2 confidence bounds that could not be determined to the overall maximum BAYFLUX bounds for comparison purposes.

Fig
Fig D. P-13 C MOMA and P-13 C ROOM work by first computing the WT flux profile distribution using BayFlux (in orange).Each dot corresponds to a flux profile (i.e.reaction fluxes for all reactions in a metabolic model) in the flux phase space.A representative set of flux profiles for this base flux profile distribution is obtained through sampling (orange dots).For each of these flux profiles, a new flux profile is computed by using MOMA (in blue) or ROOM (in green) to predict the resulting flux after a knockout.The fact that BayFlux produces a flux profile distribution results in MOMA and ROOM predicting distributions of flux profiles.