OpenFLUX: efficient modelling software for 13C-based metabolic flux analysis

Abstract

Background: The quantitative analysis of metabolic fluxes, i.e., in vivo activities of intracellular enzymes and pathways, provides key information on biological systems in systems biology and metabolic engineering. It is based on a comprehensive approach combining (i) tracer cultivation on 13C substrates, (ii) 13C labelling analysis by mass spectrometry and (iii) mathematical modelling for experimental design, data processing, flux calculation and statistics. Whereas the cultivation and the analytical part is fairly advanced, a lack of appropriate modelling software solutions for all modelling aspects in flux studies is limiting the application of metabolic flux analysis.

Results: We have developed OpenFLUX as a user friendly, yet flexible software application for small and large scale 13C metabolic flux analysis. The application is based on the new Elementary Metabolite Unit (EMU) framework, significantly enhancing computation speed for flux calculation. From simple notation of metabolic reaction networks defined in a spreadsheet, the OpenFLUX parser automatically generates MATLAB-readable metabolite and isotopomer balances, thus strongly facilitating model creation. The model can be used to perform experimental design, parameter estimation and sensitivity analysis either using the built-in gradient-based search or Monte Carlo algorithms or in user-defined algorithms. Exemplified for a microbial flux study with 71 reactions, 8 free flux parameters and mass isotopomer distribution of 10 metabolites, OpenFLUX allowed to automatically compile the EMU-based model from an Excel file containing metabolic reactions and carbon transfer mechanisms, showing it's user-friendliness. It reliably reproduced the published data and optimum flux distributions for the network under study were found quickly (<20 sec).

Conclusion: We have developed a fast, accurate application to perform steady-state 13C metabolic flux analysis. OpenFLUX will strongly facilitate and enhance the design, calculation and interpretation of metabolic flux studies. By providing the software open source, we hope it will evolve with the rapidly growing field of fluxomics.

Figure 1

Workflow of OpenFLUX. The software consists of 2 components (shaded box) – a JAVA-based metabolic model parser and a set of MATLAB-based algorithms for parameter estimation and sensitivity analysis. The user defines the reaction network in the model definition file, which contains the metabolic model and the experimental data input. The parser reads the metabolic model, and subsequently generates a metabolite balance model (stoichiometric matrix), and an isotopomer balance model using the EMU framework. The MATLAB-based computational algorithms use these models, in conjunction with the experimental data, to perform least-square parameter estimation and parameter sensitivity analysis. The optimized flux distributions and confidence intervals are obtained as results output.

Figure 2

Algorithm for weighted least-square parameter estimation. OpenFLUX uses MATLAB's FMINCON to minimize the residual error () between the experimental MIDs () and the simulated MIDs (), by optimizing the free fluxes (), which are subjected to lower and upper boundary value constraints. is weighted by the variance matrix (D). The inverse compactification function is used to transform reverse flux from a numerical parameter (v^{←, [0,1]}) into a physical flux parameter (v^←). The null-space matrix (Ns) is then used to calculate the flux vector () from the free fluxes (). Finally, is calculated using and known MIDs of the input substrates (). FMINCON terminates when a local minimum is found, revealing the optimum free fluxes.

Figure 3

Flux distributions for a simplified TCA cycle model as described in the text and Appendix. Fluxes were calculated using OpenFLUX (top values) and the cumomer package 13C-FLUX [22] (bottom values). The absolute forward fluxes are displayed in the circles. The free fluxes (R2, R8, R9) and their 95% CI (generated by OpenFLUX) are displayed in the rectangular boxes. Full arrowhead: forward flux. Line arrowhead: reverse flux. All reactions are unidirectional except for the SUC↔OAA reaction. Suffix: _B, biomass drain; _EX, exo-metabolites. Metabolites: PYR, pyruvate; ACCOA, acetyl-CoA; CIT, citrate/isocitrate; AKG, α-ketoglutarate; SUC, succinate; OAA, malate/oxaloacetate; GLU, glutamate. Reactions: R1, pyruvate uptake; R2, glutamate take; R3, pyruvate dehydrogenase; R4, citrate synthase; R5, iso-citrate dehydrogase; R6, α-ketoglutarate dehydrogenase & succinyl-CoA hydrolase; R7/R8; fumarate hydratase & malate dehydrogenase; R9, malic enzyme; R10, pyruvate carboxylase; R11/R12, oxidative phosphorylation; R13, oxygen uptake; R14, CO2 evolution; R15, ATP maintenance; R16, pyruvate biomass drain; R17, α-ketoglutarate biomass drain; R18, oxaloacetate biomass drain.

Figure 4

Overview of the network model definition text file. The model definition file consists of one table and five lists. The metabolic reaction network is defined in the table, and is organized under the following headings: reaction ID (rxnID), reaction equation (rxnEQ), reaction atom transition (rxnCTrans), reaction rate (rates), reaction type (rxnType), free flux allocation (basis) and flux value standard error (deviation). If a given reaction rate is known (e.g., biomass drain rates), then the flux value and the corresponding measurement error can be included in the basis and deviation columns respectively. The parser requires the user to separately list down metabolites that are excluded from the stoichiometric model (excludedMetabolites), EMUs that are to be calculated by the isotopomer model (simulatedMDVs), and input substrates that contribute to the isotopomer balance (inputSubstrates). The experimental MIDs (measurements) and the associated errors (error) are listed in the same order as the EMUs in the simulatedMDVs list. The figure insert (bottom right) is an example of how pyruvate dehyrogenase reaction is described in the rxnEQ and rxnCTrans columns. CO₂ is produced by cleaving the carboxylic end (C1) of pyruvate, leaving the acetyl moiety (C2–C3). Alphabet letters in the atom transition equation is used to represent the transfer of first carbon atom of pyruvate to CO₂, and the second and third atoms of pyruvate to the first and second atoms of acetyl moiety. An exception is the letter "x" or "X", which is used to exclude metabolite from the isotopomer balance, such as NADH.

Figure 5

The null-space matrix generated from the stoichiometric model. The structure of this NS reflects the generic NS matrix form shown in Eq. 3. The Identity matrix indicate the 1-to-1 mapping of the 9 free flux parameters – ν₈, v₂, v₉ are the unknown parameters, while ν₁, ν₁₉, ν₂₀, ν₁₆, ν₁₇ and ν₁₈ are the known parameters.

Figure 6

Coordinate systems used to describe bi-directional reactions. The schematics under the matrix equations shows how the two coordinate systems, [v^xch, v^net] and [v^←, v^{→, net}], are mapped to [v^←, v^→ ]. The former requires a conditional operator, which cannot be implemented using the null-space matrix. The prerequisite for the compactification transformation is that the black and grey arrows are overlapping.