Abstract
Background
Flux balance analysis (FBA) together with its extension, dynamic FBA, have proven instrumental for analyzing the robustness and dynamics of metabolic networks by employing only the stoichiometry of the included reactions coupled with adequately chosen objective function. In addition, under the assumption of minimization of metabolic adjustment, dynamic FBA has recently been employed to analyze the transition between metabolic states.
Results
Here, we propose a suite of novel methods for analyzing the dynamics of (internally perturbed) metabolic networks and for quantifying their robustness with limited knowledge of kinetic parameters. Following the biochemically meaningful premise that metabolite concentrations exhibit smooth temporal changes, the proposed methods rely on minimizing the significant fluctuations of metabolic profiles to predict the timeresolved metabolic state, characterized by both fluxes and concentrations. By conducting a comparative analysis with a kinetic model of the CalvinBenson cycle and a model of plant carbohydrate metabolism, we demonstrate that the principle of regulatory on/off minimization coupled with dynamic FBA can accurately predict the changes in metabolic states.
Conclusions
Our methods outperform the existing dynamic FBAbased modeling alternatives, and could help in revealing the mechanisms for maintaining robustness of dynamic processes in metabolic networks over time.
Background
Systems biology paradigm has provided insights in the maintenance of robustness for biological processes involving a multitude of interconnected elements (e.g., genes, proteins, metabolites) [1]. In addition, recent advances in metabolomics have provided a large amount of highly reproducible data [24], allowing reconstruction and analysis of genomescale metabolic networks [5]. These developments in metabolomics technologies have challenged computational systems biology with the need to accurately describe the dynamics of metabolic networks in order not only to glean the flux rates at different time points, representing the temporal flux (re)distribution, and the interdependent metabolic profiles, but also to identify key elements for metabolic engineering [610].
Metabolic flux analysis (MFA) has propelled the development of computational methods for analysis of metabolic networks [11,12]. Flux balance analysis (FBA), as one of the most prominent of the MFA methods, is based on linear programming (LP) whereby a given objective function (e.g., biomass yield) is optimized under the assumption that the system operates at steady state under the constraints given by the stoichiometric matrix [11,1316]. By optimizing an objective function, the linear program identifies one feasible flux distribution from the set of fluxes satisfying the constraints imposed by the massbalance equations and reaction bounds [15]. Consequently, the biological implications of the optimal flux distribution depend on the choice of the objective function [17]. Maximization of biomass is one of these functions, which is assumed particularly suitable for microbial models [18]. For eukaryotic cells (e.g., in plants) where biomass or yield may not be the primary goal, a different objective function has to be determined. For instance, cellular maintenance at minimal efforts has been proposed to be one of the alternatives [19]. Nevertheless, finding an adequate objective of a metabolic network, and especially a sub network related to particular metabolic processes, remains a problem of ongoing interest [18,20]. However, the steadystate assumption on which FBA is based precludes the analysis of the dynamics of metabolite concentrations and flux (re)distribution. Furthermore, the classical FBA ignores the possibility that perturbed metabolic networks may not immediately regulate towards the (assumed) optimal objective.
Based on the hypothesis that fluxes in metabolic networks, altered by removal of a reaction, undergo a minimal redistribution compared to those of the wild type, minimization of metabolic adjustment (MOMA) [21] and regulatory on/off minimization (ROOM) [22] have been devised as two contending alternatives for analysis of perturbed metabolic network models. MOMA predicts the flux redistribution which has the smallest Euclidean distance to the wild type flux distribution obtained by FBA, while ROOM minimizes the number of (significant) flux changes from the wild type flux distribution. As a consequence, large modifications in single fluxes are prevented in MOMA; however, such large modifications may be required for rerouting metabolic flux through alternative pathways [22], which has been observed in experiments [23]. Existing studies have demonstrated that ROOM outperforms MOMA and FBA in the flux prediction of the final metabolic steady state, albeit, in the particular case of pyruvate kinase knockout in E. coli [22].
The analysis of metabolic network dynamics has traditionally been facilitated by models based on ordinary differential equations (ODEs) which require a large amount of information for simulating the temporal metabolic changes [9,10]. To this end, the phenomenological parameters of specific enzyme kinetics (e.g., MichaelisMenten or Hill) have to be determined by accurate measurements of enzyme activities and datafitting to experimentally obtained (timecourse) data. In turn, the obtained fits are often used to make predictions and draw conclusions based on the postulated kinetic model.
On the other hand, dynamic FBA (DFBA) offers the alternative to predict timeresolved metabolic profiles with limited knowledge of enzyme kinetics [68]. Unlike, the analyses based on FBA, which focus on the steadystate behavior, DFBA offers the means to analyze transient (nonsteady) states. In addition, DFBA has been combined with MOMA, resulting in the MDFBA approach which has been employed for predicting the dynamics of photosynthetic metabolism and positing hypotheses about its robustness under different CO_{2}^{} and water conditions [7]. MDFBA extends the MOMA hypothesis from minimal redistribution of metabolic fluxes in perturbed metabolic networks to a minimal fluctuation of the profile of metabolite concentrations over time.
We point out that the mechanisms describing the temporal changes in the metabolic state, characterized by the metabolite concentrations and flux rates, are principal to the notion of robustness, as used with DFBAbased approaches. For instance, the posited mechanism in MDFBA is that metabolic networks operate to minimize the fluctuations in metabolic concentrations over time. This suggests alterations to the definition of robustness as a property that maintains system function in the case of external perturbations [24,25] to include the internal perturbations due to changes in the dynamic state of the system's elements. External perturbations arise due to changing environmental conditions, such as stress conditions (environmental perturbations), but also by changes in the structure of a metabolic networks (e.g., caused by gene knockouts, also known as genetic perturbations). In contrast, internal perturbations are due to temporal changes of the metabolic state characterized by the coupling of metabolite concentrations and flux distributions. Considering mass actions kinetics, which states that the flux rate of a reaction is proportional to the product of the concentrations of the participating reactants [26]. As a result, even in this simplest kinetic law, a change in metabolite concentrations may affect flux rates, which, in turn, as a result of the mass balances have an effect on the concentrations. Here, we investigate a set of biochemically plausible hypotheses which can be used to characterize the mechanisms responsible for maintaining the metabolic network robustness due to internal perturbations. These mechanisms can in turn be used to simulate the dynamics of a given metabolic network relying purely on the stoichiometry and a limited amount of information regarding the phenomenological constants.
Driven by the idea of mutually influencing system elements, central to the systems biology paradigm, we argue that minimal fluctuations of metabolic profiles may only represent one possible explanation of the temporal changes in metabolic state. In this study, we design the ROOMbased DFBA approach (abbreviated to RDFBA) by coupling the principle of ROOM with DFBA. RDFBA extends the premise of the ROOM approach, which relies on significant flux changes, by considering the minimization of the total number of significant changes of metabolite concentrations. Furthermore, the coupling of ROOM and DFBA renders it possible to use the advantages of ROOM compared to MOMA in a dynamic setting.
In addition, we modify and extend the proposed RDFBA and the existing MDFBA approaches to consider minimizing fluctuations not only in concentrations but also in flux levels, which result in seven proposed methods based on different optimization functions and programming formulations. Including classical DFBA and the two known variants of MDFBA [8], a total of ten different approaches are presented to predict and analyze the dynamics of metabolite concentrations and flux levels over time. Finally, the accuracy of a dynamic FBAbased approach in simulating the dynamics of metabolic networks must be established based on the comparison of the obtained results and the outcome of a welldefined kinetic model for a given metabolic network. Here, by comparing the results of applying DFBA, MDFBA, RDFBA and their extensions to those of the kinetic models for the CalvinBenson cycle and plant carbohydrate metabolism, our analysis discriminates between the different mechanisms which result in the apparent metabolic network dynamics and robustness. Our quantitative and qualitative comparative analyses demonstrate that the extensions based on RDFBA outperform the existing DFBAbased approaches.
Results and discussion
Here we describe the seven proposed methods and the comparison of their performance with two kinetic modelsof the CalvinBenson cycle and of the plant carbohydrate metabolism. Since the proposed methods build upon the DFBA approaches, we provide a brief overview of the mathematical apparatus used in their formulation. The suite of proposed methods and their relation to the existing DFBA approaches are depicted in Figure 1.
Figure 1. Overview of the ten approaches, which are used to analyze the dynamics of metabolite and flux profiles in the CalvinBenson cycle and the plant carbohydrate metabolism. The novel methods proposed in this study are framed in red.
We point out that DFBAbased methods and kinetic modeling constitute independent approaches. The DFBAbased methods involve performing a constrained optimization over a time period to approximate the dynamics of a system with stoichiometric constraints only, which include relative quantities of reactants and products of the particular reaction. In contrast, a kinetic model requires all kinetic parameters to simulate the dynamics of the system. However, kinetic modeling is often hampered by incomplete knowledge of the underlying enzymekinetics and their associated parameter values [17].
In addition, it is necessary to determine the extent to which the dynamics of both metabolite profiles as well as flux predictions can be successfully predicted from DFBAbased approaches. The outcome of testing this hypothesis will then suggest the most likely mechanism which describes the robustness of the system to internal perturbations over time. We argue that such approach goes beyond the analysis of steady states typical for metabolic control analysis (MCA) [27] and structural kinetic modeling as its extension [28].
DFBA
DFBA overcomes the main drawback of the classical FBA which precludes the analysis of the dynamic behavior of a networkthe steadystate assumption. Mahadevan and coworkers introduced two DFBA formulationsstatic and dynamic [6]. The static optimization (SOA) involves first dividing the time period of interest into uniform time intervals and then solving the instantaneous optimization problem at the beginning of each time interval, followed by integration to compute the metabolite concentration over time. In contrast, the general dynamic optimization approach (DOA) is formulated as follows:
where x and v are vectors of metabolite concentrations and reaction fluxes over time, S denotes the stoichiometric matrix, with rows corresponding to metabolites and columns to reactions of the metabolic network described by S, and t is the time. The minimum and maximum allowable fluxes of each reaction and metabolite concentrations are defined by v_{min }and v_{max }and x_{min }and x_{max}, respectively. The vector X_{0 }gives the initial concentration for the set of metabolites. The formulation of the DOA approach results in a nonlinear program (NLP) if nonlinear constraints and/or a nonlinear objective function are included.
The DOA involves optimization over the entire time period t to obtain timeresolved flux rates and metabolite concentrations. The optimization is rendered computationally feasible by parameterizing the dynamic equations with the help of orthogonal collocation on finite elements [29]. To this end the time period of interest is divided into a finite number of intervals, named finite elements. Furthermore, the metabolite concentrations x and flux levels v are parameterized at the roots of an orthogonal polynomial (e.g., Legendre polynomial) within each finite element [6]. For readers not familiar with orthogonal collocation on finite elements an instructive example of how this algorithm works is provided in Additional file 1.
Additional file 1. Tutorial Orthogonal Collocation on Finite Elements. Tutorial for the orthogonal collocation on finite elements is included. An illustrative example demonstrating the usage of this method in approximating solutions to ODEs is also presented.
Format: PDF Size: 150KB Download file
This file can be viewed with: Adobe Acrobat Reader
We point out that the number of variables in each optimization step of SOA is smaller compared to that of DOA, allowing scalability to larger networks. However, SOA does not allow dynamic formulation, and the remainder of the methods focus on DOA.
MDFBA
Luo and coworkers developed an approach called MDFBA which combines MOMA with DFBA
[7]. By employing the MOMA hypothesis, the objective function of the dynamic optimization
approach of DFBA is altered to a minimization of the Euclidean distance between metabolite
concentrations at adjacent orthogonal roots
with M representing the number of orthogonal roots, and δ, the Dirac delta function (see Additional file 1).
Extended versions of MDFBA
We extend the MDFBA method so that the objective function, subject to minimization, is given by the sum of Euclidean distances between metabolite concentrations and reaction rates at adjacent orthogonal roots:
where F represents the number of reactions in the network. This approach, denoted MDFBA_{CF}, extends the MDFBA approach by considering the hypothesis of the minimal fluctuation of flux levels in addition to minimal fluctuation of the profile of metabolite concentrations. In general, only limited knowledge of the enzymekinetic rate laws is available for an incorporation of kinetic expressions in the constraints. On the other hand, the inclusion of only stoichiometric constraints for the majority of the reactions can have the effect that flux rates and metabolite concentrations do not change as expected due to their coupling based on the underlying kinetics. Therefore, stabilizing only the profile of metabolite concentrations in the objective function, like in the basic MDFBA approach, may lead to a large variation of flux rates over time. The inclusion of both metabolite concentrations and flux rates in the objective function renders a more constrained optimization problem, which could avoid large fluctuations in the predictions.
Moreover, to facilitate discrimination of the different mechanisms yielding a particular timeresolved metabolic state, we consider a third type of the MDFBA method, whereby the objective function, subject to minimization, includes only the fluctuation of flux rates, as follows:
RDFBA
In contrast to MOMA, ROOM minimizes the total number of significant (large enough) flux changes from the wildtype flux distribution. In our approach, we combine ROOM with DFBA, naming this modeling approach RDFBA. Here, the ROOM hypothesis is extended from a minimization of significant flux changes to a minimization of fluctuation of the profile of metabolite concentrations, according to the assumption of the basic MDFBA approach. Unlike MDFBA, in RDFBA this is realized by minimizing the significant concentration changes at adjacent orthogonal roots. In contrast to MDFBA, which captures smooth changes over time (by the Euclidean distance), RDFBA prevents a large number of small (significant) changes over time. Instead, it allows a large concentration change at a few time points (orthogonal roots) and a relatively constant concentration of the metabolites elsewhere, which is bounded by thresholds determining the significant change. Finally, there are two possibilities to capture concentration changes between two time pointsthe continuous, specified by the usage of a distance measure (e.g., Euclidean distance), and the binary, based on an appropriate definition of significant change. The proposed formulation of RDFBA captures the binary case, as the remaining alternative to specify dynamic changes. From a biological perspective, the possibility for large concentration changes in a small time range, as specified in RDFBA, may be necessary to capture the difference in concentrations between two steady states in a bistable system [3032].
Two different types of the RDFBA approach are considered: The first includes the basic ROOM approach with integer constraints, whereas the second relies on relaxing the integer constraints. Due to the relaxation of the integer constraints, the programming problem in ROOM becomes linear in comparison to mixedinteger linear program (MILP) in the basic ROOM approach [22]. The first RDFBA approach, based on MILP ROOM, is described as follows:
where w^{u}, w^{l }are the thresholds determining significant changes of metabolite concentrations (u = upper bound, l = lower bound) and γ_{x}, ε_{x }are the relative and absolute ranges of tolerance, respectively.
Due to the constraints defining the thresholds, the programming problem becomes a mixedinteger nonlinear programming problem (MINLP). These problems are difficult to solve precisely as they couple the combinatorial nature of mixedinteger programming (MIP) problem and the computational complexity of solving NLP problems [33].
Consequently, we formulate a second RDFBA approach based on relaxed integer constraints, described as follows:
Compared to the first formulation of RDFBA, the relative and absolute ranges of tolerance (γ_{x}, and ε_{x}) are set to be zero, due to the relaxed integer constraints for y_{i,j}.
Extended versions of RDFBA
As for the MDFBA approach, we also extend the proposed RDFBA by considering the minimization of fluctuation of flux levels in addition to the minimized fluctuation of the profile of metabolite concentrations at adjacent orthogonal roots. We denote this approach as RDFBA_{CF}:
where b^{u}, b^{l }are the thresholds determining significance of the flux levels (u = upper bound, l = lower bound) and γ_{v}, ε_{v }are the relative and absolute ranges of tolerance, respectively. Due to the inclusion of the flux rates and metabolite concentrations in the description of the thresholds, RDFBA_{CF }is solved based on mixedinteger nonlinear programming, by the first RDFBA formulation, and based on nonlinear programming, according to the second. This extension minimizes fluctuations of both, metabolite concentrations and flux rates, over time, similarly to MDFBA_{CF}. In comparison, in RDFBA_{CF }large concentration and flux rate changes between two adjacent orthogonal roots for some metabolites and/or reactions are possible, which is precluded in MDFBA_{CF}.
Additionally, we consider a third objective function for the RDFBA approach, whereby only fluctuation of the profile of flux levels of the network over time is minimized:
The proposed optimization functions of the MINLP and NLP variants of RDFBA, fluxbased RDFBA and RDFBA_{CF }allow a comparison of the necessity of stabilizing the system's profile of concentration, flux or a combination of both to predict dynamics in a metabolic network.
Comparison of the proposed methods and a kinetic model of the CalvinBenson cycle
The CalvinBenson cycle consists of three phases relying on energy supply in form of ATP and redox elements (NADP/NADPH) and supply of CO_{2}: (1) carboxylation, during which the enzyme RuBisCo adds CO^{2 }to ribulose1,5bisphosphate (RuBP) to get two molecules of phosphoglycerate (PGA), (2) reduction, converting the obtained PGA into 1,3diphosphoglycerate (DPGA) and glyceraldehyde3phosphate (GAP), and (3) regeneration, which recovers RuBP after several intermediate steps from ribulose5phosphate (Ru5P) [34,35]. The CalvinBenson cycle is a wellstudied system of the plant metabolism due to its importance in providing precursors for the biomass synthesis which is necessary for plant growth. More than 15 kinetic models were already developed to analyze the dynamics of this important pathway [36]. In this study, we use the simplified model of Zhu et al. [37] to describe the dynamic behavior of the Calvin cycle, depicted in Figure 2 and Table 1. The concentrations of ATP, NADPH, phosphate (P), and CO_{2 }are assumed to be positive and unbounded, yielding the simplified equations shown in Table 1[38]. For the present analysis, other factors which affect photosynthetic metabolism, such as illumination and temperature, are assumed to be constant.
Figure 2. Simplified model of the Calvin cycle. The model includes six metabolites (Ru5P, RuBP, PGA, DPGA, GAP, and Sink) and seven reactions (v_{1 } v_{7}).
Table 1. Reactions in the simplified model of the Calvin cycle
With the simplified model of the Calvin cycle, we examine and compare the performance of the proposed RDFBA approach together with its extensions, the existing MDFBA approach and our proposed extensions, as well as the classical DFBA. In total, ten different approaches, based on different constraintbased formulations, are used in the analysis of the dynamics of metabolite and flux profiles in the Calvin cycle (Figure 1). With each approach, the concentration of the metabolites ribulose5phosphate (Ru5P), ribulose1,5bisphosphate (RuBP), 3phosphoglycerate (PGA), 1,3diphosphoglycerate (DPGA), glyceraldehyde3phosphate (GAP) and sink as well as the rate of the seven reactions (v_{1}v_{7}) are investigated over a time period of ten seconds, which is sufficient due to the fast settling of the system in a steady state [39].
For the analysis, the time period is divided into five finite elements, and the variables (metabolite concentrations and reaction rates) are parameterized at the roots of the fifth order Legendre polynomial. Maximizing the sink production is assumed as systemic objective for the DFBA approach, capturing the utilization of PGA and GAP to build the metabolite precursors necessary for sucrose and starch synthesis. For the remaining methods, the objective is modified to minimizing the fluctuation of the profile of metabolite concentrations/fluxes, according to the hypotheses of MDFBA and RDFBA.
The results obtained by the different methods are compared to the outcome of a kinetic model of the Calvin cycle, described in [37]. For the kinetic model as well as for the other compared methods, the same initial conditions are used (RuBP = 2.0 mmol l^{1}; PGA = 2.4 mmol l^{1}; DPGA, GAP, Ru5P, sink = 1.0 mmol l^{1}). In addition, for the DFBAbased approaches the limits for the GAP concentration over time are set to the outcome of the kinetic model with an added tolerance of +/ 0.1 mmol l^{1}. The solutions to MINLP formulation of RDFBA (and its extensions) are obtained with γ_{x }= 0.4, γ_{v }= 0.2, and ε_{x }= 0.01, ε_{v }= 0.05, representing the relative and absolute ranges of tolerance, respectively. For each approach, we determine the residual sum of squares (RSS), quantifying the correspondence between the constraintbased and kinetic modeling approaches. Furthermore, we determine the temporal evolution of Kendall τ for the metabolic states of the constraintbased and kinetic modeling approaches in order to qualitatively examine the coupling between metabolic profiles and distribution of fluxes for the two types of modeling approaches.
RDFBA yields accurate timeresolved predictions of metabolite profiles
Figure 3 shows the simulations results of two metabolites, RuBP and DPGA, for the ten compared methods. For the proposed basic RDFBA approach, minimizing the concentration fluctuations, the results differ between the variants: while the NLPbased variant predicts that the concentration of DPGA are almost constant over time, the MINLPbased formulation results in changes of the concentration of DPGA over the 10second time interval. The results of basic MDFBA show that the concentration of DPGA slightly decreases during the first three seconds, and is almost constant following this time period.
Figure 3. Modeling results of concentrations of RuBP (red) and DPGA (blue) in the simplified model of the Calvin cycle by the different approaches. The residual sum of squares values of each approach and metabolite are presented in the bottom right of the corresponding subfigure.
The effect of minimizing the flux changes is especially pronounced in the two variants of the fluxbased RDFBA. Here, the results for the DPGA concentration are similar to those of the basic MDFBA and its extensions. The behavior of the DPGA concentration predicted by both fluxbased RDFBA variants at the beginning of the simulation is very close to the kinetic modeling results, which can be explained by the coupling between fluxes and concentrations giving rise to particular temporal profiles.
However, this behavior is lost in (NLPbased) RDFBA_{CF}, which combines the basic and fluxbased variant of RDFBA. In the variants of RDFBA_{CF}, the changes in flux rates and metabolite concentrations are weighted equally, as
given by the objective function:
The results from the DFBA approach for the DPGA and RuBP concentration over time are very similar, with small changes in the first three seconds, leading to almost zero concentrations at the end of the simulation. The predicted RuBP concentrations of the basic MDFBA and its extensions are nearly constant over time; on the other hand, with the proposed NLPbased RDFBA approach, changes in concentrations are obtained in the first half of the simulated time interval, followed by a constant concentration. The concentration of this metabolite at the end of the simulation differs between the proposed MINLPbased RDFBA variants. We point out that experimental measurements in isolated cells and chloroplasts have resulted in low levels of RuBP [3941].
Comparing these results with the outcome of kinetic modeling shows that for these two metabolites MINLPbased RDFBA predicts the most accurate temporal profiles. The constant concentration in the second half of the simulated interval of the NLPbased RDFBA is very close to that of kinetic modeling, for both, the basic RDFBA and RDFBA_{CF}. To quantify these observations, we determine the residual sum of squares, by calculating the distance between kinetic modeling results and predicted results of the different approaches, see Table 2. The proposed approaches exhibit the smallest RSS values for these two metabolites, demonstrating the predictive power of RDFBA not only on MINLP but also on NLP.
Table 2. Residual sum of squares values for each metabolite
The simulated results of the concentrations of the other four metabolites are depicted in Figures S1  S4 in Additional file 2. Due to the fact that we have integrated rather strict constraints for GAP, the predicted concentrations of this metabolite over time are very similar for all approaches and also compared to the kineticbased model. For Ru5P, the results obtained by the basic MINLPbased RDFBA are in accordance with the kinetic modeling results and outperform the remaining approaches; on the other hand, the NLPbased RDFBA variants predict a constant amount of Ru5P of about 1 mmol l^{1 }except for the fluxbased variant, which also predict a nearly consumed Ru5P after four seconds similar to the kineticbased results. Surprisingly, the predictions based on the MINLPbased RDFBA_{CF }are inconsistent with all other approaches. Here, the minimization of significant fluctuation changes of flux levels and metabolite concentrations results in a short decrease followed by an increase of the Ru5P concentration. We believe that this behavior of the RDFBA_{CF }variants is due to the overconstraining of both fluxes and concentrations.
Additional file 2. Additional Figures. Modeling results of concentrations of Ru5P, PGA, GAP and sink as well as results of reaction rates of v_{1}v_{7 }in the simplified model of the Calvin cycle by the different approaches.
Format: PDF Size: 69KB Download file
This file can be viewed with: Adobe Acrobat Reader
In addition, the PGA concentration closest to the kinetic modeling results is predicted by the MINLPbased RDFBA approach which minimizes significant flux changes. Contradictory to the results from all other methods, the fluxbased RDFBA, formulated as NLP, predicts an increase of the PGA concentration.
Finally, we analyze the concentration of sink over time. Intriguingly, the sink concentrations are almost constant in the predictions of all NLPbased variants of MDFBA and RDFBA. However, for all MINLPbased RDFBA variants as well as for the DFBA approach, an increase of the sink concentration is observed, which is in line with the kinetic modeling results. Quantitatively the basic MINLPbased RDFBA approach predicts the best simulation results for the concentration of sink.
Taken together, the results of the RSSbased analysis demonstrate that MINLPbased RDFBA outperforms MDFBA in four out of six cases in predicting metabolite concentrations over time.
(De)coupling of timeresolved flux predictions
The flux levels over time, also known as timeresolved flux predictions, obtained by the different methods are depicted in Figures S5  S11 in Additional file 2. By inspecting the results, it becomes apparent that the predictions of the NLPbased variants of both MDFBA and RDFBA differ from the outcome of kinetic modeling. The reactions rates are very small for these six methods. Comparing the basic methods with the fluxbased variants of MDFBA and RDFBA show smoothing of curves, which is also appears in the results of MDFBA_{CF }and RDFBA_{CF}. Furthermore, it can be observed, that after four seconds for all NLPbased RDFBA variants the reaction rates are close to zero.
Interestingly, the reaction rates results obtained by DFBA are very close to the kinetic modeling results, as indicated by the RSS values in Table 3. The accuracy of DFBA in the prediction of reaction rates demonstrates that the nonperturbed network may operate towards maximization of sink production. However, we point out that for the reactions v_{1}, v_{3}, v_{4}, v_{6 }and v_{7 }the predicted reaction rates by DFBA are bounded by the corresponding v_{max}; the latter do not constrain the predictions in the rest of DFBAbased methods.
Table 3. Residual sum of squares values for each reaction
We argue that to maintain robustness in case of an internal perturbation, the system may opt to optimize additional objectives which are captured by the DFBAbased extensions. Indeed, we observe that the predicted rates obtained by the MINLPbased RDFBA_{CF }are closer to the kinetic results than all MDFBA approaches, as quantified by the RSS.
In the results of the MINLPbased RDFBA variants, we observe slight local fluctuation which may be attributed to the presence of a large null space in the process of searching an optimal solution of the objective function. A higherorder Legendre polynomial and, therefore, more orthogonal roots in each finite element may reduce these fluctuations. Altogether, these results indicate the decoupling of metabolite concentrations and flux rates in the determined solutions, i.e., that the metabolite concentrations and flux levels change independently of each other. However, it is possible that other optimal solution may reside in the vicinity of the results from the kinetic modeling.
To qualitatively examine the extent to which both metabolic concentrations and flux distributions are in agreement with the outcome of kinetic modeling over time, we determine the Kendall τ correlation for the corresponding metabolic states. The results in Figure 4 demonstrate that, qualitatively, the DFBA approach captures best the coupling of flux distributions and metabolic concentrations over time, followed by the fluxbased MDFBA and three of the proposed approachesthe fluxbased RDFBA and RDFBA_{CF }based on MINLP as well as MDFBA_{CF}.
Figure 4. Timeresolved Kendall correlation coefficient for the outcome of kinetic modeling and each of the constraintbased modeling approaches. The Kendall correlation coefficient at each time point integrates the vector of metabolite concentrations with that of the flux rates.
In addition, these results demonstrate that while the NLPbased RDFBA methods yield good predictions of metabolite concentrations, the predicted flux distributions over time are decoupled from the metabolite profiles. Interestingly, the same observations for the decoupling of fluxes from metabolite concentrations hold for all MDFBAbased approaches, except the fluxbased MDFBA. These findings indicate that some results of DFBAbased approaches warrant caution in biological interpretation and applications to metabolic engineering via flux control. Due to the highly coupled regulation of flux and concentrations in metabolic networks [28], it is not surprising that the MINLPbased formulations, especially the flux based RDFBA and RDFBA_{CF}, outperform the rest of the investigated approaches.
Comparison of the proposed methods and a kinetic model of the plant carbohydrate metabolism
In this section, we analyze the diurnal dynamics of central carbohydrates in leaf cells of wild type C_{3}plants with the help of the ten implemented methods. In parallel to the wild type, analyses are performed for the inv4 mutant, which is defective in the vacuolar invertase gene, AtβFruct4 (AT1G12240). Like in the previous section, the performance of the methods to accurately simulate the dynamics of the metabolite concentrations and fluxes is determined based on the RSS values between the results of the constrainedbased approaches and the kinetic models for the wild type plant and the inv4 mutant (proposed in [9]).
The models for both, wild type and mutant, include six metabolites, namely: glucose, fructose, sucrose, sugar phosphates, and starch, as well as a combined sink export, which are interconverted through seven reactions (Figure S12 in Additional file 3). These models include the most abundant sugars and sugar phosphates in the carbon metabolism of plants, without resolving subcellular compartments. Therefore, they sacrifice the complexity to increase feasibility of simulation studies. The reactions are modeled over a 24 h diurnal cycle with inclusion of limited knowledge regarding the kinetic parameters in the models of the wild type and the mutant. Therefore, these examples demonstrate that the constrainedbased formulations proposed in this study support also the inclusion of kinetic information. The kinetic models used in these comparisons depend on publicly unavailable measurements. However, among others, the rate of net photosynthesis was approximated by a smoothing spline interpolation of measurements over whole diurnal cycles. Furthermore, unknown model parameters were identified by minimizing the sum of squared errors between simulated and measured states [9]. As a result, these case scenarios of comparing the performances with the results from the kinetic models can be viewed as a proxy for the datadriven comparison.
Additional file 3. Additional results Photosynthesis Model. Results from the RSS analysis of the proposed methods for a lumped model of C_{3}plant carbohydrate metabolism are presented. The model includes 6 metabolites and 7 reactions, as depicted in Figure S12. In addition, the ODEs for the kinetic model used in establishing the RSS results are presented.
Format: PDF Size: 112KB Download file
This file can be viewed with: Adobe Acrobat Reader
The simulation day begins at 6 o'clock and has a 16 hours light and 8 hours dark phase. The time period is divided into twelve finite elements, and orthogonal collocation is applied to a time interval of two hours. This results in at least 780 variables which need to be optimized based on over 1000 constraints, depending on the different approaches. For the DFBA approach, maximal sinks production is supposed to be the systemic objective to be optimized.
RSS for metabolite concentrations and reaction fluxes
Inspection of the RSS results for the metabolite concentrations, presented in Table S1 and S3 in Additional file 3 for the wild type and the inv4 mutant, respectively, support the statements made in the previous section, dealing with the model of the CalvinBenson cyclethe NLPbased RDFBA_{CF }variant outperforms not only the three MDFBA variants, but also the classical DFBA based on maximal sink production for the wild type as well as the mutant. In addition, the same claim holds for majority of the metabolites for the basic and the fluxbased NLP variants of RDFBA.
With respect to the reaction fluxes for the wild type as well as for the inv4 mutant, the fluxbased MDFBA outperforms all other methods; nevertheless, the basic NLPbased RDFBA ranks second. Inspection of the RSS results, presented in Table S2 and S4 in Additional file 3, are in line with our observation that constraintbased approaches result in decoupling of the metabolite concentrations and reaction fluxes. Nevertheless, taken together, these findings indicate that the proposed constraintbased approaches for simulating the dynamics of metabolic networks may be a suitable tool in analyzing models of metabolic processes for which little information, aside from stoichiometry, is available.
Challenges in comparison of experimental data with the findings from method applications
The timedependent concentration and flux profiles predicted by the proposed methods are in good agreement with the respective profiles generated from the kinetic models of the CalvinBenson cycle and carbohydrate metabolism. Nevertheless, the ultimate biological validity of the minimization principle(s), used to obtain the aforementioned predications, remains to be established with the help of experimental data of in vivo metabolite concentrations and intracellular reaction fluxes. However, designing and conducting experiments from which nonstationary metabolite concentrations and reaction fluxes are to be measured is an extremely challenging task. For instance, decades of experimental work for quantifying the absolute concentrations of photosynthetic intermediates [42] and endproduct pathways have resulted in methods applicable for obtaining the in vivo steadystate levels for 40 metabolites [39]. Moreover, virtually all existing methods for determining intracellular fluxes rely on model estimates generated with the help of labeled metabolomics data and the steadystate assumption [4345]. Finally, the methods for monitoring the nonstationary state are currently not suited for systemlevel probing. Given the stateoftheart, we employed synthetic data, generated from the considered models, as a first step in validating the claims.
Application to genomescale metabolic networks
The analysis of genomescale networks is the focus of current modeling strategies with many different applications, including metabolic engineering and drug targeting [46]. Furthermore, in the last years, the quality of genomescale models has improved due to the incorporation of extensive experimental data, resulting in a large number of included reactions and metabolites. For instance, the new reconstruction of Escherichia coli includes 2,251 metabolic reactions, and 1,136 unique metabolites [47]. An accurately simulation of the behavior of this network over 10 seconds using DFBAbased methods with the same setting as for the CalvinBenson cycle model, would result in at least 56,275 unknown variables for the flux rates and 28,400 for the metabolite concentrations over time. It is important to note, that while solving some special classes of NLP can be performed in polynomial time, MINLP problems are often not tractable (i.e., they belong to the class of NPhard problems [33,48]). Therefore, an application of DFBAbased methods to genomescale models is currently hampered by the size of the resulting instances as well as the lack of optimization platforms which scale well. In addition, a comparison of the predictive power of the different methods is only possible by including experimental data, or at least kinetic model predictions, as it is the case for our analysis. To our knowledge, kinetic model based predictions, and particularly experimental data of the profile of metabolite concentrations and flux levels, are not available for genomescale models. Nevertheless, as shown above, the predictive power demonstrates that DFBAbased methods are useful for an accurate prediction of timeresolved metabolite concentrations and flux rates of metabolic networks for which kinetic parameters cannot be obtained but solving at least the NLP problems are possible.
Conclusions
In the present work, we proposed a new approach to analyze the dynamic adjustment of metabolic networks, called RDFBA. RDFBA combines the dynamic FBA with regulatory on/off minimization by minimizing the total number of significant metabolite concentration changes in comparison to flux changes in the classical approach. Additionally, we extended this method and the MDFBA approach by considering not only the fluctuation of the metabolite concentration profiles but also those of the flux levels in the network. We introduced seven new approaches and, in total, ten constraintbased approaches were implemented and their accuracy was compared with the outcome of two kinetic modelsfor the CalvinBenson cycle and for the plant carbohydrate metabolism.
For the models of the CalvinBenson cycle and the plant carbohydrate metabolism, we demonstrated that our proposed approaches yielded results in accordance with predictions from kinetic modeling, specifically for the case of metabolite concentrations. In addition, we demonstrated that the (MI)NLPbased RDFBA approach captures the dynamic coupling of reaction fluxes and metabolite concentrations. Therefore, constraintbased approaches in combination with collocation on finite elements offer a promising framework for analysis of metabolic network dynamics without specifying the details of enzyme kinetics. The proposed approaches outperform the existing variants in several cases and are suitable for positing modelbased hypotheses for the dynamics of metabolic pathways when little enzymatic details are available. Finally, our findings suggest that minimizing the combination of flux and metabolite concentration fluctuations is the mechanism most likely responsible for maintaining the metabolic network robustness due to internal perturbations.
Methods
Statistical analysis
The results of applying the described methods are compared with the outcome from kinetic modeling with the help of residual sum of squares (RSS). The RSS for each system element (i.e., metabolite and reaction flux) is defined as follows:
where y_{j }is the result of the kinetic modeling and f(x_{j}) of the compared approaches at the orthogonal root j.
The Kendall rank correlation coefficient, denoted by τ, evaluates the degree of similarity between two sets of ranks given to a same set of objects [49]. This coefficient depends upon the number of inversions of pairs of objects which would be needed to transform one rank order into the other. We use the Kendall τ to qualitatively discriminate between the constraintbased approaches with respect to their correspondence to the outcome of kinetic modeling by using the temporal distribution of both metabolic concentrations and reaction fluxes.
Implementation
All mathematical programming approaches are implemented in MATLAB 7.11.0, R2010b with the optimization platform TOMLAB v7.7 [50]. We use CPLEX to solve LP problems, the SNOPT solver for NLP problems and MINLPbb for MINLP problems. The kinetic modeling of the CalvinBenson cycle is performed by using the biochemical network simulator COPASI 4.6 [51], while the kinetic modeling of the plant carbohydrate metabolism is conducted with the Systems Biology Toolbox2 [52].
Authors' contributions
SK and ZN designed the study and developed the mathematical formulation of the methods. SK carried out the implementation and computational analysis. SK and ZN performed the statistical and comparative analyses, and were both involved in writing the manuscript. All authors read and approved the final manuscript.
Acknowledgements
Z.N. would like to acknowledge the support from the GoFORSYS project funded by the German Federal Ministry of Education and Research, Grant no. 0313924. Z.N. and S.K. thank the MaxPlanck Society for financial support.
References

Kitano H: Systems biology: a brief overview.
Science 2002, 295(5560):16621664. PubMed Abstract  Publisher Full Text

Dettmer K, Aronov PA, Hammock BD: Mass spectrometrybased metabolomics.
Mass Spectrom Rev 2007, 26:5178. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Moco S, Bino RJ, Vorst O, Verhoeven HA, de Groot J, van Beek TA, Vervoort J, de Vos CHR: A Liquid ChromatographyMass SpectrometryBased Metabolome Database for Tomato.
Plant Physiol 2006, 141:12051218. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lisec J, Schauer N, Kopka J, Willmitzer L, Fernie AR: Gas chromatography mass spectrometrybased metabolite profiling in plants.
Nat Protoc 2006, 1:387396. PubMed Abstract  Publisher Full Text

Fiehn O: Combining genomics, metabolome analysis, and biochemical modelling to understand metabolic networks.
Comp Funct Genomics 2001, 2(3):155168. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Mahadevan R, Edwards JS, Doyle JF: Dynamic Flux Balance Analysis of Diauxic Growth in Escherichia coli.
Biophys J 2002, 83:13311340. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Luo RY, Liao S, Tao GY, Li YY, Zeng S, Li XY, Luo Q: Dynamic analysis of optimality in myocardial energy metabolism under normal and ischemic conditions.
Mol Syst Biol 2006, 2:2006 0031. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Luo RY, Wei H, Ye L, Wang K, Chen F, Luo L, Liu L, Li Y, Crabbe MJC, Jin L, Li Y, Zhong Y: Photosynthetic metabolism of C3 plants shows highly cooperative regulation under changing environments: A systems biological analysis.
Proc Natl Acad Sci USA 2009, 106(3):847852. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Nägele T, Henkel S, Hörmiller I, Sauter T, Sawodny O, Ederer M, Heyer AG: Mathematical Modeling of the Central Carbohydrate Metabolism in Arabidopsis Reveals a Substantial Regulatory Influence of Vacuolar Invertase on Whole Plant Carbon Metabolism.
Plant Physiol 2010, 153:260272. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

RiosEstepa R, Lange BM: Experimental and mathematical approaches to modeling plant metabolic networks.
Phytochemistry 2007, 68(1618):23512374. PubMed Abstract  Publisher Full Text

Libourel IGL, ShacharHill Y: Metabolic Flux Analysis in Plants: From Intelligent Design to Rational Engineering.
Annu Rev Plant Biol 2008, 59:625650. PubMed Abstract  Publisher Full Text

Fischer E, Zamboni N, Sauer U: Highthroughput metabolic flux analysis based on gas chromatographymass spectrometry derived 13 C constraints.
Anal Biochem 2004, 325(2):308316. PubMed Abstract  Publisher Full Text

Varma A, Palsson BO: Metabolic Flux Balancing: Basic Concepts, Scientific and Practical Use.
Nat Biotechnol 1994, 12:994998. Publisher Full Text

Gianchandani EP, Chavali AK, Papin JA: The application of flux balance analysis in systems biology.
Wiley Interdiscip Rev Syst Biol Med 2010, 2(3):372382. PubMed Abstract  Publisher Full Text

Orth JD, Thiele I, Palsson BO: What is flux balance analysis?
Nat Biotechnol 2010, 28:245248. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Raman K, Chandra N: Flux balance analysis of biological systems: applications and challenges.
Brief Bioinform 2009, 10(4):435449. PubMed Abstract  Publisher Full Text

Steuer R, Gross T, Selbig J, Blasius B: Structural kinetic modeling of metabolic networks.
Proc Natl Acad Sci USA 2006, 103(32):1186811873. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Feist AM, Palsson BO: The biomass objective function.
Curr Opin Microbiol 2010, 13(3):344349. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Holzhütter HG: The principle of flux minimization and its application to estimate stationary fluxes in metabolic networks.
Eur J Biochem 2004, 271:29052922. PubMed Abstract  Publisher Full Text

Schuetz R, Kuepfer L, Sauer U: Systematic evaluation of objective functions for predicting intracellular fluxes in Escherichia coli.
Mol Syst Biol 2007, 3:119. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Segrè D, Vitkup D, Church GM: Analysis of optimality in natural and perturbed metabolic networks.
Proc Natl Acad Sci USA 2002, 99(23):1511215117. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Shlomi T, Berkman O, Ruppin E: Regulatory on/off minimization of metabolic flux changes after genetic perturbations.
Proc Natl Acad Sci USA 2005, 102(21):76957700. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Emmerling M, Dauner M, Ponti A, Fiaux J, Hochuli M, Szyperski T, Wuthrich K, Bailey JE, Sauer U: Metabolic Flux Responses to Pyruvate Kinase Knockout in Escherichia coli.
J Bacteriol 2002, 184:152164. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kitano H: Biological robustness.
Nat Rev Genet 2004, 5(11):826837. PubMed Abstract  Publisher Full Text

Kitano H: Towards a theory of biological robustness.
Mol Syst Biol 2007, 3:137. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Waage P, Gulberg CM: Studies concerning affinity.
J Chem Educ 1986, 63(12):1044. Publisher Full Text

Heinrich R, Schuster S: The Regulation Of Cellular Systems. 1st edition. Springer; 1996.

Grimbs S, Selbig J, Bulik S, Holzhütter HG, Steuer R: The stability and robustness of metabolic states: identifying stabilizing sites in metabolic networks.
Mol Syst Biol 2007, 3:146. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Cuthrell JE, Biegler LT: On the optimization of differentialalgebraic process systems.
AICHE J 1987, 33(8):12571270. Publisher Full Text

Wilhelm T: The smallest chemical reaction system with bistability.
BMC Syst Biol 2009, 3:90. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Eissing T, Conzelmann H, Gilles ED, Allgöwer F, Bullinger E, Scheurich P: Bistability analyses of a caspase activation model for receptorinduced apoptosis.
J Biol Chem 2004, 279(35):3689236897. PubMed Abstract  Publisher Full Text

Veening JW, Smits WK, Kuipers OP: Bistability, epigenetics, and bethedging in bacteria.
Annu Rev Microbiol 2008, 62:193210. PubMed Abstract  Publisher Full Text

Bussieck MR, Pruessner A: MixedInteger Nonlinear Programming.

Berg J, Tymoczko J, Stryer L: Biochemistry. 5th edition. New York: WH Freeman and Company; 2002.

Nelson DL, Cox MM: Lehninger Principles of Biochemistry. 4th edition. WH Freeman; 2004.

Arnold A, Nikoloski Z: A quantitative comparison of CalvinBenson cycle models.
Trends in plant science 2011, 16(12):676683. PubMed Abstract  Publisher Full Text

Zhu XG, Alba R, de Sturler E: A simple model of the Calvin cycle has only one physiologically feasible steady state under the same external conditions.
Nonlinear Anal Real World Appl 2009, 10(3):14901499. Publisher Full Text

Grimbs S, Arnold A, Koseska A, Kurths J, Selbig J, Nikoloski Z: Spatiotemporal dynamics of the Calvin cycle: Multistationarity and symmetry breaking instabilities.
Biosystems 2011, 103(2):212223. PubMed Abstract  Publisher Full Text

Arrivault S, Guenther M, Ivakov A, Feil R, Vosloh D, Dongen JTV, Sulpice R, Stitt M: Use of reversephase liquid chromatography, linked to tandem mass spectrometry, to profile the Calvin cycle and other metabolic intermediates in Arabidopsis rosettes at different carbon dioxide concentrations.
Plant J 2009, 59(5):826839. PubMed Abstract  Publisher Full Text

Collatz GJ: The interaction between photosynthesis and ribuloseP2 concentrationeffects of light, CO2 and O2.

Giersch C, Heber U, Kaiser G, Walker DA, Robinson SP: Intracellular metabolite gradients and flow of carbon during photosynthesis of leaf protoplasts.
Arch Biochem Biophys 1980, 205:246259. PubMed Abstract  Publisher Full Text

Benson A: Following the path of carbon in photosynthesis: a personal story.
Photosynth Res 2002, 73:2949. PubMed Abstract  Publisher Full Text

Nöh , Wahl A, Wiechert W: Computational tools for isotopically instationary 13 C labeling experiments under metabolic steady state conditions.
Metabolic Engineering 2006, 8(6):554577. PubMed Abstract  Publisher Full Text

Yuan J, Bennett BD, Rabinowitz JD: Kinetic flux profiling for quantitation of cellular metabolic fluxes.
Nat Protoc 2008, 3(8):13281340. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Antoniewicz MR, Kelleher JK, Stephanopoulos G: Elementary metabolite units (EMU): A novel framework for modeling isotopic distributions.
Metab Eng 2007, 9:6886. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kim TY, Sohn SB, Kim YB, Kim WJ, Lee SY: Recent advances in reconstruction and applications of genomescale metabolic models.

Orth JD, Conrad TM, Na J, Lerman JA, Nam H, Feist AM, Palsson BO: A comprehensive genomescale reconstruction of Escherichia coli metabolism2011.
Mol Syst Biol 2011, 7:535. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Murty KG, Kabadi SN: Some NPcomplete problems in quadratic and nonlinear programming.
Math Program 1987, 39(2):117129. Publisher Full Text

Abdi H: Kendall rank correlation. In Encyclopedia of Measurement and Statistics. Salkind NJ. Thousand Oaks (CA): Sage; 2007:508510.

Holmström K: The TOMLAB Optimization Environment in MATLAB. [http://tomopt.com/tomlab/] webcite

Hoops S, Sahle S, Gauges R, Lee C, Pahle J, Simus N, Singhal M, Xu L, Mendes P, Kummer U: COPASIa COmplex PAthway SImulator.
Bioinformatics 2006, 22(24):30673074. PubMed Abstract  Publisher Full Text

Schmidt H, Jirstrand M: Systems Biology Toolbox for MATLAB: A computational platform for research in Systems Biology.
Bioinformatics 2006, 22(4):514515. PubMed Abstract  Publisher Full Text