OpenFLUX2: 13 C-MFA modeling software package adjusted for the comprehensive analysis of single and parallel labeling experiments

Steady-state (13)C-based metabolic flux analysis ((13)C-MFA) is the most powerful method available for the quantification of intracellular fluxes. These analyses include concertedly linked experimental and computational stages: (i) assuming the metabolic model and optimizing the experimental design; (ii) feeding the investigated organism using a chosen (13)C-labeled substrate (tracer); (iii) measuring the extracellular effluxes and detecting the (13)C-patterns of intracellular metabolites; and (iv) computing flux parameters that minimize the differences between observed and simulated measurements, followed by evaluating flux statistics. In its early stages, (13)C-MFA was performed on the basis of data obtained in a single labeling experiment (SLE) followed by exploiting the developed high-performance computational software. Recently, the advantages of parallel labeling experiments (PLEs), where several LEs are conducted under the conditions differing only by the tracer(s) choice, were demonstrated, particularly with regard to improving flux precision due to the synergy of complementary information. The availability of an open-source software adjusted for PLE-based (13)C-MFA is an important factor for PLE implementation.The open-source software OpenFLUX, initially developed for the analysis of SLEs, was extended for the computation of PLE data. Using the OpenFLUX2, in silico simulation confirmed that flux precision is improved when (13)C-MFA is implemented by fitting PLE data to the common model compared with SLE-based analysis. Efficient flux resolution could be achieved in the PLE-mediated analysis when the choice of tracer was based on an experimental design computed to minimize the flux variances from different parts of the metabolic network. The analysis provided by OpenFLUX2 mainly includes (i) the optimization of the experimental design, (ii) the computation of the flux parameters from LEs data, (iii) goodness-of-fit testing of the model's adequacy, (iv) drawing conclusions concerning the identifiability of fluxes and construction of a contribution matrix reflecting the relative contribution of the measurement variances to the flux variances, and (v) precise determination of flux confidence intervals using a fine-tunable and convergence-controlled Monte Carlo-based method.The developed open-source OpenFLUX2 provides a friendly software environment that facilitates beginners and existing OpenFLUX users to implement LEs for steady-state (13)C-MFA including experimental design, quantitative evaluation of flux parameters and statistics.

section_content: Metabolic flux analysis (MFA) plays a key role in systems biology because intracellular fluxes, i. e., in vivo reaction rates through different pathways within an intact living cell [1], are the functional output of all conventional genetic and metabolic regulatory systems and determine the physiological phenotype of the cell [2]. 

In recent decades, a metabolic steady-state version of 13 C-labeling-based MFA ( 13 C-MFA) has become the best developed and most powerful method for quantifying intracellular fluxes when all fluxes and metabolite concentrations can be considered to be (at least approximately) constant [2] [3] [4]. 13 C-MFA, applied to microbial, plant and mammalian systems, has been increasingly used in systems biology and metabolic engineering [5] [6] [7], biotechnology and medicine [8] [9] [10]. 

Due to the high complexity of native metabolic networks, 13 C-MFA typically involves the use of a simplified stoichiometric model in which only the key pathway reactions of the central carbon metabolism and the set of lumped targeted biosynthetic reactions are parameterized before the assumed model-based fluxes are inferred from measurable quantities [11]. 

Concerning the experimental data applied in 13 C-MFA, the physiological/extracellular fluxes or effluxes (e. g., biomass precursor drain, substrate uptake, and product excretion rates) and 13 C-labeling patterns (i. e., isotopomer distributions) of metabolic products resulting from feeding partially 13 C-labeled substrates (tracers) are used. These effluxes are determined from the time courses of cellular dry weight and extracellular metabolite concentrations during cultivation [12]. The 13 C-isotopomers generated due to the metabolic conversion of tracers are detected through nuclear magnetic resonance (NMR) spectroscopy [13], mass spectrometry (MS) [14], and/or tandem MS (MS/MS) [15]. 

Prof. W. Wiechert and co-workers significantly contributed towards formalizing the framework for 13 C-MFA: from measured effluxes and intracellular labeling information the intracellular fluxes could be computed [16] [17] [18] [19]. On this basis, several mathematical models have been developed that can simulate a unique profile of isotopomer abundance for the fluxes with assigned parameters by describing the propagation of labeled atoms from the tracer through an assumed metabolic network according to the known atom rearrangements for each reaction [18, [20] [21] [22] [23] [24] [25] [26]. All simulations are providing under the essential assumption that the possible isotopic mass effects [27] are negligible, i. e., that the labeling states of the metabolites do not influence the rate of their enzymatic conversion [16]. The goal of 13 C-MFA is to determine the set of initially unknown flux parameters that minimizes the differences between experimentally observed and simulated measurements. In mathematical essence, this set is a solution of a large-scale non-linear parameter estimation problem [28]. Analytical solutions of this problem are available only for the simplest systems. Therefore, the values of the assumed fluxes must generally be inferred from the experimental datasets through computer model-based interpretation using an iterative least-squares fitting procedure [2, 3, 26]. 

Several high-performance computational software suites for performing flux calculations have been developed and described, e. g., 13CFLUX [29] and its reinforced version -13CFLUX2 [30], METRAN [31, 32], OpenFLUX [33], FIA [25], influx_s [34], OpenMebius [35]. 

These software toolboxes most often automatically generate metabolite and isotopomer balance models relying on an initially user-defined simple notation of metabolic networks and the known atom transitions occurring in biochemical reactions. Then, starting from the generated models and from measured effluxes that must be constrained within the obtained error ranges, semi-random guesses regarding intracellular fluxes are used to simulate in silico 13 C-labeling patterns of targeted metabolites, which, in turn, are compared with the measured patterns. This process is repeated until a satisfactory match to the measurable quantities is achieved, i. e., the constrained non-linear least-squares minimization problem (NLLSP) is solved [29]. According to the rules of regression analysis, providing a statistical goodness-of-fit test of the adequacy of the applied flux model is required, at a minimum, after determining the optimized fluxes [28, 36]. Then, linearized statistics [17, 37, 38], a non-linear-based search algorithm [28], and/or the Monte Carlo approach [39, 40] are used to estimate the precise flux resolution, i. e., the uncertainty of the determined fluxes. The optimized parameters of the fluxes and their confidence intervals in the statistically adequate user-made metabolic model must be obtained as the concerted results of these computations. 

When the 13 C-labeling data were obtained from NMR in the early stages of 13 C-MFA development, each analysis was typically performed on a single labeling experiment (SLE), primarily for cost reasons [41]. Implementing highly sensitive MS- [42] [43] [44] and MS/MS-mediated [15, 45] measurements, which are development approaches that involve 13 C-tracer experiments at a miniaturized scale [46, 47], led to a significantly increased accessibility and decreased cost of labeling experiments (LEs). Thus, it has become possible to realize the advantages of parallel labeling experiments (PLEs), in which two or more LEs are initiated from the same seed culture and conducted in parallel under the same experimental conditions differing only in the set of 13 C-tracers applied [48] [49] [50] [51] [52] [53]. 

SLE-based 13 C-MFA remains a widely used method and can be implemented with the application of a single labeled substrate as a tracer or using a mixture of isotopomers of the same compound or multiple labeled substrates [32, [54] [55] [56] [57]. Studies have shown that achieving optimal resolution of fluxes from different parts of the central carbon metabolism requires different 13 C-tracers [19, 48, 58]. Several sophisticated experimental design strategies have been adopted to improve the desired flux precision [19, 32, 36, 49, 54, 56, [59] [60] [61] [62] [63] [64]. Therefore, the use of only one set of tracers will likely not maximize the resolution of all fluxes in SLEs, particularly when a large-scale metabolic model is employed [28, 65]. 

According to previous studies [53, 58, 66, 67], there are several advantages of using PLEs for 13 C-MFA compared with an SLE-based approach. In general, the data from each LE are integrated to achieve an improved flux resolution, primarily due to the synergy of the complementary information used for fitting to the single metabolic model [48, 53, 58, 66]. Indeed, the latest applications of the COMPLETE (short for COMplementary Parallel Labeling Experiments TEchnique [58] ) MFA approach employing all six singly labeled glucose tracers to evaluate metabolic fluxes resulted in the most accurate and precise flux parameters obtained thus far for wild-type E. coli as well as for some metabolically engineered strains of the bacterium [58, 67]. However, for laboratories lacking in-house experience, one crucial factor in the implementation of the PLE approach is the availability of a free, ready-to-use software package allowing the successful manipulation of the complex data obtained in PLEs, which is necessary for comprehensive flux analysis. 

In the present study, the open-source software Open-FLUX [33], which uses an elementary metabolic unit (EMU) decomposition-based algorithm to generate an isotopomer balance model [26] and was initially developed for SLE analysis, has been extended for the computation of PLE data (see, Additional file 1: SF-1. 3. The methodology of PLE data implementation is rather clear, and one of the possible algorithms has been earlier schematically described in [53]. The expertized investigators have already adjusted their home-made 13 C-MFA software by PLEs-mediated data (see, [66] for review). Currently, additional MATLAB-based scripts have been appeared on the OpenFLUX homepage (http://openflux. sourceforge. net) that demonstrated to users how the data of two labeling experiments conducting in parallel could be implemented in the already existing software. The presented open-source OpenFLUX2 provides a friendly software environment that facilitates beginners and existing OpenFLUX users to manipulate with SLE-and PLE-based data, for experimental design, determination of flux parameters, and for broaden evaluation of flux statistics. 

Using OpenFLUX2, direct in silico simulation confirmed that the flux resolution was improved when 13 C-MFA was provided with PLE data that were fitted to and integrated with the common metabolic model as compared with the individual analysis of each LE. Additionally, the best flux resolution was achieved in the analysis of PLE results when the choice of tracer for each provided LE was based on a computed experimental design targeted to minimize the approximated variances of several fluxes from the different parts of the assumed metabolic network. The statistical methods of analysis of the obtained experimental and simulated data, followed by a goodness-offit test of the adequacy of the applied metabolic model, have been extended in OpenFLUX2, including the statistical conclusions concerning the feasibility of the obtained flux parameters in the final report and the flux confidence intervals estimated at the desired significance level. In turn, the flux confidence intervals could be computed in OpenFLUX(2) using different methods, but up today the most dependable and precise approach is a fine-tunable Monte Carlo-based determination of the flux variances, which are dependent on the randomly corrupted measured data [33, 39] and that has been modified in OpenFLUX2 due to implementation of a convergence control and visualization of computation results. 

Following the original position of the OpenFLUX authors [33], OpenFLUX2, which is an extended version of the already available software, has been developed as open-source software. We hope that OpenFLUX2 will be useful to research groups applying 13 C-MFA, particularly for beginners not yet experienced in fluxomics analyses. Additionally, the availability of the OpenFLUX2 code could promote further improvement of the software based on the experiences of different researchers. OpenFLUX2 can be downloaded from SourceForge (http://sourceforge. net/ projects/openflux2).
section_name: Background

section_name: Results

section_content: OpenFLUX2 was developed as an extension of the OpenFLUX software. New calculation facilities were added, mostly as extensions of the initial options, without dramatic changes in the parent content. Moreover, the initial forms of the model and experimental data setup, together with results representation, were maintained as much as possible during the development of OpenFLUX2 to facilitate the transition from one version to the other. In the present study, the procedures that were developed previously in OpenFLUX and retained in OpenFLUX2 without essential modifications are indicated as "OpenFLUX(2)", and only the added/modified elements are indicated as being implemented in OpenFLUX2. 

To clarify the essence of the modifications implemented at the stage of OpenFLUX2 software development, the following items are schematically described in Additional file 1: (SF-1. 1. ) the assignment of free fluxes, followed by (SF-1. 2. ) flux variability analysis; (SF-1. 3. ) the calculation of optimized fluxes through iterative fitting; (SF-1. 4. ) a goodness-of-fit analysis of the adequacy of the metabolic model; (SF-1. 5. ) local linearized statistical approximations; and non-linear-search of the optimal flux confidence intervals (SF-1. 6. ), where, in particular, the introducing a convenient concept, "the normalized flux precision" function, is described, as well; (SF-1. 7. )a fine tunable and convergence-controlled Monte Carlo-based approaches for precise determination of the optimized flux confidence intervals according to "discarding" strategy at the predetermined confidence level, significantly modified at the stage of implementation in OpenFLUX2. The main aim of this description is to demonstrate that the individual procedures are essential interconnected parts of a unique solution to a complex optimization problem, where the statistical significance of the calculated modelbased parameters must be verified via the comprehensive goodness-of-fit of the model's adequacy. As an auxiliary aim of this part, it is a rather short, but slightly (in comparison with the excellent review [4] ) mathematicallyenriched introduction in the 13 C-MFA background that could be helpful, especially for beginners, to repair their knowledge by essential parts of linear algebra and statistics. 

The workflow of OpenFLUX software consists of two components: (i) the automated set-up of stoichiometric and isotopomer balance models from user-supplied data, performed by Java PARSER, and (ii) the application of the generated models to the 13 C-MFA of SLEs for flux parameter estimation and sensitivity analysis, performed by a set of specially developed MATLAB functions. The modifications introduced in OpenFLUX2 can be divided into two types (Figure 1 ). First, a set of new options related to statistical analysis of the flux evaluated from SLE data and experimental design facilities were added. Second, the new software was extended to perform metabolic flux analysis of PLE-based data. Here, the PLE concept is considered as determined in [66] ; PLE-based design means that several tracer experiments that are started from the same seed culture to minimize biological variability, and must be performed independently under identical conditions on the same substrate(s), using different tracer(s). 

Because an 13 C-MFA PLE-based approach requires the simultaneous fitting of several datasets obtained from independent LEs to a single model, there is no major difference in the spread-sheet model set up and the consequent automated generation of stoichiometric and isotopomer balance models by the Java PARSER for PLEs and SLEs in OpenFLUX2. The set of substrates is fixed during the model generation step, and individual substrate tracer configurations are then defined by the user for each LE constituting the PLE together with the corresponding measured data (Figure 1 ). The option to use either a single label or a labeling mixture for each substrate in the PLE is provided by OpenFLUX2, as was previously provided in OpenFLUX for SLEs. Thus, all of the introduced modifications were finally concentrated in the MATLAB-based portion of the computational algorithm. 

Comprehensive flux analysis of a Corynebacterium glutamicum model, as an example, using OpenFLUX2 software The C. glutamicum model

The specific features of the developed OpenFLUX2 software were illustrated via the in silico simulation of a set of SLEs and PLEs followed by their comprehensive 13 C-MFA. The Corynebacterium glutamicum metabolic model initially proposed for an L-lysine-overproducing C. glutamicum strain by Becker et al. [68] and later used to demonstrate the possible applications of OpenFLUX software [33] was chosen (with minor changes; see Methods) as the object of the computational simulation and 13 C-MFA. This model, schematically presented in Figure 2, accounted for the main central metabolic pathway, the biosynthetic pathways for several excreted products, including lysine as the main product, and anabolic demand. CO 2 -mediated carbon transfer was accounted for using expression reactions, accompanied by CO 2 production/consumption, in an explicit manner. The bi-directional reactions were represented as nonnegative forward and reverse fluxes. Finally, the metabolic model was composed of 51 unknown fluxes balancing 36 metabolites, resulting in the stoichiometric matrix S, dim(S) = (36 × 51), with a residual 15 degrees of freedom. The degrees of freedom were associated with the following free fluxes: experimentally determined effluxes, including biomass biosynthesis (1 flux), effluxes of secreted products (5 fluxes), and the glucose uptake rate (1 flux); reverse branches of five bi-directional reactions; and free irreversible reactions associated with branch points. An isotopomer model was generated automatically under the assumption that the mass isotopomer distributions of the following metabolites were measured: alanine, valine, threonine, aspartate, glutamate, serine, phenylalanine, glycine, tyrosine, and trehalose. Application of the EMU approach offered by OpenFLUX(2) led to a major simplification of the isotopomer model (see Methods).
section_name: Key features of OpenFLUX2 software

section_content: The full set of fluxes necessary for the simulation of in silico SLEs and PLEs for the defined model was unfortunately not directly available from previously published materials [33, 68]. Thus, the evaluation of fluxes from published experimental data was repeated in the present study. To this end, a new constrained NLLSP was designed based on a modified metabolic model, with assigned free fluxes, previously measured [68] effluxes with their variances, and MIDs with variances assumed to be 0. 15 mol% (see Methods). Solutions of the designed NLLSP were obtained by OpenFLUX2 using a gradient-based local optimization method with the assistance of MATLAB's FMINCON function (see Additional file 1: SF-1. 3., for details). To determine the global minimum, 100 independent iterative trials were performed with an arbitrary set of initial free fluxes, θ k ∈ ℜ p, k = 1, 2, …, 100, where ℜ p indicates a feasible constrained domain in the free flux variation space (see Additional file 1: SF-1. 2. ). These trials resulted in a set of flux estimates, uð ⌢ θ k Þ; k ¼ 1; 2; …; 100, corresponding to the minimum values of the objective function SSR SLE f À Á k in Eq. (S -1. 3. 7), which were reached in each k-th trial. According to the optimization report, all of the provided trials were terminated successfully due to the achievement of the given termination tolerance (ΤΤ = 1 × 10 -4 ) within the constraint violation tolerance. The χ 2 -based goodnessof-fit confirmed the model adequacy (see Additional file 1: SF-1. 4., for details) for 82 (at a significance level of α = 0. 05) and 34 (at α = 0. 32) of the 100 obtained values of

where W is the number of independent measurements, and p is the number of estimated free fluxes ((Wp = 21) in this case). All essential information concerning the obtained solutions is presented in The flux identifiability analysis, which was based on model linearization and on the computation of Null J f ð ⌢ θÞ in Eq. (S -1. 5. 9) (see Additional file 1: SF-1. 5. ), was provided for the obtained statistically available solution as described by Yang et al. [69]. This analysis resulted in the conclusion that only a unique set of flux parameters for the global minimum of the constrained NLLSP could be computed numerically; because the calculated null space matrix was empty. The goodness-of-fit analysis was finalized by confirming the Ν(0, 1) distribution of the individual variance-weighted residuals in Ξð ⌢ θÞ (see Additional file 1: SF-1. 4. for details). The obtained statistically acceptable solution of the constrained NLLSP, uð ⌢ θÞ, primarily coincided with the previously published [33, 68] values of fluxes in the range of the earlier evaluated flux confidence intervals (Additional file 2: Figure SF-2. 2). Thus, all flux parameters, including previously unavailable, were evaluated and assigned as the true values for the assumed metabolic model, u(θ true ). 

Other χ 2 -statistically acceptable trials resulted in larger values of Ξ(θ), which corresponded to the points of Figure 2 Assumed metabolic model of an L-lysine-producing Corynebacterium glutamicum strain grown on glucose. Single (solid) and double (solid and dotted) lines with arrows correspond to the assumed irreversible reactions and to bi-directional (forward and reverse) reactions, respectively. The letter "θ" indicates fluxes assigned as free, and the letter "v" indicates dependent fluxes. The free and dependent fluxes are numbered continuously throughout the applied model. Amino acid biosynthetic pathways are shown schematically as the drain of precursors of the corresponding amino acid. Carbon sources are indicated with ovals. The following abbreviations are used for metabolites: 3PG -3-phosphoglycerate; ACCOAacetyl-CoA; ACETALacetaldehyde; ACETAL_EXextracellular acetaldehyde; AKG -α-ketoglutarate; AKG_EXextracellular α-ketoglutarate; CITlumped pool of citrate and isocitrate; CO2 -carbon dioxide; CO2_EXextracellular carbon dioxide; DAPdiaminopimelate; DHAPdihydroxyacetone phosphate; E4P -D-erythrose 4-phosphate; F16BPfructose 1,6-bisphosphate; F6P -D-fructose 6-phosphate; G3Pglyceraldehyde 3-phosphate; GLC_EXextracellular glucose; GLC6P -D-glucose 6-phosphate; GLY_EXextracellular glycine; LAC -lactate; LAC_EXextracellular lactate; LYS_ EXextracellular lysine; MALmalate; OAAoxaloacetate; P5Ppentose-5-phosphates (lumped pool of D-ribulose 5-phosphate, D-ribose 5-phosphate, D-xylulose 5-phosphate); PEPphosphoenolpyruvate; PYRpyruvate; S7P -D-sedoheptulose 7-phosphate; SUCsuccinate; TREMtrehalose constituent (D-glucose 6-phosphate and UDP-glucose); TRE_EXextracellular trehalose. Amino acids are expressed using the standard three-letter code. To determine the global minimum of the NLLSP, 100 independent iterative trials were provided with an arbitrary set of initial free fluxes for all represented in silico LEs. All iterative trials were terminated successfully for each independent NLLSP. More detailed characteristics for all performed in silico experiments are listed in Additional file 3. 

(a)

In the case of the SLE, the vector of the measured data consists of both measured MID values and measured effluxes. In the case of the PLE, the vector of the measured data consists of the measured MIDs and effluxes that were available for each SLE. 

(b)

The measured effluxes are included in the set of free fluxes, constrained according to Eq. (S -1. 2. 2) with experimentally determined parameters q mea eff ; σ mea eff. 

(c)

The values of Ξ k that were obtained during different iterative trials were assumed to be equal to the Ξ min if the values were different from the Ξ min by a value less than or equal to the TT level, i. e., Ξ k -Ξ min ≤ TT. 

(d)

The analysis was performed for the NLLSP solution, which provided the minimal achieved value of Ξ, i. e., Ξ min. 

(e)

The individual variance-weighted residuals were recognized to be not distributed as N(0, 1) if both the Kolmogorov-Smirnov test, provided by MATLAB, rejected this hypothesis and the individual residual plot revealed unpredictable high values for the individual residual/group of the individual residuals or, in contrast, unpredictable small values for the primary part of the individual residuals. Thus, additional solutions of the constrained NLLSP may be obtained in independent trials starting from the random initial values of the free fluxes, resulting in the premature termination of the computer search for the global minimum due to the discovery of a termination tolerance (TT). Indeed, repeated numerical solutions of the same constrained NLLSP under conditions with a decreasing TT of the objective function value sloped from the default value of 1 × 10 -4 up to the more precise value of ΤΤ = 1 × 10 -6 (see Additional file 1: SF-1. 3., for details), resulting in an increased number of trials that reached the global minimum (77 of 100, in contrast to 26 in the previous calculations; see Table 1 Thus, a unique optimal set of fluxes for the statistically adequate metabolic model was obtained: uð ⌢ θÞ≡u θ true ð Þ. Then, new GC-MS-based "experimental data," i. e., new MIDs, could be generated through direct in silico simulation describing the propagation of 13 C atoms from different tracers through a metabolic network with known flux parameters (see Methods).
section_name: Circumstantiation of flux parameters for the computer simulations

section_content: Initially, the simulation was provided for 99% of [1-13 C]glucose as a tracer. The new set of MIDs was generated as described in Methods, excluding the last step of data corruption. To confirm that flux estimations could be unambiguously inferred from the obtained non-noised data, these simulations were used as the experimental data, with assumed MID variances equal to 0. 4 mol%, and the assignment of all effluxes as variable free fluxes constrained in the range of the 95% confidence intervals with the standard deviations determined in [68] (see Additional file 3: Table SF-3. 3 ). The solution to the corresponding constrained NLLSP was obtained using OpenFLUX2 software according to the standard procedure described above with details presented in Table 1 (Exp_2. 1). In total, 98 values of Ξ i, which corresponded to solutions from the total obtained set of uð ⌢ θ i Þ; i ¼ 1; 2; …; 100, were smaller than the upper, and even the lower, critical threshold values, χ 2 0:975 W -p ð Þand χ 2 0:025 W -p ð Þ, respectively, at the 95% confidence level and with (Wp) = (38 + 7) -(8 + 7) = 30 degrees of freedom. Moreover, the group of trials (6 solutions at ΤΤ = 1 × 10 -4 and 51 solutions in the case of TT decreased to 1 × 10 -6 ; see Table 1 (Exp_2. 1; Exp_2. 2)) had a minimal value of 1. 5 ⋅ 10 -3, which was significantly less than the minimal threshold. Such a questionably small value for the objective function generally indicates possible overfitting of the applied model. However, in our case, the cause stems from the solution of the inverse task without corrupting the "experimental" data generated at the stage of direct simulation. The same cause resulted in the negative evaluation of the individual weighted residuals in Ξ i according to the normality test, and the null hypothesis (concerning Ν(0, 1) distribution of residuals) was rejected. The flux estimates that corresponded to the solutions of the group with the minimal Ξ i that was reached were assigned as reference The next stage of the simulation was the generation of "experimental data" that were closed to results that could be obtained in real experiments. The values of the previously simulated MIDs and the measured effluxes were corrupted with the Ν(0, σ mea ) normally distributed random errors using the Statistics Toolbox of MATLAB. Five

2; …; L were used for calculation of mean, and unbiased estimator of the variance, respectively:

for each from w MID MIDs and w eff effluxes, (w = w MID + w eff ), and to assign the SSR SLE f objective function that finally determines the constrained NLLSP (S -1. 3. 8). Statistically acceptable solution (corresponding to the value 9. 65 of Ξ(θ) function (for the group of 27 from 100 performed trials) that was smaller even the lower threshold, χ 2 0:025 30 ð Þ ¼ 16:79 À Á ) at the 95% confidence level, with Ν(0, 1) distributed weighted residuals, was found using OpenFLUX2 (Table 1 (Exp_3. 6 )).
section_name: C-MFA for in silico SLEs with [1-C]-glucose as a tracer

section_content: At the final stage of the 13 C-MFA performed for the considered above constrained NLLSP, accurate 68% and 95% confidence intervals (CI MC 0:68 and CI MC 0:95, respectively) were computed for the optimized fluxes initially by fine-tunable and convergence-controlled Monte Carlo-based approach implemented in OpenFLUX2 for this purpose (see Additional file 1: SF-1. 7. for details). According to [39], Monte Carlo-mediated analysis of flux statistics was carried out on the basis of a discrete approximation of optimized flux estimation distributions obtained in L multi-trials when the experimental data for each trial, comprising measured MIDs and effluxes, were artificially generated by corrupting of real initial data with normally distributed random errors. So, consecutive providing of L optimization trials finally resulted in L-set of optimized flux estimations:

where ⌢ Usually, the Monte Carlo search of CI MC γ u i ð Þ bounds is based on an assumption that their estimated values have to converge in case of significant increasing of the total number of trials. So, the number of trials, L, is one of the most important parameters of Monte Carlo-based procedure, and it has to be optimally chosen (L = L MAX ) for CI MC γ u i ð Þ precise determination in a reasonable computation time. The special procedure was implemented in OpenFLUX2 for a control of an essential number of optimization trials that was performed during target flux CI MC γ u i ð Þ bound estimation. This control finalizing in determination of all CI MC γ u i ð Þ, could be realized ON-LINE, with a help of specially preset control parameters (the second group of parameters implemented for fine tuning of the Monte Carlo-based search procedure), or according to direct user's decision based on visualization of estimated bound plots in dependence on the current L, and flux estimation histograms that could be presented after the predetermined L MAX trials were performed. 

The developed Monte Carlo-based approach was used for determination of CI MC γ u i ð Þ for the solution of the constrained NLLSP described above (Table 1 (Exp_3. 6)). The optimized flux parameter distributions were obtained for L MAX = 500 trials using the default set of tunable control parameters that facilitate search of global minimum in each trial of the applied "multi runs per trial" fitting strategy (N AS =3, K NR =50, ΤΤ ¼ 1 Â 10 -4 ; ε ¼ 1 Â 10 -4, see Additional file 1: SF-1. 7., for details). Then, for all 51 fluxes of assumed metabolic network the sliding control was performed for upper and lower bounds of CI MC γ u i ð Þ evaluated initially according to "discarding" strategy (i. e., bounds of CI MC-1 u i ð Þ were closed to zero. So, these bounds, manifesting the obvious convergence, had to be restricted at the level of absolute (not relative) value of spreading (defined in Eq. (S -1. 7. 6)), e. g., Δ MC-1 LB γ ð Þ≤2. In any case, all bounds of CI MC-1 γ u i ð Þ could be obviously determined after the Monte Carlo-based optimized flux estimation had been obtained after L MAX =500 performed trials in case of "discarding" strategy was used. In case of "mean-varianced" strategy was applied, the main part of CI MC-2 0:95 6A ). As could be seen from the (Figure 6C ), the reason of this rather slowspeed convergence of several flux bounds was appearance of a small quantity of flux estimations significantly exceeded the ordinary values: "mean-varianced" strategy could not discard these outstanding values and had to accumulate a lot of usual estimations for significant decreasing the totally calculated "mean". 

Summing up, it could be concluded that obtaining a proper approximation of the optimized flux parameters was the most important part of the Monte Carlo based search of CI MC γ u i ð Þ: even rather small quantity of significantly differed "outstanding" values in the flux estimations distribution that were obtained when the global minimum was not achieved in the fitting procedure, could significantly increase the width of the evaluated CI MC γ u i ð Þ. In the tested cases, the convergence of CI MC γ u i ð Þ bounds was achieved faster (at the smaller number of the performed trials) if the "discarding", but not "mean-varianced" strategy of the bound determination was used. So, the "discarding" strategy manifested significantly higher resistance for these occasional "outstanding" estimations. Moreover, this strategy was preferable to demonstrate the asymmetric character of optimized flux estimation distribution resulting in the asymmetric locations of the upper and lower bounds presented simultaneously for CI MC-1 0:95 u i ð Þ, and CI MC-1 0:68

in comparison with symmetric locations of their bounds obtained in case of "mean-varianced" strategy was used (Figure 4 ). 

The later feature of the CI MC-1

well with the parameters of CI n-lin γ u i ð Þ independently obtained by non-linear search developed by Antoniewicz et al. in [28] and implemented in OpenFLUX (see Additional file 1: SF-1. 6. for details). All obtained data were presented in the Additional file 3: Table SF-3. 6. As could be seen, the most part of CI n-lin γ u i ð Þ parameters coincided rather well with Monte Carlo-based results, especially with CI MC-1 γ u i ð Þ bounds (this was true from 39 of 51 fluxes). Nevertheless, in several cases evaluations of the CI γ (u i ) bounds given by Monte Carlo and non-linear approaches differed, and some times significantly (e. g., UB MC-1 0:95 ≤UB n-lin 0:95 and, on the contrary, UB MC-1 0:68 ≥UB n-lin 0:68 for the θ 14, and θ 22 fluxes; both determined upper bound values for the θ 31 -flux confidence interval, UB n-lin 0:95 and UB n-lin 0:68, was lower in case on non-linear computing). Absolutely statistically incorrect result was obtained for 9 fluxes (e. g., θ 32, θ 33, v 2, v 5 and etc. ): estimated upper bound of CI n-lin 0:68 u i ð Þ had higher values than their upper bound of CI n-lin 0:95 u i ð Þ. It seems that these "mistakes" appeared as a result of low accuracy of numerous calculations performed according to nonlinear search-based algorithm in the used software. It finally resulted in termination of computation even when necessary optimality conditions were not satisfied and the real global minimum was not reached in the optimization procedure. The proposed modifications targeted to improvement of the calculation efficiency are at the final stages of testing and implementation in new release of OpenFLUX2 software (see, Additional file 1: SF-1. 6., for details). Up today, the current version of OpenFLUX2 contains the initial variant of non-linear algorithm of flux confidence intervals search. Keeping in mind incorrect results that could be computed now for CI n-lin γ of some fluxes, that are very difficult to recognize without information for comparison, but, on the other hand, very high speed of all flux CI n-lin γ u i ð Þ computation (about one hour for estimation flux statistics for SLE using computers described in Methods), it could be highly recommended to use the current algorithm of CI n-lin γ u i ð Þ non-linear search mainly for quick preliminary evaluation. The accurate determination of CI γ (u i ) could be performed, e. g. as CI MC-1

according to the fine-tunable Monte Carlo based approach with automatic and/or visual control of all bounds convergence.
section_name: Monte Carlo-based and non-linear approaches for determination of flux confidence intervals

section_content: To compare the achieved efficiency of the u i flux resolution, it is convenient to use the normalized flux precision function (see, Additional file 1: SF-1. 6. ) with β as the scaling parameter:

Generally, in this function, u i ⌢ θ is the best available estimation of the unknown true value that can be computed from the SLE-or PLE-based 13 C-MFA. In the case of computer simulations, the "true" flux parameters are known: namely u i ⌢ θ ¼ u i ð Þ true were used for the calculation of η γ (u i ; β) values in the present study. In turn, the superscript employed for the η γ (u i ; β) and CI γ (u i ) functions, e. g., n-lin and MC -1, respectively, indicates the method of estimation of the flux confidence interval. So long as the Monte Carlo-based approach with application of "discarding" strategy was used in all examples of the present study as the main method for flux confidence interval determination, we would not specially indicate below the way of determination of flux confidence intervals using the superscript " MC-1 ". The variable scaling parameter β was set to 0. 1 in the present study, and the dependence of the η γ (u i ) function on β as the parameter is not directly indicated in the corresponding equations below for brevity. According to its definition, the η γ (u i ) function at each fixed β parameter is close to "1" for precisely estimated u i fluxes (with narrow confidence intervals) and is close to "0" for a poorly determined u i. 

The computed values of η γ (u i ) for the obtained solution of the NLLSP were flux specific and varied from zero to almost unity, with a sum of Σ η 0:95

6 for details). The flux specificity of the η γ (u i ) function for the determined metabolic model was subsequently shown to be primarily dependent on the applied labeled tracer (see Figures 7 and 8 ). Indeed, the high resolution of one set of fluxes (e. g., v 3 : GLC6P → FR F6P; θ 4 : F6P → R GLC6P; The summarized values of the η -function for all fluxes, (Σ η(0. 95) ([1 -13 C] SLE )) k, which were calculated for each flux from k = 1, 2, …, 5 SLEs, varied in a rather narrow range (between 37. 9 and 38. 5). Considering the total number of fluxes (equal to 51) for the assumed metabolic model and the detected measurement-dependent sensitivity of the η -function, if the difference between the Σ η(0. 95) values is no greater than 1, then the LEs can be considered to be provided with an essentially equivalent flux precision. Certainly, the value of Σ η(0. 95) could be considered only as a general, conditional measure of flux resolutions in the investigated network: some metabolic branches could be resolved better, and other branchesworse in different experiments with, perhaps, equal values of summing normalized flux precision functions. So, only the value of the η γ (u i ) -function could be considered as the absolute measure of the u i -flux resolution efficiency that could be compared with other LEs performed with the same metabolic model.
section_name: Normalized flux precision function as a measure of flux resolution efficiency

section_content: In the present study, the results of the 13 C-MFA included the set of estimated fluxes, the goodness-of-fit of the model's adequacy, and the confidence intervals of fluxes, in accord with the recently recommended publishing guidelines [70]. As shown in Additional file 1: SF-1. 4., the comprehensive goodness-of-fit analysis had to consist of the χ 2 -mediated testing of the Ξ(θ) objective function at the point of convergence, ⌢ θ, and confirmation that the individual weighted residuals used as the summands of this function were Ν 0; 1 ð Þ distributed. The solution of the constrained NLLSP was considered statistically acceptable only if all of the tests were successfully passed. 

In this report, an example is presented in which the mistaken flux parameters could be assumed when the statistical analysis of the obtained solution was partially provided but to an insufficient extent. Let us analyze the first from five earlier described in silico experiments with 99% [1-13 C]-glucose as a tracer (see Table 1
section_name: The necessity of a comprehensive statistical analysis of the NLLSP solution

section_content: The "contribution matrix" (see Additional file 1:

, was computed at the true point, u(θ true ), as an important component of the flux statistics (Additional file 3: Table SF-3. 12). It is known [28] that the (CM) ij elements of this matrix indicate the relative importance of the variance of the j-th measurement to the local variance of the i-th flux. As can be observed in Additional file 3: Table SF-3. 12, the variances of the "serine" MIDs demonstrated the high importance of the flux resolution among MS measurements; all matrix columns corresponding to the SER mass isotopomers showed rather high sums of their elements. The variances of these "serine" ("SER") MIDs significantly influenced the resolution of several fluxes (θ 8 ; θ 22 ; v 23 ; θ 31 ) for the following reactions: GLC6P → F P5P + CO2; PYR + CO2 → F OAA; OAA → F PYR + CO2; and CO2_EX → R CO2, respectively. The new set of "experimental data" was generated in the following fashion: the measured effluxes and all MIDs, except for "SER" MIDs, were considered as in the previously analyzed example. The "SER" MIDs were modified as in the case of "poor" resolution of the SER-390 (m + 0) MID, which exhibited an unknown byproduct that increased the value of the corresponding SER peak to +4. 5%. Due to the necessity of normalizing all SER-isotopomers, the applied modification resulted in a proportional decrease in the other SER MID portions; therefore, the SER-390 MIDs were modified from (m + 0)/ (m + 1)/(m + 2)/(m + 3) = 0. 443/0. 357/0. 140/0. 042 (in the previously described example) to a ratio of 0. 463/0. 344/ 0. 135/0. 040. 

This "mistakenly" modified set of "experimental data" was used to generate the newly constrained NLLSP. A unique solution, corresponding to the minimal values of the Ξ = Ξ(θ) objective function with the empty null space of the J f (θ) matrix evaluated at the new point of convergence, θ ¼ ⌢ θ new, was found. Moreover, the value of Ξð ⌢ θ new Þ was in the χ 2 -statistically acceptable interval (see Table 1 (Exp_4. 1)). Simultaneously, the parameters obtained for several fluxes did not coincide with the true flux values even in the range of their confidence intervals. Among them there were fluxes θ 8 and θ 31 (Figure 9 (I -B )) which exhibited variances that were essentially determined by SER MID measurements. The fact that the obtained flux parameters were incorrect could be established only at the stage of normality distribution testing of the weighted residuals; three variance-weighted residuals corresponding to "SER-390" MIDs and one corresponding to "alanine (ALA)-260" were shown to be located outside of the 95% confidence interval (Figure 9 (II-B)). According to the general recommendation presented in SF-1. 4, improvement of the solution statistics could be achieved due to the excision of the "prominent" residuals from the expression for the objective function followed by the repeated solution of the constrained NLLSP with the partially truncated Ξ function. Two alternative NLLSPs were generated where ALA or SER (See figure on previous page. ) Figure 7 The values of the normalized flux precision function in the different parts of the assumed metabolic network depend significantly on the applied tracer(s). The results of flux estimation: (A) SLE using 100% [1-13 C]-glucose as the tracer; and from PLEs consisting of: (B) 5 LEs using different mixtures of [U- 13 C]/[U- 12 C]-glucose in each LE, with 20%, 35%, 50%, 65%, or 80% [U- 13 MIDs were excluded from the "experimental data" set. Parameters obtained after solution of the NLLSP with ignored "ALA-260" MIDs did not coincide with the true flux values in the range of their confidence intervals (Figure 9 (I-C )). Moreover several variance-weighted residuals still located outside of the 95% confidence interval including that of "SER-390" (Figure 7 (II-C)). In contrast, the proper statistically acceptable solution was obtained for NLLSP with ignored "SER-390" MIDs, which coincided with true fluxes in the range of the determined flux confidence intervals (Figure 9 (I-D ) and (II-D)). Thus, "incorrect flux parameters" could be assumed if test for normality of variance-weighted individual residual distribution is ignored. Certainly, higher negative influence on the proper flux estimation could be assumed if χ 2 -statistically unacceptable solution is obtained. Analysis of variance-weighted residuals plot followed by step by step exclusion substances with prominent residuals, as described in the example, could help to identify "incorrect measurements" and, perhaps, to obtain statistically acceptable solution with correctly estimated parameters. 

It is obvious that measured MIDs and effluxes' parameters are completely separate categories of experimental data that could not be trivially compared. It is well established [3] that labeling experiments are performed to resolve internal fluxes, even parallel and cycle pathways, and reversible reactions, which cannot be resolved on the basis of measured effluxes. Nevertheless, the variances of the measured effluxes could provide the most significant influence on the resolution of some fluxes, as could be seen from the values of the corresponding (CM) ij elements in the Additional file 3: Table SF-3. 12. 

Thus, generally, it is desirable to execute the efflux measurements with the highest possible accuracy to decrease the variances and totally improve the flux resolution. One of the interesting step in this direction has been done recently when the authors tried to increase the measurement accuracy for an efflux corresponding to quantifying biomass composition due to exploiting of the high-precision GC-MS technique [71]. Unfortunately, in many cases, the efflux measurements as the stage of labeling experiment have received much less attention than the more glamorous stages of the subsequent highly-precised MS-based measurements.
section_name: (Exp_3.1)).

section_content: A uniformly 13 C-labeled isotopomer of glucose, [U- 13 C]glucose, is often used as a tracer in 13 C-MFA. In the present study, new sets of "experimental data" were generated for the same metabolic model when different relative amounts of [U- 13 C]-glucose (20%, 35%, 50%, 65%, and 80%) mixed with non-labeled ([U- 12 C]) glucose were used as the sole carbon source (Additional file 3: Table SF

Statistically acceptable solutions were obtained for constrained NLLSPs for five independent SLEs and for a PLE using all of the generated "experimental data" for fitting to the single metabolic model (see Table 1 SF-3. 8 ). Again, a set of tracer-specific fluxes that possessed rather high values of the η -function could be detected (at least more than in the case of exploiting [1-13 C]-glucose as a tracer), e. g., θ 22 : PYR + CO2 → F OAA; v 23 : OAA → F PYR + CO2; θ 31 : CO2_EX → R CO2. The PLE-based 13 C-MFA of experiments using [U- 13/12 C]glucose as the carbon source actually improved the resolution of all fluxes estimated in the corresponding SLEs, Σ η(0. 95) ([U - 13/12 C] PLE ) = 42. 4. Moreover, the tracerspecific behavior described in relation to SLEs was reproduced in the PLE. The colored scheme presented in Figure 7 and the diagram in Figure 8 correspond to the calculated values of the η -function for all fluxes obtained in different experiments, illustrating this fact and demonstrating that optimization of the experimental design is necessary to increase the precision of the targeted fluxes [19, 49, 54, 62].
section_name: Simulations of LEs with [U-13 C]-glucose as a tracer

section_content: To achieve the best flux resolution in an LE employing a mixture [U- 13 C]-, [1-13 C]-, and [U- 12 C]-glucose isotopomers as the carbon source, optimization of the experimental design could be performed according to the method proposed by Möllney et al. [19]. In their study, the comparison of different designs was based on an evaluation of D-factor values [72], i. e., the squared volumes of the flux confidence ellipsoids for a given confidence level. According to [19] and using the notations introduced in the present study, the squared volume of the p-dimensional confidence ellipsoid for free fluxes, which was evaluated at the point of convergence, θ ¼ ⌢ θ, could be obtained (up to a constant that does not depend on designed parameters) using the following expression:

where Σ⌢ θ is a covariance matrix of the free fluxes estimated according to linearized statistics (see Additional file 1: SF-1. 5. ) and to the approximation of this matrix by the inverse Hessian matrix presented in Eq. (A2 -10) Additional file 1: SF-1 Appendix 2. Specifically, the last equality from Eq. ( 5 ) was used for computation of the D-factor. According to a previous suggestion [19], the use of the value of the ffiffiffi ffi D 2p p parameter, which can be interpreted as an average length of the confidence interval of estimated fluxes, is more suitable. As noted previously, the confidence ellipsoid volume estimation provided by Eq. ( 5 ) holds true within some vicinity of the predefined point θ ¼ ⌢ θ and can change due to a shift in that point [19]. Because the true values of the free fluxes were known for our artificial metabolic model, the optimization of the experimental design, which was dependent on the tested 13 C-labeled tracers, was significantly simplified and was based on the computation of the determinant of the inverse Hessian matrix evaluated at the point θ = θ true. In practice, when the true flux values are not known a priori, they are assumed, for example, from data in the literature, FBA or even from 13 C-MFA performed on the basis of a rather inexpensive label [37, 63]. In some cases, a second round of experimental design optimization is necessary if the proposed flux distribution is far from the feasible point of convergence obtained in a planned labeling experiment [73, 74]. In the present study, all of the new experimental designs were compared with the reference design (SLE with 99% [1-13 C]-glucose as the tracer) using the value of the ffiffiffi ffi D 2p p parameter. addition to the "general optimization" achieved due to the minimization of the relative value of the average confidence interval length:

"partial optimization" could be obtained as a result of minimization of the standard deviation of the targeted θ i free flux, which, in turn, was approximated as the square root of the corresponding diagonal element of the Σ ⌢ θ matrix:

The results of the corresponding computations performed using the special subprograms implemented in OpenFLUX2 software are presented in Additional file 2: Figure SF-2. 6. According to the provided calculations, "general optimization" required the use of a mixture of [1 -13 C] 78 % /[U -13 C] 22 % glucose isotopomers in the SLE (Additional file 2: Figure SF-2. 6 (A)). Interestingly, these D-factor-mediated optimized tracer compositions were extremely close to the [1 -13 C] 80 % /[U -13 C] 20 % -labeled mixture of glucose that has been used without any calculations to achieve a rather high flux resolution in other experimental systems (e. g., Escherichia coli-based systems [48], in particular). In contrast, the same D-criterionbased approach resulted in another optimal mixture of the same glucose isotopomers ([1 -13 C] 48 % /[U -13 C] 40 % / [U -12 C] 12 % ) when some modifications differed in the metabolic model of the L-lysine-producing C. glutamicum strain applied in the present study, and plans were made to obtain another set of measurements [19]. 

According to the provided computations, "partial optimization," i. e., the optimal resolution of different targeted free fluxes, had to be achieved in the SLEs with the distinguished mixtures of the same glucose isotopomers (e. g., see Additional file 2: Figure SF-2. 6 (B -F)). Several of the proposed optimal designs (e. g., see Figure 10 (I) ) were realized using computer simulations, and the set of experimental data was generated, followed by a search of the solution for the corresponding constrained NLLSP. One of the simulated SLEs was performed using the "generally optimized" mixed tracers, five others employed the "partially optimized" mixtures, and five experiments were performed with the randomly chosen mixture of [1 -13 C], [U -13 C] and [U -12 C] glucose isotopomers ("randomly mixed") as the tracers. As shown in the presented results (Additional file 3: Table SF-3. 9), the SLE with "generally optimized" conditions resulted in Σ η(0. 95) ([1

Þ ∈ 36:0; 38:7 ½, and with "randomly mixed" tracers -Σ η 0:95

. It has to be especially noted, that SLEs with the higher calculated values than Σ η(0. 95) = 38. 4 was detected here, and among described earlier (see, final Additional file 3: Table SF-3. 11), and so "general" and/or "partial" optimizations does not result in maximization of the sum of normalized precision functions for all fluxes, Σ η(0. 95). Simultaneously, the free fluxes, θ i, that were used as the targets for the partial optimization generally showed the best resolutions (i. e., the narrowest confidence intervals, CI 0. 95 (θ i ), (Figure 10 (II)) and the highest values of the η 0. 95 (θ i ) function (Figure 10 (III)) among all of the estimations obtained for different SLEs that could be provided using all possible combinations of the tested tracers. The detected exceptions (see Figure 10 (II) (D)) could most likely be explained by the simplified linearized statistics used for design optimization; however, a Monte Carlo-based approach was later performed to calculate the non-symmetrical flux confidence intervals. 

It was interesting to see the possible effect of "general" and "partial" optimizations when the corresponding set of independently analyzed SLEs was considered as LEs in a PLE, followed by rigorous fitting of all of the simulated "experimental" data to the single metabolic model. The set of PLE-based experiments consisted of 5 earlier analyzed SLEs consisted of "generally" and "partially" optimized, as well as "randomly mixed" tracers were analyzed. As could be expected, the PLE consisted of "generally" and/or "partially" optimized LEs were among experiments with the most precisely resolved fluxes (with the maximal values of Σ η(0. 95) ). On the other hand, the real difference between summarized values of the η -function for all fluxes determined for these PLEs was very small (between 42. 7 and 44. 2 see, Additional file 3: Table SF-3. 9) and all of these values significantly exceeded the corresponding sum estimated for any SLEs among used in this comparison. So, the positive effect of PLE-based experiment provided with many LEs substantively differed in the used tracers was so significant, that further improvement of flux resolution due to additional optimization of experimental design could have rather marginal positive sense.
section_name: Optimal design for LEs using mixtures of [U-13 C]-, [1-13 C]-, and [U-12 C]-glucose isotopomers

section_content: C-MFA for SLEs/PLEs with singly 13 C-labeled-glucose tracers It has been repeatedly shown [19, 59, 72] that significantly improved resolutions of fluxes from different parts of a metabolic network can be achieved using commercially available or specially synthesized singly or multiply 13 Clabeled glucose isotopomers other than the previously applied and cheapest [1- 13 C] and [U- 13 C] variants. Moreover, the COMPLETE-MFA approach, using all singly labeled glucose tracers ([i- 13 C]-glucose, where i = 1, 2, …, 6) in the individual LEs of the PLE, was recently developed to evaluate metabolic fluxes for the metabolic model of wild-type Escherichia coli with a high precision [58]. Additionally, the PLE included two LEs with only [2-13 C]-glucose, and separate [3-13 C]-glucose in the medium provided the most precise fluxes for the E. coli model among all possible paired combinations of singly labeled glucose added as the sole tracer. 

It was of interest to evaluate the advantages of this approach in the case of another metabolic model used in the present study. Initially, the possible exploitation of all six singly labeled 13 C-glucose isotopomers, which were separately added to the medium as tracers for SLEs and PLEs, was investigated. The solutions of the corresponding constrained NLLSPs are summarized in Additional file 3: Table SF-3. 10 and Figure 8 : Columns 1 -6. The following conclusions were drawn from these data. Initially, the values of the normalized flux precision function varied significantly for individual fluxes in a tracer-dependent manner, and it was not the best singly labeled glucose tracer for the estimation of all fluxes in the assumed metabolic model. Indeed, according to the formally computed value of the D-factor (see (S -1. 5. 14)) for the six singly labeled tracers, the volume of the confidence ellipsoid for free fluxes had to be increased in the following order for the applied tracers:

i. e., the minimal volume of the confidence ellipsoid for free fluxes could be obtained using [3-13 C]-glucose as a singly labeled glucose tracer. Simultaneously, a different order of the tracers could be formed by increasing the summarized values of the normalized flux precision functions, Σ η(0. 95)

i. e., [4-13 C]-glucose was the best tracer according to this Σ η -based criterion; namely, [4-13 C]-glucose was the best tracer for determining 13 fluxes, and the corresponding values of the η -function for this tracer for 28 other fluxes were among the maximal values in the range of no more than 0. 02, which is a typical value for measurementspecific errors. Notably, the applicability of [4-13 C]-glucose as the tracer for the efficient resolution of primary branch points and reversibilities was previously demonstrated for a similar C. glutamicum metabolic model [49]. 

The flux resolution obtained in each of the six independent LEs could clearly be improved via PLE-based flux analysis if two or more LEs were combined in the PLE. The corresponding calculations for the possible pair combinations are provided (see Additional file 3: Table SF-3. 10), which resulted in the conclusion that the combination of [2-13 C]-glucose and [1-13 C]-glucose provided the most precise fluxes. In contrast, the use of [2-13 C]-glucose and [5-13 C]-glucose as the tracers in two independent LEs provided the least precise fluxes:

The results concerning the flux resolution obtained in experiments using singly labeled tracers demonstrated the model specificity of the experimental design optimization. Indeed, [4-13 C]-glucose was one of the best tracers for the resolution of a large portion of the fluxes for the assumed C. glutamicum model and was one of the worst tracers for determining precise fluxes in the E. coli model [58]. Furthermore, the best detected pair combination i¼1 j¼2 PLE for the resolution of fluxes in the model applied in the present study differed from the i¼2 j¼3 PLE combination, which is best for E. coli [58]. 

In full accordance with the published data on the E. coli model [58], a significant improvement of the flux resolution was detected under the COMPLETE-MFA approach, where all six singly 13 C-labeled glucose isotopomers were used separately as tracers in the individual LEs combined in the PLE (Figure 6 : Column 13). Indeed, the computed value of the normalized flux precision function for all fluxes was as follows:

The flux resolution detected in this PLE was the best among the above-described resolutions obtained in the present study. It was clear that this result was obtained through PLE-mediated 13 C-MFA due to the synergy of complementary information concerning the highly efficient resolution of the set of fluxes from the different parts of the metabolic model and the dependence of these sets on the different tracers used in individual LEs. Performing partial optimization of the experimental design for individual LEs for further PLE-based 13 C-MFA it seemed probable to improve the flux resolution. 

To test this possibility, partial optimization was conducted by searching the minimal approximation of free flux variances in the triangle composed of all possible mixed tracers, which consisted of the three singly labeled glucose isotopomers ([1- 13 It was interesting to mention, that all singly labeled glucose isotopomers used in these calculations could be considered as "partially optimized" for some free fluxes in the current metabolic model. According to this partial optimization, two PLEs were designed. The first consisted of three LEs with separately applied singly [1-, 3-, or 4- 13 C] labeled glucose. The second PLE consisted of six LEs, where, in addition to the three LEs from the previous PLE, three LEs used the optimized mixed paired combinations of glucose isotopomers. The corresponding PLEs were simulated, and 13 C-MFA was performed, resulting in solutions of NLLSPs with highly resolved fluxes (Figure 7 (D) and (E)). Concerning the PLE consisting of three LEs (Figure 8 : Column 10), the corresponding value of Σ η(0. 95) (1, 3, 4 ) PLE = 44. 6 Additional file 3: (Table SF-3. 10) significantly exceeded the best parameters obtained for the pairs of tracers, and was more than earlier obtained for PLEs from 5LEs with another tracers (see previous item). In turn, the flux precision obtained in the second PLE (Figure 8 : Column 12) was approximately identical (maybe even slightly exceeded) to the previously described result for COMPLETE-MFA performed with six different singly labeled isotopomers of glucose (Additional file 3:Table SF-

Notably, an exceptionally high flux resolution was obtained in both PLEs including 6 LEs. For all fluxes, η 0. 95 (u i ) > 0; i. e., the length of the 95% confidence interval did not exceed the value of the corresponding flux corrected by the scaling factor.
section_name: 13

section_content: The main aim of the present study was to extend the possibilities of the previously developed open-access OpenFLUX software for comprehensive 13 C-MFA. These extensions included (i) fitting the obtained experimental data, not only in SLEs but also in PLEs, to the assumed metabolic model; (ii) computing the flux parameters and providing the goodness-of-fit of the model's adequacy, followed by an estimation of the model's viability and its probable improvement; (iii) fine-tunable and convergence-controlled Monte Carlo-based approach to obtain distribution of optimized flux estimations followed by precise computing of flux confidence intervals; and (iv) conducting general and/ or partial experimental design through searching for the minimal value to characterize the average confidence interval length for all free fluxes and/or the minimal linear approximation of the targeted free flux variances, respectively. The considered examples demonstrated the specific features of these steps and their concerted essentiality for obtaining statistically verified results of the described 13 C-MFA provided with the help of OpenFLUX2. 

Introducing the normalized flux precision function, η γ (u i ; β), allowed for the quantitative characterization of the efficiency of the flux resolution at a confidence level of γ, with values of η close to "1" or "0" being obtained for efficiently or poorly resolved fluxes, respectively, depending on β as the scaling parameter, in particular. Moreover, the sum of η for all of the evaluated fluxes in the model, Σ η, could be a rather convenient parameter for conditional comparison of the flux precision achieved in different experiments, i. e., depending on the tracer(s) used. 

The goodness-of-fit test of the assumed metabolic model's adequacy is an essential and extremely important part of the statistically verified solution of the NLLSP. This test involved not only obtaining the χ 2 -statistically acceptable value of the Ξ(x input, x mea, σ mea, θ) objective function but also confirmation of the expectation concerning the Ν 0; 1 ð Þ distribution of its summands. Providing an insufficient goodness-of-fit test could result in mistaken flux estimations, and the computed contribution matrix could be helpful for improving the statistical properties of the obtained solutions. It is desirable to retain the summands corresponding to the most important measurements if some of the "outstanding" variance-weighted residuals must be deleted to satisfy the normalization criterion. 

The assumed metabolic model clearly influences the optimization of the experimental design and, ultimately, the precision of the flux estimation. Nevertheless, according to comparative calculations, the complementary parallel experimental labeling technique for metabolic flux analysis (COMPLETE-MFA [58] ) using a set of different labeled tracers, e. g., singly 13 C-labeled glucose isotopomers or combinations of mixed labeled tracer(s), chosen according to the partial optimization of the targeted fluxes, usually resulted in better flux precision compared with the resolution of fluxes computed from SLE data, which were provided according to a sophisticated general experimental design strategy. Only one simplified metabolic model, consisting of an L-lysine-producing C. glutamicum strain, was used for the simulation of experimental data and for comprehensive 13 C-MFA in the present study; however, the obtained general conclusions coincided well with the data reported by other groups working with other models (e. g., [48, 58] ). 

We hope that the developed OpenFLUX2 open-access software adjusted for comprehensive 13 C-MFA of SLE and PLE data will help to broaden investigations aimed at quantitatively estimating the cellular metabolic state with high precision. 

OpenFLUX2 is being released as open-source software. Regarding the citation of the OpenFLUX2 application, it is highly recommended that the present paper be cited, adding that this software is an extended version of the OpenFLUX open-source software described in [33].
section_name: Conclusions

section_name: Methods

section_content: The C. glutamicum metabolic model used as an example included catabolic reactions of the central metabolism, such as the Embden-Meyerhof-Parnas (EMP) and pentose phosphate (PP) pathways, the tricarboxylic acid (TCA) cycle, anaplerotic carboxylation, and the decarboxylation reaction of oxaloacetate and malate. Moreover, the pathways for the biosynthesis and transport of L-lysine and different extracellular co-products (glycine, trehalose, lactate, and α-ketoglutarate) were included. For glycine synthesis, two possible pathways, starting from serine and threonine [75, 76], were considered. The glyoxylate bypass (shunt) was assumed to be inactive [39]. Pools of pyruvate/PEP or oxaloacetate/malate were lumped, followed by the expression of reactions catalyzed by PEP/PYR carboxylase or by PEP carboxykinase/malic enzyme as an irreversible reaction. To achieve accurate accounting of CO 2 -associated carbon transfer, the reactions, accompanied by CO 2 production or consumption, were expressed in an explicit manner, including an anabolic reaction and a reaction involving CO 2 exchange with environment. Additionally, the linear consequence of the irreversible reactions was represented as a single reaction. The forward and reverse components of a bi-directional reaction were considered as two non-negative fluxes. The biosynthetic pathways for the following amino acids were expressed explicitly: (i) the amino acids whose mass isotopomer distribution (MID) was assumed to be measured and (ii) the amino acids whose synthesis was accomplished by CO 2 release. In Figure 2, for simplicity, the amino acid biosynthetic pathways are expressed schematically as a drain of the precursors to amino acids. Anabolic demand is represented in a slightly different manner than in previous studies [33, 68]. Specifically, a single biomass equation was designed by presenting the biomass composition of C. glutamicum [75] as the sum of the amino acids whose biosynthesis was expressed explicitly and the residual amounts of the precursors drained to produce biomass. The reactions involved in alanine, aspartate, and glutamate synthesis were used to map the MID of the metabolite onto the MID of the corresponding amino acid (S-type reactions according OpenFLUX(2) notation), whereas the anabolic demand to synthesize these amino acids was considered through precursors. The anabolic demand for lysine included both the lysine used in protein synthesis and the diaminopimelate used in cell wall synthesis, as previously described [75]. 

The resulting metabolic model contained 54 reactions and 36 balanced metabolites. The reactions, which only map the labeled distribution of the metabolite onto the corresponding amino acid (e. g., pyruvate to alanine), did not participate in the stoichiometric balance. Thus, the stoichiometric matrix S, dim(S) = (36 × 51), with 36 balanced metabolites, 51 unknown fluxes, and r = rank(S) = 36, was finally generated. As a result, 15 fluxes should be assigned as free. Seven fluxes were experimentally determined effluxes, which included biomass biosynthesis (1 flux), the effluxes of secreted products (5 fluxes), and the glucose uptake rate (1 flux). Unless otherwise stated, these fluxes were considered free fluxes constrained according to Eq. (S -1. 2. 2) with the experimentally determined parameters V mea eff ; σ mea eff ; thus, these fluxes subsequently formed the corresponding residuals in the SSR SLE f objective scalar-function Eq. (S -1. 3. 7), which was subjected to a least-squares minimization procedure. Five free fluxes were automatically assigned due to their accordance to the reverse reactions: (i) in the non-oxidative branch of the PP pathway (3 fluxes), (ii) catalyzed by glucose-6-phosphate isomerase (1 flux), and (iii) for the carbon dioxide intracellular exchange (1 flux). The three remaining free fluxes were previously automatically assigned using OpenFLUX software [33]. There were fluxes corresponding to irreversible reactions catalyzed by glucose-6-phosphate dehydrogenase in the PP pathway, PEP/PYR carboxylase in the PEP-PYR-oxaloacetate node, and glycine synthesis in the serine-glycine biosynthetic pathway. 

An isotopomer model was automatically built based on the EMU approach that simulated the MIDs of target compounds from the known MIDs of input substrates. The input substrates were specifically labeled glucose and naturally labeled CO 2 ( 13 C isotope abundance of 1. 07%). In total, 17 sets of matrix equations of EMU balances were used to calculate the 107 unknown EMUs from 15 known EMUs of input substrates. The application of the EMU approach reduced the number of unknown variables from 9,138 unknown scalar variables in the full isotopomer model to 360 scalar MID variables, corresponding to 107 unknown EMUs.
section_name: In silico experiments Metabolic and isotopomer balance models

section_content: Initially, the calculation of "true" flux values for the assumed metabolic model was performed using experimental data on effluxes and MIDs, together with their variances, available from the literature [33, 68]. The published labeling patterns of the proteinogenic amino acids were determined due to GC-MS-mediated selective ion monitoring of selected ion clusters, representing [M-57] fragments with the complete carbon skeletons of the amino acids. The simulated MIDs of a compound with n carbon backbone atoms was represented by a mass distribution column vector (MDV), whose elements corresponded to the fractional abundances (x m 0 þi ) of mass isotopomer m 0 + i

x m 0 þi ¼ 1, where m 0 is the molecular weight of an unlabeled compound. As the published MIDs were uncorrected in accordance with the natural mass isotopomer abundance, all the simulated EMU variables were modified for mass interference from non-carbon backbone isotopes in dependence on the chemical structure of the derivatized substance according to the method developed in [77]. Then for each fragment, the first (n + 1) simulated mass isotopomers (m 0, m 1, …, m n ) were normalized to unit followed by truncation of the isotopomers number to the dimension of the x mea MID vector before performing least-square analysis. The measured MIDs dim x mea MID À Á ¼ 3 Â 1 ð Þ, i. e., only (m 0, m 1, m 2 ) were presented for derivatized mass isotopomers fragments of Ala-260, Val-288, Thr-404, Asp-418, Glu-432, Ser-390, Phe-336, Tyr-466, Tre (trehalose)-361 isotopomers, and MIDs dim

for (m 0, m 1 ) of Gly-246. At the same time, an extremely small MID variance of 0. 15% for mass isotopomer fractional abundances was indicated in [68] and accepted in [33], which ultimately led to rather large and statistically nonacceptable values of the objective function, i. e., the variance-weighted sum of squared residuals, calculated at the point of convergence (Ξð ⌢ θÞ≡2⋅ SSR SLE f ð ⌢ θÞ, according to the notations of the present publication (see Additional file 1: SF-1. 3. ). These MID variances were assumed to be significantly underestimated. Indeed, the minimal published GC-MS measurement errors are usually approximately 0. 2-0. 4 mol% [78], and values of this magnitude have been used for variance weighting in studies where 13 C-MFA has been applied [53, 58, 79]. To achieve statistically acceptable model fitting for the MIDs measured in [68] and to demonstrate all of the essential stages of the statistical analysis of the solution to the constrained NLLSP provided by the designed OpenFLUX2 software, the MID variances were assumed to be equal to 0. 15 mol% for mass isotopomer fractional abundances (instead of 0. 15%, as in [68] ).
section_name: Data used for estimating true flux values in the assumed model

section_content: To generate new "experimental data" for in silico LEs, the OpenFLUX(2) forward simulation option was used, which calculated the MIDs of selected EMU variables from the assumed metabolic and isotopomer networks as well as the known label states of input substrates ( 13 C-tracer(s)) and the assumed values of the fluxes. Unless otherwise stated, the true flux values, u true (θ true ), were used for direct "experimental" data simulation in the present study. An isotopic purity of 99% was assumed for all of the specifically labeled carbons in the input substrates, and natural enrichment (1. 07%) was assumed for other carbons. Again, the simulated MIDs of a substance with n carbon backbone atoms were represented by MDV with (n + 1) elements, x m 0 þi ; i¼1; 2; …; n X n i¼0 x m 0 þi ¼ 1. The natural mass isotopomer abundance of non-carbon-backbone atoms was generated according to the method [77] that resulted in significant increasing of the elements number in MDV for each tested substance, but OpenFLUX2 automatically truncated these mass isotopomers up to the first (n + 1) (with or without their normalization according to the user's choice) (Additional file 3: Table SF-3. 2). Then, mass isotopomer fractional abundances of lower than 0. 04 (i. e., ≤ 4 mol%) from the total of 1 (i. e., 100 mol%) for each derivatized fragment were excluded from the set of simulated "experimental data," thus modeling the limited sensitivity of the MS equipment. It is important to mention, that this procedure of the experimental data generation resulted in dependence of x mea MID vector dimension for each "measured" substance not only on its specific number of carbon backbone atoms, but on the different 13 C-labeled tracers used in the simulated experiments. At the final stage of "experimental" MID generation, mass isotopomer fractional abundances were corrupted with N 0; σ mea MID À Á distributed noise, where σ mea MID È É i ¼ 0:4 mol% ½ ; i ¼ 1; 2; …; w MID. As a result, a set of non-normalized, noisy MS "experimental data," m mea MID, that were corrected for the presence of natural isotopes, was generated for each LE according to the described procedure. 

Another type of "experimental data" used during the in silico experiments was measured effluxes, V mea eff. To generate these data, the "true" value of each efflux was corrupted with noise distributed as Nð0; fσ mea eff j Þ, where fσ mea eff g j was the experimentally determined standard deviation for each j-th efflux [68]. As a result, the full set of the "experimental data," m mea ¼ m mea MID V mea eff, was generated for the in silico-simulated LE.
section_name: Generation of new "experimental data"

section_content: Flux estimation, identifiability, and goodness-of-fit analyses were performed using OpenFLUX2 software. To determine the global minimum of the constrained NLLSP, as a rule, 100 independent iterative trials, starting from randomly selected points from a feasible constrained domain in the free flux variation space, ℜ p, were applied. In special cases, e. g., when the minimal value of the Ξ function was achieved in only a few trials, the number of iterative trials was increased to 300 to verify the detected global minimum. Only those iterative trials were used for analysis, which were terminated by satisfaction of the termination criteria without constraint violation (see Additional file 1: SF-1. 3. ). Monte Carlo-based approach included previously in OpenFLUX, but significantly modified in the OpenFLUX2 software for convenience of the computation process tuning and control of convergence, was used for estimation of the flux confidence intervals. Usually the most reliable results could be obtained if the variant when "multi runs per trials" approach was chosen for estimation of the optimized flux parameter distributions followed by determination of flux confidence interval borders according to "discarding" strategy after confirmation of the borders convergence (see Additional file 1: SF-1. 7. for details).
section_name: Flux estimation and statistical analysis

section_content: OpenFLUX2 (http://sourceforge. net/projects/openflux2) requires Java and MATLAB, including the Optimization and Statistics Toolboxes. The current version of the OpenFLUX2 software was tested using Java 6 (Sun Microsystems, Inc., Santa Clara, CA, USA) and MATLAB 7. 12. 0 (R2011a; MathWorks Inc., Natick, MA, USA) software, together with the Optimization Toolbox (version 6. 0) and the Statistics Toolbox (version 7. 5), on the Microsoft Windows 7 Professional (2009; Microsoft Corp, Redmond, WA, USA) platform on a PC equipped with a 3. 2 GHz CPU and 4 Gb of RAM. Using these computation facilities, the calculation of 51 fluxes for the applied metabolic model (see Methods) requires approximately 20 minutes and 2 hours for SLEs and a PLE consisting of 5 LEs, respectively. The confidence interval estimations for the 51 optimized fluxes using the Monte Carlo approach (see Methods) takes approximately 30-70 and 60-120 hours for the SLEs and the PLE noted above, respectively. The computation time significantly depends on the used tunable control parameters and linearly depends on the number of the provided trials. In general, the flux estimation and its statistical analysis via OpenFLUX2 can be performed within one to four days, depending on the type of LE. OpenFLUX2 was also tested on the Microsoft Windows XP (Professional × 64 edition, 2003; Microsoft Corp, Redmond, WA, USA) platform. The following additional software packages were used: Windows Microsoft Excel 2003 (Microsoft Corp, Redmond, WA, USA) for metabolic model formulation and for the generation of Additional file 3 and OriginPro 9. 1 (Originlab Corp, Northampton, MA, USA) for the generation of ternary plots, and box charts. Furthermore, several in-house MATLAB scripts were employed to generate "experimental data" and to visualize data related to the NLLSP solving process.
section_name: Software requirements

section_name: Acknowledgments

section_content: First, we are extremely grateful to all authors of [33], particularly to Dr. LE Quek for programming the software and to Dr. JO Krömer for supervising the study that resulted in the release of OpenFLUX as open-source software. We are thankful to Dr. Yousuke Nishio ( Ajinomoto Co., Inc. ) for helpful discussions, particularly upon the initiation of the 13 C-MFA study at our institute. MSS and LIG are thankful to Profs. E Heinzle, M Reuss, and C Wittmann for their supervision and their excellent introduction to the basic problems of metabolic flux analysis in the 2012 Braundwald Course on Biosystems Engineering -Bioreactors and Cell Factories. Moreover, MSS and LIG are thankful to all of the organizers and lecturers, particularly to Prof. W Weichert, Dr. K Nöh, and S Niedenführ, of the 3 rd Advanced Course on 13 C-based MFA (Jülich, 2013), for the creative atmosphere and intensive training in all major topics relevant to the design, execution and analysis of the labeling experiments using the high-performance cutting-edge simulator 13CFLUX2 as well as for helpful discussions related to recent and future developments in 13 C-MFA.
section_name: Acknowledgments

section_content: http://www. microbialcellfactories. com/content/13/1/152

section_name: Additional files

section_content: Additional file 1: General description of 13 C-MFA steps performed by OpenFLUX2. The PDF-format file contains the schematic description of the mathematical and statistical approaches used for comprehensive metabolic flux analysis via OpenFLUX2 software. 

Additional file 2: Additional figures illustrated results of performed in silico LEs. The PDF-format file contains additional illustrations of the in silico experimental results, denoted as Figure x in the text. 

Additional file 3: C. glutamicum metabolic model. The XLS-format file includes the C. glutamicum metabolic model in a format that is readable by OpenFLUX(2); the "experimental data" (MIDs and efflux measurements with their errors) used in the in silico experiments; the main results of the in silico experiments, including the estimated optimal flux values, confidence intervals and normalized flux precision function values; and the results of the "general" and "partial" optimization of the experimental design.
section_name: Additional files

section_content: The authors declare that they have no competing interests. 

Authors' contributions MSS contributed his expertise to the design of the final version of OpenFLUX2 software and programmed the primary portion of the modifications. He also contributed to the writing of the manuscript, particularly to Additional file 1. LIG provided her expertise in the mathematical modeling of the fluxes and contributed to the writing of the paper. SSR provided the primary part of the OpenFLUX2 simulations, programmed visualization of Monte Carlo-based determination of flux confidence intervals, and contributed to the preparation of the primary portion of the figures and tables in the paper. DAP programmed the special subprograms in OpenFLUX2 to compute the inverse Hessian matrix at the stage of optimization of the experimental design. SI contributed his expertise to the 13 C-MFA and mathematical simulations of mass isotopomer distributions as well as to the writing of the paper. SVM designed and supervised the study and coordinated writing of the paper. All authors have approved the final manuscript.
section_name: Competing interests