 Methodology
 Open Access
 Published:
Parallel isotope differential modeling for instationary 13C fluxomics at the genome scale
Biotechnology for Biofuels volume 13, Article number: 103 (2020)
Background
A precise map of the metabolic fluxome, the closest surrogate to the physiological phenotype, is becoming progressively more important in the metabolic engineering of photosynthetic organisms for biofuel and biomass production. For photosynthetic organisms, the stateoftheart method for this purpose is instationary 13C fluxomics, which has arisen as a sibling of transcriptomics or proteomics. Instationary 13C data processing requires solving highdimensional nonlinear differential equations and leads to large computational and time costs when its scope is expanded to a genomescale metabolic network.
Result
Here, we present a parallelized method to model instationary 13C labeling data. The elementary metabolite unit (EMU) framework is reorganized to allow treating individual mass isotopomers and breaking up of their networks into strongly connected components (SCCs). A variable domain parallel algorithm is introduced to process ordinary differential equations in a parallel way. 15fold acceleration is achieved for constantstepsize modeling and ~ fivefold acceleration for adaptivestepsize modeling.
Conclusion
This algorithm is universally applicable to isotope granules such as EMUs and cumomers and can substantially accelerate instationary 13C fluxomics modeling. It thus has great potential to be widely adopted in any instationary 13C fluxomics modeling.
Background
With the arrival of the postgenome era, 13C fluxomics has matured as a stateoftheart approach to derive in vivo metabolic flux information in parallel with transcriptomics, proteomics and metabolomics [1,2,3]. This method captures very important and unique information reflecting intracellular physiology that could never be achieved by other omics techniques and is pivotal to metabolic engineering of microbes for biofuel and bioproduct [4,5,6,7]. Its irreplaceability qualifies it as an extraordinarily powerful tool for exploring metabolic flux and has encouraged progressively wide application to a broad variety of organisms, such as photosynthetic organisms, fungi or mammalian cells [8,9,10].
13C fluxomics can be divided into two categories. One is stationary 13C fluxomics, which assumes a steadystate metabolic and isotopic labeling environment and deals with algebraic balance equations of mass and isotopes [3, 11]. Model construction and flux estimation can be performed with different isotope granules, such as isotopomers, cumomers and elementary metabolite units (EMUs) [12,13,14,15]. This capability has excited a wave of tool development, such as 13CFLUX2, FiatFlux, and WUFlux [16,17,18,19,20].
However, the isotopic steadystate assumption does not hold true for many circumstances, for example: (1) fedbatch conditions where the isotopic steady state cannot readily be reached in a short time [21, 22]. (2) Continuous cultivation conditions with a pulsed substrate supply [23,24,25]. In particular, steadystate 13C fluxomics becomes inefficient where the substrate is a onecarbon compound, such as CO_{2} for photoautotrophic cultivation for biofuel and biomass conversion, since the steadystate isotopic distribution is binomial and independent of the flux values [26, 27].
To circumvent the limits of steadystate 13C fluxomics, the second type of 13C fluxomics—isotopically instationary 13C fluxomics (INSTFluxomics) has been invented to be applied to systems that are in a metabolic steady state while being isotopically instationary [23, 24, 26,27,28,29]. This technique has been implemented in software such as INCA, OpenMebius [30, 31], and extensively adopted to investigate the intracellular physiology of cyanobacteria, microalgae and plants [6, 32,33,34].
Unlike the stationary case, instationary 13C fluxomics is required to model a large set of ordinary differential equations (ODEs) instead of algebraic equations. In particular, the equations contain highdimensional derivative variables over the flux values and pool sizes and inevitably incur large computing and time costs [10, 24]. For example, modeling a realistic central carbon model of E. coli with constant time integration method requires several hours [30]. This problem grows significantly when a genomescale metabolic network is considered [35, 36]. Such models usually have hundreds of metabolites and reactions, resulting in a steep increase in the number of equations [35, 36]. The motivation of this paper is to find a way to more quickly model a large set of isotope ordinary differential equations, especially those for a genome scale metabolic network. Here, we present the first parallel method for instationary fluxomics analysis. We reorganized the mass isotopomers of the EMU to model them individually. This treatment facilitates downstream isotopic network decomposition and simplification. An algorithm for parallelization in the variable domain is used to model the ODE systems in a parallel way. Implementing this method has achieved tens of fold improvement with constantstepsize ODEs and ~ fivefold improvement for adaptivestepsize ODEs. Since the method is universally applicable to labeling frameworks such as EMUs and cumomers, we expect that this method will benefit all 13C fluxomics communities.
Results
SCC decomposition and EMU reorganization in a mass isotopomer network
Here, we use a toy network previously reported to show the framework (Fig. 1 and Table 1) [13]. According to the generation relationship, we can draw the mass isotopomer labeling network in a directed graph G_{m}. Then, the unnecessary part that could not contribute to the measured mass isotopomer distribution (MID) is removed from this mass isotopomer labeling network. Thus, the computation complexity can be significantly reduced.
If two mass isotopomers are interdependent, they should be modeled together and simultaneously. The interdependency is that the nodes belong to one strongly connected component (SCC). To model the mass isotopomers separately and in parallel, the classic Tarjan algorithm is employed to decompose the whole network G_{m} into small pieces of SCCs (Fig. 2a) [37, 38].
The SCCs are organized by their topological sort to reveal their dependency, as shown in Fig. 2a (see Method description section). As such, the mass isotopomers are transformed into different SCCs, which still contain all original mass isotopomers, as shown in Fig. 2b. The motivation to perform such reorganization is that the whole network can be split into SCCs with smaller scale and large quantity and similar manipulation is universally applicable to other isotope granules like cumomer and isotopomer [12, 14]. In addition to the toy network, this study also uses an E. coli metabolic network and a genomescale metabolic network of Synechococcus 2973, which is modified from imSyn593 reconstruction [36]. The number of SCCs is 14, 45 and 98 for the 3 networks, respectively.
Additionally, the adjacent SCCs with the same weight can be combined in a headtotail way, to form a larger unit that could be modeled simultaneously similarly. We also call such a unit as SCC. To assess the impact of such combination, we set a parameter λ representing the minimum mass isotopomer quantity of an SCC with each of the new SCCs being larger than λ. As λ goes up, the number of SCCs declines.
Parallelization for a constantstep ODE at the genome scale
The parallelization strategy is shown in Fig. 3. In simulating the ODE, the value of each time step of one SCC depends on not only its own value in the previous step, but also the values of its parent SCCs in the previous step. Once the values of one SCC for time Ti have been computed, these values are delivered immediately to the threads where their downstream SCCs are waiting. The delivery paths of these values (Fig. 3a) are the same as the cascade relationship in the mass isotopomer network. With this strategy, all SCCs can be computed simultaneously step by step, and the corresponding time is significantly shorter than the sum of times for each SCC individually (Fig. 3b).
There are at least two kinds of ODE methods applicable according to the stepsize choice: a constantstepsize method and an adaptivestepsize method [39, 40]. The constantstepsize method is well suited for parallelization. A 4thorder Runge–Kutta method with a constant step size is implemented with a nonparallel technique and a parallel technique on the genomescale carbon mapping models (Fig. 4) [36, 39]. The evolution modeling is carried out with tensorbased and vectorbased method. The following data are from the vectorbased method since it got a significantly faster speed.
The model encompasses 78 free fluxes in 174 total fluxes and 189 free pool size in 266 metabolites, resulting in 407 carbon transition reactions (Additional file 1). The bicarbonate uptake rate was set as 10 mmol/gDW/h before scaling, while the growthrelated dilution was not considered. The pool sizes were randomly set within the physiological range from micromole per liter to millimole per liter. The MID of whole molecules of 14 metabolites was set as the measured MID for network simplification as previously reported (2phosphodglycerate, 6phosphodgluconate, ribose5phosphate, dglucose6phosphate, phosphoenolpyruvate, pyruvate, sedoheptulose7phosphate, succinate, malate, 3phosphodglycerate, dfructose6phosphate, citrate, sedoheptulose1,7bisphosphate and dfructose1,6bisphosphate, 5methyltetrahydrofolate, and shikimate) [36].
Then, a total of 819 mass isotopomers were identified in the transitive closure of 14 metabolites. They were transformed into 98 pieces of SCCs with up to 2.16 × 10^{5} ODE equations. The ordinary equations were simulated from 0 to 10 s. The step size for constantstepsize method is set to be 0.005 s. The absolute and relative tolerance was set as 10^{−9} and 10^{−7}. 3 different sets of the flux distribution were used for the S.2973 network. The ODEs parameters are set in Table 2 and flux values are documented in the Additional file 2. The speed of the nonparallel and parallel methods is compared in Fig. 5. The parallel method obtains an average 15fold acceleration over the nonparallel method, as shown in Fig. 5. The λ is 10, generated by a grid search for best speedup.
Parallelization for the adaptivestepsize ODE at the genome scale
In contrast, the adaptivestepsize method is potentially hard to parallelize since the step size depends on the current value and may differ for different SCCs [40]. Therefore, a universal step size requires computation on all SCCs simultaneously.
Fortunately, the labeling curve of any mass isotopomer is not complicated and can be mimicked by a cubic spline [41]. In particular, the sum of the mass isotopomers of the same EMU equals one while the sum of the mass isotopomers’ derivatives relative to particular flux or pool size equals zero. So, the dynamic curves of the mass isotopomers and the derivatives belonging to the same metabolite have similar steepness and require similar sampling density, even though there is a wide distribution in the concentrations of metabolites. We calculated the trajectory of all mass isotopomers and the derivatives of EMUs of different sizes in a mediumscale carbon metabolic network of E. coli (the Additional file 3). The calculation was repeated for 3 sets of flux rate and pool size. The result indicates the evolution curve of different mass isotopomers and their derivatives from the same EMU have the similar changing amplitude in the same timescale. This phenomenon is exemplified by the curves of all EMUs of AcCoA and isocitrate in Fig. 6a–c, from one dataset. More data are provided in the Additional file 4.
Figure 6d–f shows the mass isotopomers of mass 0 (m_{0}) and their derivatives of different EMUs of the same metabolite for 18 metabolites, according to above parameter settings. The curves belonging to the same metabolites are indeed similar to each other, implying that the step size calculated for m_{0} of one EMU from one metabolite is suitable for other mass isotopomers of any EMU of the same metabolite. More data are provided in the Additional file 4.
This property enables us to carry out parallelization even in the adaptivestepsize method. First, a set of m_{0} values is constructed from the first few SCCs, which contain at least one EMU of each metabolite. This set is used to calculate the step size with an adaptive method, which is simultaneously delivered to the SCCs to carry out the parallelized constantstepsize modeling.
The speedup in the adaptivestepsize method is shown in Fig. 5. The parallel method obtains an average speedup of 5 times. λ does not exhibit an observable effect because the timelimiting step is the process of calculating the step size for the systems, which is not correlated with λ. The adaptivestepsize method can identify and use the maximal step size in each step while producing an error of an allowed magnitude. This results in a significant reduction in the total number of steps and a decrease in the running time.
Stiffness problem is an interesting problem for instationary 13C fluxomics. As shown in Eq (7), the ODEs for each SCC are one order and their derivatives are determined by the value of SCC at right side and the SCCs with lower topological sort. Since the value of the SCCs with lower topological sort are known, the derivatives from these components are constant over time. The truncate error will be fixed once the stepsize fixed. As such, the stiffness of SCC_{i} is major resulted from the square matrix M_{i} in Eq (7), which is multiplying with SCC_{i}. The stiffness is characterized by this matrix’s stiffness ratio, which is the ratio between maximal value and minimal value of the real components of the eigenvalues of the matrix.
We have calculated the eigenvalues of the corresponding matrix of different SCCs with 3 sets of flux distribution on E. coli network. All the eigenvalues are negative real number, which shows that the isotope differential equation is stable. The maximal stiffness ratio is less than 10^{4} for all the SCCs with λ ranging from 5 to 1000, which means the stiffness of these SCCs is not a serious problem.
Discussion
Modelling a large set of isotope ordinary differential equations in a parallel and faster way have a potentiality to be highly used by the researchers who develop new 13C fluxomics method and will help to utilizes 13C fluxomics to explore more species. To this end, a new framework for modeling the isotope is generated. A parallel strategy is realized for isotope ordinary equations modeling and achieved observable speedup. Some factors affecting the speed of isotope ordinary equations modeling are discussed.
Specifically, we have put forward a new framework for modeling the isotope fate within the atom transition network. This framework is in parallel with current isotope granule EMU, cumomer and bondmer and easy to be implemented. The modeling is based upon a new combination of mass isotopomer. The mass isotopomer generation relationship has been abstracted as a mass isotopomer network. Mass isotopomer network is an important concept here. The equations can be significantly simplified by preserving the transitive closures of the measured mass isotopomers and cutting out the rest parts. This is like what have been done for cumomer and EMU. The left mass isotopomer network is then decomposed into different SCCs, which is the basic unit for modeling. The modeling can be implemented in two ways, one of which is tensorbased and the other vectorbased. The tensorbased method is analogue to the cumomer algorithm from Weitzel’s work [38]. The vectorbased method is analogue to the EMU block algorithm from Young’s work [26].
A parallel strategy of modeling isotope ordinary equations has been realized with the constantstepsize and adaptivestepsize integration methods. Their comparison is as shown in Fig. 5. This strategy is essentially universally applicable and can be adapted in different isotope modeling method and benefit to the 13C fluxomics community. For the test cases, the vectorbased method shows systems advantage in modeling speed over the tenserbased method. This is reasonable because the vectorbased method requires only the product of two vectors, while the tensor algorithm requires an additional component. These results confirm the wide applicability of our parallel method.
Many other factors also impact the modeling speed. In addition to the hardware performance such as the number of CPU cores, the speed of the modeling is also affected by the number of EMU reactions and the flux distribution value. When the difference between the reaction values within one flux distribution goes down, the size of time point calculated by the adaptive step size method become smaller and the modeling can be quickly ended. This information is also instructive for modeling based upon EMU and cumomer.
Conclusion
Isotopically instationary 13C fluxomics has emerged as the gold standard to obtain a precise picture of the fluxome for photosynthetic microorganisms utilizing a onecarbon substrate. The increased computational and storage demands may hinder easy adoption when a genomescale fluxome is required. The parallel algorithm for isotope ODEs here is found to be a successful strategy to promote the speed of instationary 13C fluxomics, regardless of whether constantstepsize or adaptivestepsize modeling is used. This parallel strategy utilizes the cascade relationship of the isotope granules balance equations. This property also holds true for any other isotope granules, such as EMUs or cumomers. It is essentially universally applicable and will benefit to the 13C fluxomics community.
Method description
Reorganization of EMU mass isotopomers
An EMU can be defined as a specific subset of metabolite atoms [13]. A mass isotopomer is a group of isotopomers that is classified according to the number of heavy isotopes rather than the position [42]. The mass isotopomers here are considered individually, unlike those previously treated in one EMU as a whole. One mass isotopomer is treated as one node in the graph. When a mass isotopomer A is a precursor of mass isotopomer B, we define a directed arrow starting from node A and ending on node B.
Subsequently, according to EMU equations, a directed graph G_{m} (V_{m}, E_{m}) of the mass isotopomer network can be drawn [38]. The vortex V_{m} represents the mass isotopomers, while the directed edge E_{m} connects the reactant isotopomers to the product mass isotopomers in an EMU reaction. To remove the unnecessary mass isotopomers, the graph G_{m} is first transformed into its transposed graph G ^{T}_{m} by reversing each edge of G_{m}. Then, the transitive closures of the measured mass isotopomers are identified by a Floyd–Warshall algorithm [43].
Then, the Tarjan algorithm has been implemented to decompose G_{m} into different SCCs, which is a subgraph where there exists at least one directed path for any pair of its vertices [37]. A solitary mass isotopomer shall be treated as an SCC consisting of a single node. If SCC A contains a node with an edge directed toward another node in SCC B, then SCC A has a prioritized topological sort over B, and SCC A is called the parent SCC of SCC B. Thus, for all edges from A to B, node A appears before node B. The SCCs with higher sort rely on the SCCs with lower sort, while the latter do not rely on the former.
As a heavier mass isotopomer will never contribute to a lower mass isotopomer in terms of an EMU, one SCC contains mass isotopomers of the same weight. If one directed edge connects SCC A to SCC B, then we define SCC A as having a topological sort less than that of SCC B [38]. The content of an SCC is as follows:
where i represents the weight of the mass isotopomer, while j represents the topological sort of the SCC for the same weight.
Tensorbased modeling of the mass isotopomers of SCCs
The instationary framework presented here is based upon the mass isotopomers of SCCs. Generally, any matrix manipulation of isotopes can be represented by a combination of SCC reactions. One SCC reaction with n reactants and one product has the following form:
The generation of the MID of SCC_{P} can be calculated as the product of the transition tensor of SCCs of all reactants, such as
Here, \(SCC_{p} \cdots SCC_{{R_{n} }}\) is the vector of mass isotopomers of product and reactant SCCs sorted by their weights and topological sort. Q is a transition tensor with order equal to the number of reactant and product SCCs and dimensions equal to the dimensions of the SCCs. The precise definition of the transition tensor can be described as follows:
The consumption of the mass isotopomer can be formulated by an eliminating matrix. Its content can be defined as follows:
The dimension is equal to the number of mass isotopomers for the corresponding SCC. The ODEs of mass isotopomer SCCs have the following form:
SCC^{ik,jk} is an SCC with the j_{k}th topological sort and the i_{k}th weight. v_{p} is the flux value of reactions producing certain mass isotopomers of SCC^{ik,jk}. v_{c} is the flux value of reactions consuming certain mass isotopomers of SCC^{ik,jk}. v_{inp} is the flux value of the input reactions. SCC^{i1,j1} and so on are the SCCs whose element is involved in v_{p} to produce SCC^{ik,jk}. i_{p} is the number of reactants of v_{p}, while j_{p} is the topological sort. Q ^{(i1,j1) (i2,j2),…, (ip,jp)}_{p} is the transition tensor of v_{p} as defined in Eq. (6) to adapt the SCC instead of the EMU, which produces SCC^{ik,jk} from SCC^{i1,j1} and so on. E_{c} is the eliminating matrix of v_{c}. Q_{inp} is the transition tensor of input reactions. C_{ik, jk} is a diagonal matrix whose diagonal elements are the metabolite pool size associated with SCC^{ik,jk}.
Implicit differentiation of Eq. (3) with respect to the free fluxes and pool size generates the differential equations of firstorder derivatives of the measured mass isotopomer equations. The initial conditions for these derivatives are set as zero to solve these equations.
Vectorbased modeling of mass isotopomers of SCCs
Like EMU [13], mass isotopomer SCCs can also be modeled directly in a vector way as Eq (7)
where \(SCC_{{{\text{j}}_{1} }}\) and \(SCC_{{{\text{j}}_{2} }}\) are the SCCs which contribute to SCC_{i} through a single molecule transition. The single molecule transition from certain SCCs to target SCC_{i} can be characterized in the matrix \(M.\)
\(SCC_{{{\text{j}}_{k1} }}\) and \(SCC_{{{\text{j}}_{k2} }}\) are the SCCs which contribute to SCC_{i} through a double molecule transition. \({ \circledast }\) is a userdefined vector product specific to the set of \(SCC_{{{\text{j}}_{k1} }} , SCC_{{{\text{j}}_{k2} }}\) and \(SCC_{i}\). For different set of \(SCC_{{{\text{j}}_{k1} }} , SCC_{{{\text{j}}_{k2} }}\) and \(SCC_{i}\), the content of \({ \circledast }\) is different. \(M_{{{\text{j}}_{n} }}\) is set as a unit matrix to be compatible with such configuration. The sensitivities equations can be organized in a similar way.
Parallel algorithm for the ODEs
The decomposed SCCs are a cascade system, and heavier SCCs rely on the value of lower SCCs. Generally, the ODEs are solved sequentially to obtain the timeseries data of mass isotopomers. Therefore, parallel implementation of the SCC differential equations substantially accelerates the forward simulation. The parallel strategy is schematically shown in Fig. 4.

1.
Each SCC evolution is carried out in one individual thread.

2.
For the first thread, a 4thorder Cash–Karp Runge–Kutta scheme is employed to solve the equations about the first SCC^{1,1}.

2.1
Each time step of h_{k} and each value of SCC ^{1,1}_{k} of this process are preserved and communicated to the downstream threads.

2.1

3.
For the thread of SCC^{i,j}, the current value SCC ^{i,j}_{k} is calculated based upon SCC ^{i,j}_{k1} together with h_{k} and SCC ^{im,jm}_{k} received from previous threads where i_{m} and j_{m} are less than i and j, respectively. The algorithm is a typical explicit 4thorder Runge–Kutta method as described in Eq. (9) [36, 39].

3.1
Each time step of h_{k} and each value of SCC ^{i,j}_{k} of this process are preserved and communicated to the downstream threads:

3.1

4.
Repeat (3).
Algorithm implementation and parameter setting
The parallel algorithm is performed via ExecutorService. The data communication between threads is executed efficiently through ConcurrentHashMap. Constraintcompatible initial flux distributions are generated by the Pythonbased sampler OptGPSampler in COBRApy [44]. A clear protocol to guide the method development has been included in the Additional file 5. The toy network is adopted from a previous report [13]. The central carbon metabolism network of E. coli is described in the Additional file 3. The genomescale metabolic model of S.2973 was modified from imSyn593 from a previous study [36]. The computer is equipped with 2 Xeon E7 4820 V2 2.2G CPUs with 8 cores and 512G memory card, which can support 32 threads at most.
Availability of data and materials
All data generated or analyzed during this study are included in this published article and its additional files.
Abbreviations
 AKG:

Alphaketoglutarate
 AcCoA:

ActylcoA
 Ery4P:

Erythrulose4phosphate
 EMU:

Elementary metabolite unit
 Fru6P:

Fructose6phosphatase
 Glc6P:

Glucose6phosphate
 GAP:

Glyceraldehyde 3phosphate
 ICit:

Isocitrate
 Mal:

Malate
 OAA:

Oxaloacetate
 ODE:

Ordinary differential equation
 SCC:

Strongly connected component
 PEP:

Phosphoenolpyruvate
 Pyr:

Pyruvate
 Rib5P:

Ribose 5phosphate
 Rul5P:

Ribulose5bisphosphate
 Ser:

Serine
 Suc:

Succinate
 Xul5P:

Xylose5phosphate
 MID:

Mass isotopomer distribution
References
 1.
Zamboni N, Sauer U. Novel biological insights through metabolomics and 13Cflux analysis. Curr Opin Microbiol. 2009;12:553–8.
 2.
Crown SB, Antoniewicz MR. Publishing 13C metabolic flux analysis studies: a review and future perspectives. Metab Eng. 2013;20:42–8.
 3.
Niedenführ S, Wiechert W, Nöh K. How to measure metabolic fluxes: a taxonomic guide for 13C fluxomics. Curr Opin Biotechnol. 2015;34:82–90.
 4.
Hollinshead WD, Rodriguez S, Martin HG, Wang G, Baidoo EE, Sale KL, Keasling JD, Mukhopadhyay A, Tang YJ. Examining Escherichia coli glycolytic pathways, catabolite repression, and metabolite channeling using Deltapfk mutants. Biotechnol Biofuels. 2016;9:212.
 5.
Liu N, Qiao K, Stephanopoulos G. (13)C metabolic flux analysis of acetate conversion to lipids by yarrowia lipolytica. Metab Eng. 2016;38:86–97.
 6.
Wu C, Xiong W, Dai J, Wu Q. Genomebased metabolic mapping and 13C flux analysis reveal systematic properties of an oleaginous microalga Chlorella protothecoides. Plant Physiol. 2015;167:586–99.
 7.
Zhao L, Zhang H, Wang L, Chen H, Chen YQ, Chen W, Song Y. (13)Cmetabolic flux analysis of lipid accumulation in the oleaginous fungus Mucor circinelloides. Bioresour Technol. 2015;197:23–9.
 8.
Christen S, Sauer U. Intracellular characterization of aerobic glucose metabolism in seven yeast species by 13C flux analysis and metabolomics. FEMS Yeast Res. 2011;11:263–72.
 9.
Metallo CM, Walther JL, Stephanopoulos G. Evaluation of 13C isotopic tracers for metabolic flux analysis in mammalian cells. J Biotechnol. 2009;144:167–74.
 10.
Qian X, Zhang Y, Lun DS, Dismukes GC. Rerouting of metabolism into desired cellular products by nutrient stress: fluxes reveal the selected pathways in cyanobacterial photosynthesis. ACS Synth Biol. 2018;7:1465–76.
 11.
Heux S, Berges C, Millard P, Portais JC, Letisse F. Recent advances in highthroughput (13)Cfluxomics. Curr Opin Biotechnol. 2017;43:104–9.
 12.
Wiechert W, Mollney M, Isermann N, Wurzel M, de Graaf AA. Bidirectional reaction steps in metabolic networks: III. explicit solution and analysis of isotopomer labeling systems. Biotechnol Bioeng. 1999;66:69–85.
 13.
Antoniewicz MR, Kelleher JK, Stephanopoulos G. Elementary metabolite units (EMU): a novel framework for modeling isotopic distributions. Metab Eng. 2007;9:68–86.
 14.
Schmidt K, Carlsen M, Nielsen J, Villadsen J. Modeling isotopomer distributions in biochemical networks using isotopomer mapping matrices. Biotechnol Bioeng. 1997;55:831–40.
 15.
Zupke C, Stephanopoulos G. Modeling of isotope distributions and intracellular fluxes in metabolic networks using atom mapping matrixes. Biotechnol Prog. 1994;10:489–98.
 16.
Zamboni N, Fischer E, Sauer U. FiatFlux–a software for metabolic flux analysis from 13Cglucose experiments. BMC Bioinform. 2005;6:209.
 17.
Quek L, Wittmann C, Nielsen LK, Kromer JO. OpenFLUX: efficient modelling software for 13Cbased metabolic flux analysis. Microb Cell Fact. 2009;8:25.
 18.
Weitzel M, Nöh K, Dalman T, Niedenführ S, Stute B, Wiechert W. 13CFLUX2—highperformance software suite for 13Cmetabolic flux analysis. Bioinformatics. 2012;29:143–5.
 19.
He L, Wu SG, Zhang M, Chen Y, Tang YJ. WUFlux: an opensource platform for 13 C metabolic flux analysis of bacterial metabolism. BMC Bioinformatics. 2016;17:444.
 20.
Birkel GW, Ghosh A, Kumar VS, Weaver D, Ando D, Backman TW, Arkin AP, Keasling JD, Martín HG. The JBEI quantitative metabolic modeling library (jQMM): a python library for modeling microbial metabolism. BMC Bioinform. 2017;18:205.
 21.
Antoniewicz MR, Kraynie DF, Laffend LA, Gonzalezlergier J, Kelleher JK, Stephanopoulos G. Metabolic flux analysis in a nonstationary system: Fedbatch fermentation of a high yielding strain of E. coli producing 1,3propanediol. Metab Eng. 2007;9:277–92.
 22.
Ahn WS, Antoniewicz MR. Metabolic flux analysis of CHO cells at growth and nongrowth phases using isotopic tracers and mass spectrometry. Metab Eng. 2011;13:598–609.
 23.
Noh K, Wahl A, Wiechert W. Computational tools for isotopically instationary 13C labeling experiments under metabolic steady state conditions. Metab Eng. 2006;8:554–77.
 24.
Noh K, Wiechert W. Experimental design principles for isotopically instationary 13C labeling experiments. Biotechnol Bioeng. 2006;94:234–51.
 25.
Toya Y, Ishii N, Nakahigashi K, Hirasawa T, Soga T, Tomita M, Shimizu K. 13Cmetabolic flux analysis for batch culture of Escherichia coli and its pyk and pgi gene knockout mutants based on mass isotopomer distribution of intracellular metabolites. Biotechnol Prog. 2010;26:975–92.
 26.
Young JD, Walther JL, Antoniewicz MR, Yoo H, Stephanopoulos G. An elementary metabolite unit (EMU) based method of isotopically nonstationary flux analysis. Biotechnol Bioeng. 2008;99:686–99.
 27.
Young JD, Shastri AA, Stephanopoulos G, Morgan JA. Mapping photoautotrophic metabolism with isotopically nonstationary (13)C flux analysis. Metab Eng. 2011;13:656–65.
 28.
Nakajima T, Yoshikawa K, Toya Y, Matsuda F, Shimizu H. Metabolic flux analysis of the Synechocystis sp. PCC 6803 Δ nrtABCD mutant reveals a mechanism for metabolic adaptation to nitrogenlimited conditions. Plant Cell Physiol. 2017;58:537–45.
 29.
Abernathy MH, Yu J, Ma F, Liberton M, Ungerer J, Hollinshead WD, Gopalakrishnan S, He L, Maranas CD, Pakrasi HB. Deciphering cyanobacterial phenotypes for fast photoautotrophic growth via isotopically nonstationary metabolic flux analysis. Biotechnol Biofuels. 2017;10:273.
 30.
Kajihata S, Furusawa C, Matsuda F, Shimizu H. OpenMebius: an open source software for isotopically nonstationary 13Cbased metabolic flux analysis. Biomed Res Int. 2014;2014:627014.
 31.
Young JD. INCA: a computational platform for isotopically nonstationary metabolic flux analysis. Bioinformatics. 2014;30:1333–5.
 32.
Jazmin LJ, O’Grady JP, Ma F, Allen DK, Morgan JA, Young JD. Isotopically nonstationary MFA (INSTMFA) of autotrophic metabolism. Methods Mol Biol. 2014;1090:181–210.
 33.
Sake CL, Metcalf AJ, Boyle NR. The challenge and potential of photosynthesis: unique considerations for metabolic flux measurements in photosynthetic microorganisms. Biotechnol Lett. 2019;41:35–45.
 34.
Ma F, Jazmin LJ, Young JD, Allen DK. Isotopically nonstationary 13C flux analysis of changes in Arabidopsis thaliana leaf metabolism due to high light acclimation. Proc Natl Acad Sci USA. 2014;111:16967–72.
 35.
Gopalakrishnan S, Pakrasi HB, Maranas CD. Elucidation of photoautotrophic carbon flux topology in Synechocystis PCC 6803 using genomescale carbon mapping models. Metab Eng. 2018;47:190–9.
 36.
Hendry JI, Gopalakrishnan S, Ungerer J, Pakrasi HB, Tang YJ, Maranas CD. Genomescale fluxome of Synechococcus elongatus UTEX 2973 using transient (13)Clabeling data. Plant Physiol. 2019;179:761–9.
 37.
Tarjan RE. Depthfirst search and linear graph algorithms. SIAM J Comput. 1972;1:146–60.
 38.
Weitzel M, Wiechert W, Noh K. The topology of metabolic isotope labeling networks. BMC Bioinform. 2007;8:315.
 39.
Dormand JR, Prince PJ. A family of embedded RungeKutta formulae. J Comput Appl Math. 1980;6:19–26.
 40.
Lourakis MI. A brief description of the LevenbergMarquardt algorithm implemented by levmar. Foundation Res Technol. 2005;4:1–6.
 41.
Horl M, Schnidder J, Sauer U, Zamboni N. Nonstationary (13)Cmetabolic flux ratio analysis. Biotechnol Bioeng. 2013;110:3164–76.
 42.
Christensen BB, Nielsen J. Isotopomer analysis using GCMS. Metab Eng. 1999;1:282–90.
 43.
Hougardy S. The FloydWarshall algorithm on graphs with negative cycles. Inform Process Lett. 2010;110:279–81.
 44.
Megchelenbrink W, Huynen MA, Marchiori E. optGpSampler: an improved tool for uniformly sampling the solutionspace of genomescale metabolic networks. PLOS ONE. 2014;9:e86587.
Acknowledgements
Not applicable.
Funding
This work was supported by the National Science Foundation of China NSFC (Grant 31760254), the Joint Fund of the Natural Science Foundation of China and the Karst Science Research Center of Guizhou Province (Grant No. U1812401) and the Doctoral Scientific Research Foundation of Guiyang University (Grant GYUZRD [2018]018).
Author information
Affiliations
Contributions
TS conceived and designed the study. TS, ZL, YW designed the algorithm. ZZ, ZL wrote code. TS, ZZ wrote the manuscript. YM, ZZ, ZC, JH prepared the figures. XX, YY supervised the study. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
All authors consented on the publication of this work.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Additional file 1. A genomescale metabolic network of
Synechococcus2973 modified fromimSyn593.
Additional file 2
. 3 different sets of the flux distribution of theS.2973 network.
Additional file 4. The trajectory of all mass isotopomers and the derivatives of EMUs of different sizes in
E.coli for 3 set of flux distribution.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Zhang, Z., Liu, Z., Meng, Y. et al. Parallel isotope differential modeling for instationary 13C fluxomics at the genome scale. Biotechnol Biofuels 13, 103 (2020). https://doi.org/10.1186/s13068020017375
Received:
Accepted:
Published:
Keywords
 Instationary metabolic flux analysis
 Parallel differential equations modeling
 Genomescale metabolic flux analysis
 13C fluxomics
 Mass isotopomer network