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 time-resolved metabolic state, characterized by both fluxes and concentrations. By conducting a comparative analysis with a kinetic model of the Calvin-Benson 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 FBA-based 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 [2–4], allowing reconstruction and analysis of genome-scale 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 [6–10].

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, 13–16]. By optimizing an objective function, the linear program identifies one feasible flux distribution from the set of fluxes satisfying the constraints imposed by the mass-balance 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 steady-state 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., Michaelis-Menten or Hill) have to be determined by accurate measurements of enzyme activities and data-fitting to experimentally obtained (time-course) 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 time-resolved metabolic profiles with limited knowledge of enzyme kinetics [6–8]. Unlike, the analyses based on FBA, which focus on the steady-state behavior, DFBA offers the means to analyze transient (non-steady) states. In addition, DFBA has been combined with MOMA, resulting in the M-DFBA 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]. M-DFBA 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 DFBA-based approaches. For instance, the posited mechanism in M-DFBA 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 knock-outs, 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 ROOM-based DFBA approach (abbreviated to R-DFBA) by coupling the principle of ROOM with DFBA. R-DFBA 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 R-DFBA and the existing M-DFBA 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 M-DFBA [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 FBA-based approach in simulating the dynamics of metabolic networks must be established based on the comparison of the obtained results and the outcome of a well-defined kinetic model for a given metabolic network. Here, by comparing the results of applying DFBA, M-DFBA, R-DFBA and their extensions to those of the kinetic models for the Calvin-Benson 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 R-DFBA outperform the existing DFBA-based approaches.

Results and discussion

Here we describe the seven proposed methods and the comparison of their performance with two kinetic models--of the Calvin-Benson 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.

We point out that DFBA-based methods and kinetic modeling constitute independent approaches. The DFBA-based 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 enzyme-kinetics 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 DFBA-based 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 network-the steady-state assumption. Mahadevan and coworkers introduced two DFBA formulations--static 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 time-resolved 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.

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.

M-DFBA

Luo and coworkers developed an approach called M-DFBA 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 , where N is the number of metabolites in the network and x_{
i,j
}represents the concentration of metabolite i at the time point given by the orthogonal root j. Hence, the general M-DFBA optimization problem is defined as follows:

with M representing the number of orthogonal roots, and δ, the Dirac delta function (see Additional file 1).

Extended versions of M-DFBA

We extend the M-DFBA 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 M-DFBA_{
CF
}, extends the M-DFBA 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 enzyme-kinetic 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 M-DFBA 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 time-resolved metabolic state, we consider a third type of the M-DFBA method, whereby the objective function, subject to minimization, includes only the fluctuation of flux rates, as follows:

R-DFBA

In contrast to MOMA, ROOM minimizes the total number of significant (large enough) flux changes from the wild-type flux distribution. In our approach, we combine ROOM with DFBA, naming this modeling approach R-DFBA. 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 M-DFBA approach. Unlike M-DFBA, in R-DFBA this is realized by minimizing the significant concentration changes at adjacent orthogonal roots. In contrast to M-DFBA, which captures smooth changes over time (by the Euclidean distance), R-DFBA 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 points--the 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 R-DFBA 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 R-DFBA, may be necessary to capture the difference in concentrations between two steady states in a bi-stable system [30–32].

Two different types of the R-DFBA 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 mixed-integer linear program (MILP) in the basic ROOM approach [22]. The first R-DFBA 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 mixed-integer nonlinear programming problem (MINLP). These problems are difficult to solve precisely as they couple the combinatorial nature of mixed-integer programming (MIP) problem and the computational complexity of solving NLP problems [33].

Consequently, we formulate a second R-DFBA approach based on relaxed integer constraints, described as follows:

Compared to the first formulation of R-DFBA, 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 R-DFBA

As for the M-DFBA approach, we also extend the proposed R-DFBA 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 R-DFBA_{
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, R-DFBA_{
CF
}is solved based on mixed-integer nonlinear programming, by the first R-DFBA 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 M-DFBA_{
CF
}. In comparison, in R-DFBA_{
CF
}large concentration and flux rate changes between two adjacent orthogonal roots for some metabolites and/or reactions are possible, which is precluded in M-DFBA_{
CF
}.

Additionally, we consider a third objective function for the R-DFBA 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 R-DFBA, flux-based R-DFBA and R-DFBA_{
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 Calvin-Benson cycle

The Calvin-Benson 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 ribulose-1,5-bisphosphate (RuBP) to get two molecules of phosphoglycerate (PGA), (2) reduction, converting the obtained PGA into 1,3-diphosphoglycerate (DPGA) and glyceraldehyde-3-phosphate (GAP), and (3) regeneration, which recovers RuBP after several intermediate steps from ribulose-5-phosphate (Ru5P) [34, 35]. The Calvin-Benson cycle is a well-studied 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.

Table 1

Reactions in the simplified model of the Calvin cycle

Symbol Enzyme name

Biochemical reaction

Simplified reaction

v_{
max
}

v_{1}

RuBisCo

RuBP + CO_{2} → 2 PGA

RuBP → 2 PGA

3.78

v_{2}

PGA kinase

PGA + ATP → DPGA + ADP

PGA → DPGA

11.75

v_{3}

GAP dehydrogenase

DPGA + NADPH → GAP + P + NADP

DPGA → GAP

5.04

v_{4}

GAP → 0.6 Ru5P

GAP → 0.6 Ru5P

3.05

v_{5}

Ru5P kinase

Ru5P + ATP → RuBP + ADP

Ru5P → RuBP

8

v_{6}

Sink capacity

PGA → Sink

PGA → Sink

3

v_{7}

Sink capacity

GAP → Sink

GAP → Sink

0.1

With the simplified model of the Calvin cycle, we examine and compare the performance of the proposed R-DFBA approach together with its extensions, the existing M-DFBA approach and our proposed extensions, as well as the classical DFBA. In total, ten different approaches, based on different constraint-based 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 ribulose-5-phosphate (Ru5P), ribulose-1,5-bisphosphate (RuBP), 3-phosphoglycerate (PGA), 1,3-diphosphoglycerate (DPGA), glyceraldehyde-3-phosphate (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 M-DFBA and R-DFBA.

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 DFBA-based 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 R-DFBA (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 constraint-based and kinetic modeling approaches. Furthermore, we determine the temporal evolution of Kendall τ for the metabolic states of the constraint-based 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.

R-DFBA yields accurate time-resolved 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 R-DFBA approach, minimizing the concentration fluctuations, the results differ between the variants: while the NLP-based variant predicts that the concentration of DPGA are almost constant over time, the MINLP-based formulation results in changes of the concentration of DPGA over the 10-second time interval. The results of basic M-DFBA show that the concentration of DPGA slightly decreases during the first three seconds, and is almost constant following this time period.

The effect of minimizing the flux changes is especially pronounced in the two variants of the flux-based R-DFBA. Here, the results for the DPGA concentration are similar to those of the basic M-DFBA and its extensions. The behavior of the DPGA concentration predicted by both flux-based R-DFBA 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 (NLP-based) R-DFBA_{
CF
}, which combines the basic and flux-based variant of R-DFBA. In the variants of R-DFBA_{
CF
}, the changes in flux rates and metabolite concentrations are weighted equally, as given by the objective function: . A different weighting in the objective function of R-DFBA_{
CF
}could avoid the loss of the observed behavior. This suggest that the objective function may have a more intricate form (e.g., , where α is a parameter weighting the contribution of reaction rates and metabolite concentrations).

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 M-DFBA and its extensions are nearly constant over time; on the other hand, with the proposed NLP-based R-DFBA 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 MINLP-based R-DFBA variants. We point out that experimental measurements in isolated cells and chloroplasts have resulted in low levels of RuBP [39–41].

Comparing these results with the outcome of kinetic modeling shows that for these two metabolites MINLP-based R-DFBA predicts the most accurate temporal profiles. The constant concentration in the second half of the simulated interval of the NLP-based R-DFBA is very close to that of kinetic modeling, for both, the basic R-DFBA and R-DFBA_{
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 R-DFBA not only on MINLP but also on NLP.

Table 2

Residual sum of squares values for each metabolite

Method

Metabolite

RuBP

PGA

DPGA

GAP

Ru5P

Sink

DFBA

20.07

2.21

16.15

2.55

3.55

224.86

M-DFBA

28.02

7.58

13.33

2.54

3.20

567.42

Flux-based M-DFBA

36.66

2.32

17.28

2.22

3.57

534.20

M-DFBA_{
CF
}

15.45

28.17

18.98

2.15

3.54

564.26

R-DFBA NLP

1.35

76.20

8.93

3.91

22.24

582.23

Flux-based R-DFBA NLP

14.38

345.25

13.68

2.35

5.66

558.03

R-DFBA_{
CF
}NLP

1.86

14.97

7.28

3.57

23.53

582.23

R-DFBA MINLP

3.80

67.25

11.72

3.17

0.61

129.49

Flux-based R-DFBA MINLP

18.17

0.66

9.76

2.50

18.85

497.71

R-DFBA_{
CF
}MINLP

9.80

5.17

5.84

3.01

55.48

420.66

The discrepancy between the kinetic modeling results and the results obtained from the DFBA-based approaches is quantified based on the residual sum of squares

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 kinetic-based model. For Ru5P, the results obtained by the basic MINLP-based R-DFBA are in accordance with the kinetic modeling results and outperform the remaining approaches; on the other hand, the NLP-based R-DFBA variants predict a constant amount of Ru5P of about 1 mmol l^{-1} except for the flux-based variant, which also predict a nearly consumed Ru5P after four seconds similar to the kinetic-based results. Surprisingly, the predictions based on the MINLP-based R-DFBA_{
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 R-DFBA_{
CF
}variants is due to the over-constraining of both fluxes and concentrations.

In addition, the PGA concentration closest to the kinetic modeling results is predicted by the MINLP-based R-DFBA approach which minimizes significant flux changes. Contradictory to the results from all other methods, the flux-based R-DFBA, 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 NLP-based variants of M-DFBA and R-DFBA. However, for all MINLP-based R-DFBA 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 MINLP-based R-DFBA approach predicts the best simulation results for the concentration of sink.

Taken together, the results of the RSS-based analysis demonstrate that MINLP-based R-DFBA outperforms M-DFBA in four out of six cases in predicting metabolite concentrations over time.

(De)coupling of time-resolved flux predictions

The flux levels over time, also known as time-resolved 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 NLP-based variants of both M-DFBA and R-DFBA differ from the outcome of kinetic modeling. The reactions rates are very small for these six methods. Comparing the basic methods with the flux-based variants of M-DFBA and R-DFBA show smoothing of curves, which is also appears in the results of M-DFBA_{
CF
}and R-DFBA_{
CF
}. Furthermore, it can be observed, that after four seconds for all NLP-based R-DFBA 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 non-perturbed 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 DFBA-based methods.

Table 3

Residual sum of squares values for each reaction

Method

Reaction

v_{
1
}

v_{
2
}

v_{
3
}

v_{
4
}

v_{
5
}

v_{
6
}

v_{
7
}

DFBA

5.93

57.44

7.23

6.47

38.49

19.14

0.07

M-DFBA

64.73

175.32

154.96

149.77

76.70

16.70

0.06

Flux-based M-DFBA

58.60

142.74

127.75

131.11

67.76

16.47

0.03

M-DFBA_{
CF
}

53.91

143.06

125.41

126.05

67.09

16.85

0.06

R-DFBA NLP

41.13

130.15

113.79

111.48

53.72

16.85

0.06

Flux-based R-DFBA NLP

54.21

187.37

173.83

164.65

82.47

16.30

0.04

R-DFBA_{
CF
}NLP

51.68

144.16

143.96

145.26

70.02

16.85

0.06

R-DFBA MINLP

63.95

169.37

158.54

158.91

73.16

17.78

0.06

Flux-based R-DFBA MINLP

47.99

108.16

99.58

102.80

61.20

16.61

0.05

R-DFBA_{
CF
}MINLP

25.91

61.38

59.78

60.24

44.04

14.36

0.02

The discrepancy between the kinetic modeling results and the results obtained from the different DFBA-based approaches is quantified based on the residual sum of squares

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 DFBA-based extensions. Indeed, we observe that the predicted rates obtained by the MINLP-based R-DFBA_{
CF
}are closer to the kinetic results than all M-DFBA approaches, as quantified by the RSS.

In the results of the MINLP-based R-DFBA 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 higher-order 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 flux-based M-DFBA and three of the proposed approaches--the flux-based R-DFBA and R-DFBA_{
CF
}based on MINLP as well as M-DFBA_{
CF
}.

In addition, these results demonstrate that while the NLP-based R-DFBA 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 M-DFBA-based approaches, except the flux-based M-DFBA. These findings indicate that some results of DFBA-based 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 MINLP-based formulations, especially the flux based R-DFBA and R-DFBA_{
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 constrained-based 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 constrained-based 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 data-driven comparison.

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 Calvin-Benson cycle-the NLP-based R-DFBA_{
CF
}variant outperforms not only the three M-DFBA 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 flux-based NLP variants of R-DFBA.

With respect to the reaction fluxes for the wild type as well as for the inv4 mutant, the flux-based M-DFBA outperforms all other methods; nevertheless, the basic NLP-based R-DFBA ranks second. Inspection of the RSS results, presented in Table S2 and S4 in Additional file 3, are in line with our observation that constraint-based approaches result in decoupling of the metabolite concentrations and reaction fluxes. Nevertheless, taken together, these findings indicate that the proposed constraint-based 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 time-dependent concentration and flux profiles predicted by the proposed methods are in good agreement with the respective profiles generated from the kinetic models of the Calvin-Benson 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 non-stationary 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 end-product pathways have resulted in methods applicable for obtaining the in vivo steady-state 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 steady-state assumption [43–45]. Finally, the methods for monitoring the non-stationary state are currently not suited for system-level probing. Given the state-of-the-art, we employed synthetic data, generated from the considered models, as a first step in validating the claims.

Application to genome-scale metabolic networks

The analysis of genome-scale 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 genome-scale 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 DFBA-based methods with the same setting as for the Calvin-Benson 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 NP-hard problems [33, 48]). Therefore, an application of DFBA-based methods to genome-scale 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 genome-scale models. Nevertheless, as shown above, the predictive power demonstrates that DFBA-based methods are useful for an accurate prediction of time-resolved 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 R-DFBA. R-DFBA 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 M-DFBA 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 constraint-based approaches were implemented and their accuracy was compared with the outcome of two kinetic models--for the Calvin-Benson cycle and for the plant carbohydrate metabolism.

For the models of the Calvin-Benson 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)NLP-based R-DFBA approach captures the dynamic coupling of reaction fluxes and metabolite concentrations. Therefore, constraint-based 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 model-based 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 constraint-based 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 Calvin-Benson 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].

Declarations

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 Max-Planck Society for financial support.

Authors’ Affiliations

(1)

Max-Planck Institute of Molecular Plant Physiology

(2)

Institute of Biochemistry and Biology, University of Potsdam

References

Kitano H: Systems biology: a brief overview.Science 2002,295(5560):1662–1664.PubMedView Article

Moco S, Bino RJ, Vorst O, Verhoeven HA, de Groot J, van Beek TA, Vervoort J, de Vos CHR: A Liquid Chromatography-Mass Spectrometry-Based Metabolome Database for Tomato.Plant Physiol 2006, 141:1205–1218.PubMedView Article

Lisec J, Schauer N, Kopka J, Willmitzer L, Fernie AR: Gas chromatography mass spectrometry-based metabolite profiling in plants.Nat Protoc 2006, 1:387–396.PubMedView Article

Fiehn O: Combining genomics, metabolome analysis, and biochemical modelling to understand metabolic networks.Comp Funct Genomics 2001,2(3):155–168.PubMedView Article

Mahadevan R, Edwards JS, Doyle JF: Dynamic Flux Balance Analysis of Diauxic Growth in Escherichia coli.Biophys J 2002, 83:1331–1340.PubMedView Article

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.PubMedView Article

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):847–852.PubMedView Article

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:260–272.PubMedView Article

Rios-Estepa R, Lange BM: Experimental and mathematical approaches to modeling plant metabolic networks.Phytochemistry 2007,68(16–18):2351–2374.PubMedView Article

Libourel IGL, Shachar-Hill Y: Metabolic Flux Analysis in Plants: From Intelligent Design to Rational Engineering.Annu Rev Plant Biol 2008, 59:625–650.PubMedView Article

Fischer E, Zamboni N, Sauer U: High-throughput metabolic flux analysis based on gas chromatography-mass spectrometry derived 13 C constraints.Anal Biochem 2004,325(2):308–316.PubMedView Article

Varma A, Palsson BO: Metabolic Flux Balancing: Basic Concepts, Scientific and Practical Use.Nat Biotechnol 1994, 12:994–998.View Article

Gianchandani EP, Chavali AK, Papin JA: The application of flux balance analysis in systems biology.Wiley Interdiscip Rev Syst Biol Med 2010,2(3):372–382.PubMedView Article

Orth JD, Thiele I, Palsson BO: What is flux balance analysis?Nat Biotechnol 2010, 28:245–248.PubMedView Article

Raman K, Chandra N: Flux balance analysis of biological systems: applications and challenges.Brief Bioinform 2009,10(4):435–449.PubMedView Article

Steuer R, Gross T, Selbig J, Blasius B: Structural kinetic modeling of metabolic networks.Proc Natl Acad Sci USA 2006,103(32):11868–11873.PubMedView Article

Holzhütter HG: The principle of flux minimization and its application to estimate stationary fluxes in metabolic networks.Eur J Biochem 2004, 271:2905–2922.PubMedView Article

Schuetz R, Kuepfer L, Sauer U: Systematic evaluation of objective functions for predicting intracellular fluxes in Escherichia coli.Mol Syst Biol 2007, 3:119.PubMedView Article

Segrè D, Vitkup D, Church GM: Analysis of optimality in natural and perturbed metabolic networks.Proc Natl Acad Sci USA 2002,99(23):15112–15117.PubMedView Article

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):7695–7700.PubMedView Article

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:152–164.PubMedView Article

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 Calvin-Benson cycle models.Trends in plant science 2011,16(12):676–683.PubMedView Article

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):1490–1499.View Article

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):212–223.PubMedView Article

Arrivault S, Guenther M, Ivakov A, Feil R, Vosloh D, Dongen JTV, Sulpice R, Stitt M: Use of reverse-phase 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):826–839.PubMedView Article

Collatz GJ: The interaction between photosynthesis and ribulose-P2 concentration-effects of light, CO2 and O2.Year B Carnegie Inst Wash 1978, 77:248–251.

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:246–259.PubMedView Article

Benson A: Following the path of carbon in photosynthesis: a personal story.Photosynth Res 2002, 73:29–49.PubMedView Article

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):554–577.PubMedView Article

Yuan J, Bennett BD, Rabinowitz JD: Kinetic flux profiling for quantitation of cellular metabolic fluxes.Nat Protoc 2008,3(8):1328–1340.PubMedView Article

Antoniewicz MR, Kelleher JK, Stephanopoulos G: Elementary metabolite units (EMU): A novel framework for modeling isotopic distributions.Metab Eng 2007, 9:68–86.PubMedView Article

Kim TY, Sohn SB, Kim YB, Kim WJ, Lee SY: Recent advances in reconstruction and applications of genome-scale metabolic models.Curr Opin Biotechnol 2011.

Orth JD, Conrad TM, Na J, Lerman JA, Nam H, Feist AM, Palsson BO: A comprehensive genome-scale reconstruction of Escherichia coli metabolism--2011.Mol Syst Biol 2011, 7:535.PubMedView Article

Murty KG, Kabadi SN: Some NP-complete problems in quadratic and nonlinear programming.Math Program 1987,39(2):117–129.View Article

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

Holmström K: The TOMLAB Optimization Environment in MATLAB.
[http://tomopt.com/tomlab/]
Advanced Modeling and Optimization 1999, 1:47–69.

Hoops S, Sahle S, Gauges R, Lee C, Pahle J, Simus N, Singhal M, Xu L, Mendes P, Kummer U: COPASI-a COmplex PAthway SImulator.Bioinformatics 2006,22(24):3067–3074.PubMedView Article

Schmidt H, Jirstrand M: Systems Biology Toolbox for MATLAB: A computational platform for research in Systems Biology.Bioinformatics 2006,22(4):514–515.PubMedView Article