Metabolic flux analysis (MFA) is an important component of metabolic engineering (Stephanopoulos, 1999). MFA measures the distribution of flux in the different metabolic pathways of an organism as well as identifies the intracellular flux changes associated with altered cellular response such as overproduction of a desired metabolite for biotechnological applications or understanding the altered pathways in diseases. The flux values from MFA can be used to constrain the solution space of the linear optimization problem in flux balance analysis to give a physiologically meaningful flux distribution (Kim & Reed, 2012; Desai & Srivastava, 2015).
Metabolic flux analysis methodology has evolved over years. Earlier studies employed an approach termed stoichiometric MFA which balanced the fluxes around the intracellular metabolites in a stoichiometric network (vanGulik et al., 2001; Srivastava & Chan, 2008). While it provided important insights into altered metabolism under different conditions, there are several limitations of stoichiometric MFA such as the inability to determine fluxes in reversible or parallel reactions (Wiechert et al., 2001). 13C-MFA, which employs 13C-labeled substrates can resolve fluxes through such reactions. For a typical 13C-MFA, the cells are fed with a labeled substrate or mixture of labeled substrates, labeled at one more carbon atoms. The cells are required to be at metabolically steady state, i.e., when the intracellular concentrations of metabolites do not change over time (Zamboni et al., 2009). This ensures that the differences in the labeling patterns are not due to temporal variations in the fluxes, but due to the differential activity of the pathways. Once the labeled substrate attains steady state distribution throughout the metabolic network (isotopic steady state), the relative abundances of the mass isotopomers of intracellular metabolites and proteinogenic amino acids are measured using NMR (Szyperski et al., 1999) or mass spectrometry to derive their respective mass isotopomer distributions (MIDs). The gas chromatography/mass spectrometry technique for measurement of MIDs of tert-butyldimethylsilyl derivatives of amino acids is well established. The fragmentation pattern produced by the derivatives is also well characterized and these ions are used to calculate the MIDs which are eventually used to calculate the fluxes (Antoniewicz, Kelleher & Stephanopoulos, 2007b; Leighty & Antoniewicz, 2013; Crown, Long & Antoniewicz, 2015; Gonzalez, Long & Antoniewicz, 2017). Zamboni et al. provide a detailed protocol for 13C-MFA (Zamboni et al., 2009). A suitable metabolic model and the MID data are then fed in software tools which estimate the flux distributions in reactions included in the MFA model.
Along with the development of experimental techniques for stationary MFA, the methodology to calculate flux distribution from the labeling data has also evolved over the years. An important development in this direction was the introduction of the elementary metabolite unit (EMU) framework (Antoniewicz, Kelleher & Stephanopoulos, 2007a). The EMU decomposition algorithm detects minimum number of isotopomer balances needed to simulate the observed labeling pattern which greatly reduces the number of variables without any loss of information, thus making it much more efficient than the isotopomer mapping matrices approach (Schmidt et al., 1997). Based on the measured fractional enrichments of metabolites, the metabolic network is decomposed into elementary reactions of the EMUs. These elementary reactions form the basis of system of equations representing the relation between metabolic fluxes and the observed isotopic labeling. Recent MFA software such as WUFlux (He et al., 2016), Metran (Antoniewicz, Kelleher & Stephanopoulos, 2007a) and OpenFLUX (Quek et al., 2009) employ the EMUs framework. Most software for 13C-MFA such as OpenFLUX (Quek et al., 2009), WUFlux (He et al., 2016), INCA (Young, 2014) and Metran (Antoniewicz, Kelleher & Stephanopoulos, 2007a) work either on the proprietary MATLAB platform or are not free and open source. Many also use the MATLAB statistics and optimization toolboxes which are proprietary. Thus, a truly free and open-source software is not available yet. Understanding the requirement for open source tools, Birkel et al. have recently developed jQMM, an open source Python library to perform 13C-MFA which uses the general algebraic modeling system (GAMS) to solve the optimization problem (Birkel et al., 2017). However, GAMS is a proprietary system for mathematical optimization. We believe that a truly open-source software in which all the components and packages are free and open source would make MFA more accessible to researchers and also encourage further development of 13C-MFA and related tools.
Materials and Methods
FluxPyt is written in Python (version 3.6.1) because Python is an open-source platform and has free, well-developed data handling packages with ability to perform the various tasks in 13C-MFA. In FluxPyt, the numerical and matrix operations are conducted using the packages NumPy and SciPy (van der Walt, Colbert & Varoquaux, 2011), and SymPy (Meurer et al., 2017). The non-linear optimization (see Eq. 3 below) is conducted using the package lmfit (Newville et al., 2014). Python was installed as Anaconda3 (version 188.8.131.52) distribution. Anaconda3 also acts as a package manager to handle the dependency and compatibility issues between the different packages necessary for the functioning of FluxPyt. FluxPyt was developed and tested using the Spyder integrated development environment (IDE) with the IPython console. FluxPyt was tested in Windows 7, Windows 10 and Linux Mint (18.2) operating systems (OS). Detailed information of packages in the Anaconda distribution is provided in Data S1. While FluxPyt is designed to be used by people with limited programming skills, the open source nature allows further modification by advanced users. The general framework of FluxPyt is shown in Fig. 1. Table 1 describes different modules present in FluxPyt. The MFA model, isotope correction data and substrate input data are written in a csv file which can be easily edited in a spreadsheet. The reactions are divided into different types, namely, irreversible (“F”), reversible forward (“FR”), reversible reverse (“R”), metabolite balancing reactions (“B”) and isotopomer balancing reactions (“S”). The metabolites which do not take part in isotopomer balancing are marked as “X” in the atom transition equations. The basis flux, which typically is the substrate uptake flux in a network, is marked with a “*”. The upper bounds of other reactions are limited to 15 times this basis flux. The reverse free fluxes cannot be determined as they have no natural upper limit (Quek et al., 2009). The substrate information and the natural correction data are read from separate csv files which must be named substrate_input.csv and corr_file.csv, respectively. The correction file also specifies the number of mass isotopomer measurements for every amino acid provided for fitting. This is used to extract the MID vectors of length equal to the length of MID vectors provided in the measurement data. For detailed information on creating an MFA model for FluxPyt, a user guide is provided in the supplemental material (Data S2).
|main.py||Main driver module. Asks user inputs|
|build_model.py||Parses csv model file. Generates metabolite and isotopomer models, stoichiometric matrix. Reads measured MIDs|
|create_atm_transition_equations.py||Creates list of elementary reactions from atom transition equations|
|input_substrate_emu.py||Reads substrate input file and calculates substrate MIDs required in the emu network|
|mid_vec_gen.py||Generates correction vector for a given natural isotope distribution vector and number of atomic elements|
|priliminary_fba.py||Generates initial feasible flux distribution using glpk as linear programming solver and sets feasible reaction bounds|
|glpk_solve.py||Solves linear programming problems using glpk|
|make_emu_networks.py||Makes EMU networks from elementary reactions and measured MIDs|
|free_flux_null.py||Calculates null matrix and free fluxes|
|solve_mfa_lm.py||Generates and solves non-linear optimization problem using lmfit|
|solve_mid_networks.py||Solves EMU network to find calculated MIDS with respect to parameters of the non-linear problem|
|mid_corr.py||Corrects calculated MIDs|
|draw_flux_map.py||Draws flux map from optimization result and flux map template|
|monte_carlo.py||Performs Monte-Carlo analysis|
|bootstrap.py||For bootstrapping results from Monte-Carlo analysis|
The model is automatically parsed into metabolite and isotopomer models. The metabolite model is a mass balance model that satisfies the metabolic steady state condition, (1) where, S is the stoichiometric matrix and v is the flux distribution vector. For a given stoichiometric matrix, the flux distribution vector can be written as a function of free fluxes (vfree) as follows: (2) where N is the null space of stoichiometric matrix.
The isotopomer balance model is generated based on the EMU decomposition algorithm (Antoniewicz, Kelleher & Stephanopoulos, 2007a). It assumes an isotopic steady state and is used to solve the non-linear optimization problem using the free fluxes as the optimization parameters to minimize the sum of squared residues (SSR) as given in Eq. (3).(3)
Levenberg–Marquardt algorithm implemented in the lmfit package is used to optimize the free flux parameters in the metabolite model. The free flux assignments are done based on the reduced row echelon form (rref) of the stoichiometric matrix. Before calculating the rref, the reactions are arranged in following order: and are the forward and reverse reactions of reversible reactions; are irreversible reactions; are reactions that the user chooses to act as free fluxes; are reactions with measured or predetermined values. The arrangement of reactions in a specified order ensures that the measured fluxes are calculated as free fluxes. In addition, the user may provide a preference for certain reactions to be calculated as free fluxes (). The default number of iterations is 10, but this value can be changed by the user. The iterations are run with different random initial parameters so as to increase the chances of finding global minimum.
Flux standard deviations are calculated using the Monte-Carlo method. Multiple datasets of the measured MID values are generated using the mean and standard deviations of the respective MIDs. The default value for the number of datasets is 500, but can be changed by the user. Each dataset of the measured MID values is then used to solve the non-linear optimization problem (He et al., 2016). This results in multiple optimal solutions (500 solutions by default) for each reaction. From these values the 68% and 95% confidence intervals (CIs) are determined from the upper and lower 16% and 2.5% values, respectively, using the bootstrap method (He et al., 2016) as follows. One thousand samples of 500 values each (if the default number of samples for Monte-Carlo analysis is used) are used to calculate the required percentiles (2.5, 16, 50, 84 and 97.5 percentiles). The calculated percentiles from each sample are stored as a separate set. The median of each percentile value is reported as the upper or lower boundary of the respective CI. For example, the median 2.5 and 97.5 percentile values represent the lower and upper boundaries of the 95% CI and the median 16 and 84 percentile values represent the lower and upper boundaries of the 68% CI, respectively. The median of the 50 percentile values represents the median of the estimated fluxes.
Goodness of fit
Assuming that the model is correct and the data are without gross errors, the SSR has a χ2 distribution with number of degrees of freedom equal to number of measurements to be fitted (n) minus the number of independent parameters (p) (Antoniewicz, Kelleher & Stephanopoulos, 2006; Leighty & Antoniewicz, 2013). The acceptable range of SSR is between and , where α is the chosen significance level, for example, α = 0.05 for 95% CI.
Output files and visualization of the results
FluxPyt generates several result files during the MFA, a list of which is provided in Table S1. The data results, containing the calculated MIDs, the optimal flux distribution values and CIs are stored as separate .csv and .pckl files. FluxPyt provides graphical representation of results. The differences between the measured and the calculated MIDs of the best fit solution are automatically plotted. Flux map is automatically generated using a template flux diagram created using Inkscape. The template for flux map is a svg file that can be created using Inkscape. Templates for tricarboxylic acid cycle (TCA) and Corynebacterium glutamicum networks are provided in the additional data (Data S1) or can be downloaded from the FluxPyt project site. The user can easily modify the template based on their model of interest using the same drawing tool to make flux maps for their organism of interest. The users need to write the ID of the reactions at the appropriate position in the template where they want the flux value of that particular reaction to be printed (Data S1). The bootstrap results are plotted as box plots. The user can choose the specific reactions to be plotted. Additionally, the bootstrap output is stored as a pandas dataframe (McKinney, 2010).
We employed the models and the MID data given in OpenFLUX (Quek et al., 2009) for initial evaluation of FluxPyt. Both the TCA cycle and the C. glutamicum central metabolic network (Data S1) were compared.
Preliminary validation with the TCA cycle model
The TCA cycle metabolic network (see Data S1 for a detailed description of the model as well as a list of metabolites) consisted of 18 metabolite reactions and three isotopomer balancing reactions (type “S” reactions). The network has two substrate input reactions, pyruvate (V01) and glutamate (V02) and one reversible reaction (represented as two separate reactions, V07 and V08). Three of the reactions (V16, V17 and V18) have measured values relative to pyruvate flux fixed at 1. The stoichiometric matrix has seven free variables. Rearranging the reactions as mentioned in the implementation section above enables the measured flux values to be calculated as free fluxes. The three measured free flux parameters were allowed to vary within three standard deviations of their means while the remaining three free flux parameters were allowed to vary between the maximum bounds as explained in the implementation section. The only remaining parameter was the flux through pyruvate uptake reaction which was fixed to 1. A total of 18 data points were fitted to six flux parameters, therefore the maximum acceptable SSR was 21. The number of iterations was set to 10. Flux analysis results from both software converged to a similar solution. The MIDs calculated by FluxPyt were the same as those calculated by OpenFLUX (Table 2). Very minor differences in the values of two fluxes were observed (Table 3). These could be due to the different solvers and algorithms used to solve the least squares fitting problem and linear system of equations (He et al., 2016). The flux map, which was automatically generated (see Methods for details) based on the template, is shown in Fig. 2. The Monte-Carlo analysis suggested tight CIs for the TCA model (Fig. 3). The reversible reaction spanned the entire flux range, as expected. The optimal values for all fluxes lay within 2.5 and 97.5 percentile range, i.e., the 95% CI (Table 4).
|Sum of squared residual (SSR)||0.0002||0.0002|
|Reaction ID||Reaction formula||FluxPyt||OpenFLUX|
|V01||PYR_EX -> PYR||1.00||1.00|
|V02||GLU_EX -> AKG||0.20||0.20|
|V03||PYR -> ACCOA + CO2 + NADH||0.62||0.62|
|V04||ACCOA + OAA -> ICI||0.62||0.62|
|V05||ICI -> AKG + CO2 + NADH||0.62||0.62|
|V06||AKG -> 0.5 SUC + 0.5 SUC + CO2 + NADH + ATP||0.59||0.59|
|V07||SUC -> OAA + FADH2 + NADH||2.62||2.61|
|V08||OAA + FADH2 + NADH -> 0.5 SUC + 0.5 SUC||2.02||2.02|
|V09||OAA -> PYR + CO2||0.30||0.30|
|V10||PYR + CO2 + ATP -> OAA||0.48||0.48|
|V11||2 NADH + O2 -> 4 ATP||1.21||1.21|
|V12||2 FADH2 + O2 -> 2 ATP||0.30||0.30|
|V13||O2_EX -> O2||1.51||1.51|
|V14||CO2 -> CO2_EX||1.65||1.65|
|V15||ATP -> ATPM||5.56||5.54|
|V16||PYR -> PYR_B||0.20||0.20|
|V17||AKG -> AKG_B||0.23||0.23|
|V18||OAA -> OAA_B||0.15||0.15|
|Reaction ID||−95% CI||−68% CI||Median||+68% CI||+95% CI|
MFA of C. glutamicum metabolic network
The C. glutamicum model (Quek et al., 2009) consisted of 67 reactions, 42 metabolites and 25 free fluxes (see Data S1 for a detailed description of the model as well as a list of metabolites). The flux measurements for 18 reactions were included in the model. The 24 variable parameters were fit to the 29 data points of amino acid labeling pattern. Here too, the calculated fluxes were similar to those calculated by OpenFLUX (Table S2).
The SSR values for the measured vs. calculated MIDs were similar for FluxPyt (690) and OpenFLUX (701) (Table 5). However, the SSR values for both the software were higher than the maximum acceptable SSR value (11.07), suggesting that either the model used was incorrect or there were measurement errors. Use of a recent MID data (Shupletsov et al., 2014), reduced the SSR to 10.33 (Table 6), which is within the acceptable range (1.15–11.07). This showed that the MID data previously used had measurement errors. The Fig. 4 shows the flux map for the C. glutamicum model for this data (see Table S3) for the list of flux values. The CIs calculated with Monte-Carlo method are shown in Table S4.
|Measured MIDs (Quek et al., 2009)||FluxPyt||OpenFLUX|
Our results show that FluxPyt can be used for performing MFA with real metabolic models. It is capable of flux estimation when multiple labeled substrates are used (e.g., the TCA model uses three labeled substrates, namely 1-13C pyruvate, U-13C pyruvate and 1-13C Glutamate).
Metabolic flux analysis is a powerful method to measure intracellular fluxes. However, all of the currently used software have some component that is proprietary and thus makes it impossible to conduct 13C-MFA without obtaining/purchasing a license. Most licenses do not allow altering the code, restricting further development. In designing this open-source software, our effort was to keep the interface simple and easy to use for scientists somewhat familiar with other MFA software. Therefore, the files for labeling-data preparation have been kept similar to the popular software OpenFLUX. The software can be run through different interfaces. The Anaconda command prompt may be convenient for people familiar with Python, while the Spyder IDE would appear more familiar to scientists familiar with MATLAB interface. The IPython console in Spyder offers a more interactive environment than the command prompt. A print of FluxPyt for performing MFA of TCA network running in the Spyder console is provided in Data S3.
The fluxes calculated by MFA are dependent on the network chosen and the measurements used. In case of the C. glutamicum metabolic network, it was found that the MID data used by Quek et al. (2009) had measurement errors giving unacceptably high SSR for the MIDs calculated. Use of more recent MID data (Shupletsov et al., 2014) significantly reduced the SSR to acceptable values.
At this point the software is released for use in Microsoft Windows and Linux OS. Similarly, the CI calculations currently use the Monte-Carlo method. The open-source nature makes it possible for advanced users to adapt it to other OS as well as add other methods for calculating CIs.
FluxPyt is free and open-source software for conducting stationary 13C-MFAanalyses with performance comparable to other available software. It can be used for flux estimation with experimental data of practical scale. Automatic generation of flux maps allows better visualization of results. This tool will be useful to researchers interested in stationary 13C-MFA studies.