INTRODUCTION

TOP ABSTRACT INTRODUCTION DYNAMIC FLUX BALANCE ANALYSIS DIAUXIC GROWTH OF E.... RESULTS AND DISCUSSION REFERENCES

Recent developments in genomics, such as genome sequencing, microarrays, and GeneChips have provided detailed information into the genetic networks of several microorganisms (De Saizieu et al., 1998; Tao et al., 1999; Selinger et al., 2000; Oh and Liao, 2000; Wei et al., 2001). The next logical step is to use this information to study the integrated behavior of the cellular networks. One of the areas of research has been the study of metabolic networks (Oh and Liao, 2000; Tao et al., 2001; Ideker et al., 2001). The analytical and experimental methods for understanding the nature of flux distribution in a metabolic network, along with molecular biology techniques for genetic engineering, assist in the design of the metabolic reaction networks (Stephanpoulos, 1999). Mathematical analysis of metabolism can guide the metabolic engineering process; for example, Hatzimanikatis et al. (1998) have addressed the problem of determining the optimal regulatory structure in terms of gene overexpression or deletion. In that study, the regulatory structure was represented by binary variables, and the objective was to maximize a desired steady-state objective through the solution of a mixed integer linear programming formulation. Several other quantitative approaches have been proposed to study metabolic networks. These approaches include metabolic control analysis (Fell, 1996), biochemical systems theory (Savageau et al., 1987a,b), cybernetic modeling (Kompala, 1984; Dhurjati et al., 1985; Varner and Ramkrishna, 1999), and flux balance analysis (FBA) (Varma and Palsson, 1994b). With the exception of FBA, these approaches require a functional form for the kinetics of the cellular reactions.

FBA is an approach to constrain the metabolic network based on the stoichiometry of the metabolic reactions and does not require kinetic information (Varma and Palsson, 1994a). Optimization of an objective function, such as growth rate, is used to obtain a metabolic flux distribution that satisfies the constraints, and FBA has been shown to provide meaningful predictions in Escherichia coli (Varma and Palsson, 1994b; Edwards et al., 2001). van Riel and coworkers have proposed a modified FBA approach, where the flux balance analysis problem was solved along with constraints on the rate of change of metabolite levels at specific time instants (van Riel et al., 2000; Giuseppin and van Riel, 2000).

Diauxic growth represents the classical reprogramming of a metabolic network and has been extensively studied with mathematical modeling (Varma and Palsson, 1994b; Wong et al., 1997; Lendenmann and Egli, 1998; Guardia and Calvo, 2001). Ramkrishna and coworkers have also used the cybernetic modeling approach to model the diauxic growth of E. coli on mixtures of glucose and organic acids such as pyruvate, succinate, and fumarate (Ramakrishna et al., 1996; Narang et al., 1997). In cybernetic modeling, the bacterial cell is viewed as an optimal strategist that maximizes the utility of the resources provided to it. The regulation of the genes and the activity of the enzymes are obtained as a solution to an optimal resource allocation problem. Because these variables are obtained as a function of the kinetic rate equations, the result is a closed-form dynamic model of the network.

Herein, we describe dynamic flux balance analysis (DFBA), which incorporates rate of change of flux constraints. We show that DFBA can predict the dynamics of diauxic growth. Classical FBA has also been used to study diauxic growth on glucose and acetate (Varma and Palsson, 1994b). However, classical FBA incorrectly predicted the reutilization of acetate. Furthermore, classical FBA cannot predict the metabolite concentrations, which is possible with DFBA. DFBA also allows the incorporation of kinetic expression when the kinetics are well characterized. In this paper, two different formulations for DFBA are presented. These two approaches are used to analyze the diauxic growth of E. coli on glucose and acetate. The sensitivity of the approaches to the rate of change of flux constraints, the functional form of the objective function, and the parameters in the model are examined. Using this formalism, we have characterized the different phases of batch growth in terms of the active constraints during each phase. Thus, DFBA provides a significant improvement over the classical FBA and will find utility as a quantitative analysis tool in the basic sciences and biotechnology.

	DYNAMIC FLUX BALANCE ANALYSIS

TOP ABSTRACT INTRODUCTION DYNAMIC FLUX BALANCE ANALYSIS DIAUXIC GROWTH OF E.... RESULTS AND DISCUSSION REFERENCES

FBA is a modeling approach that constrains the metabolic network by the balance of the metabolic fluxes (reactions) around each node (metabolite). When the metabolic network is operating in a steady state, the mass balances are described by a set of linear equations,

(1)

where A is the m × n stoichiometric matrix of the reactions, m is the number of the metabolites, n is the number of fluxes, and v is the flux vector of the network. Because the system of linear equations is underdetermined (more unknown fluxes than equations), an objective function is used to obtain a solution using linear programming (LP). Typically, the maximization of the growth flux is used as the objective function (Varma and Palsson, 1994a; Bonarius et al., 1997; Pramanik and Keasling, 1997; Edwards et al., 2001), where the growth flux is defined in terms of the biosynthetic requirements. For details on the LP formulation, see Varma and Palsson (1994a), and Edwards et al. (1999). FBA only identifies the metabolic flux distribution, and there is no information on the metabolite concentrations or on the dynamic characteristics of the metabolic fluxes. In simulations where there is a transition between two steady states, the FBA solution will indicate an instantaneous change of the metabolic fluxes (Varma and Palsson, 1994b). Therefore, constraints on the rate of change of the fluxes must be explicitly incorporated in the problem. The dynamic extension to FBA can be formulated in the following two ways. The two formulations are discussed in detail in the following sections.

Dynamic Optimization Approach (DOA): This involves optimization over the entire time period of interest to obtain time profiles of fluxes and metabolite levels. The dynamic optimization problem was transformed to a nonlinear programming (NLP) problem and the NLP problem was solved once. The details of the objective function and the constraints in the formulation are presented in Eq. 3.

Static Optimization Approach (SOA): This approach involves dividing the batch time into several time intervals and solving the instantaneous optimization problem at the beginning of each time interval, followed by integration over the interval. The optimization problem was solved using LP repeatedly during the course of the batch to obtain the flux distribution at a particular time instant. The SOA formulation is presented in detail in Eq. 4. The objective used in the optimization problem can be similar to the objective in FBA. Varma et al. (1994b) have used FBA to obtain dynamic prediction for diauxic growth in a manner similar to this approach. However, they did not incorporate rate-of-change constraints on the metabolic fluxes.

Dynamic optimization-based DFBA approach

Consider a metabolic network with m metabolites and n fluxes. The set of conservation of mass equations, for each metabolite, results in a set of ordinary differential equations,

(2)

where z is the vector of metabolite concentrations, X is the biomass concentration, v is the vector of metabolic fluxes per gram (DW) of the biomass, A is the stoichiometric matrix of the metabolic network, µ is the growth rate obtained as a weighted sum of the reactions that synthesize the growth precursors, and w_i are the amounts of the growth precursors required per gram (DW) of biomass.

Along with the system of dynamic equations, several additional constraints must be imposed for a realistic prediction of the metabolite concentrations and the metabolic fluxes. These include non-negative metabolite and flux levels, limits on the rate of change of fluxes, and any additional nonlinear constraints on the transport fluxes. A general dynamic optimization problem can be formulated as

(3)

where z₀ and X₀ are the initial conditions for the metabolite concentration and the biomass concentration, respectively, c(v, z) is a vector function representing nonlinear constraints that could arise due to consideration of kinetic expressions for fluxes, t₀ and t_f are the initial and the final times,  is the terminal objective function that depends on the end-point concentration, L is the instantaneous objective function,  is the Dirac-delta function, t_j is the time instant at which L is considered, _ins and _end are the weights associated with the instantaneous and the terminal objective function, respectively, and v(t) is the time profile of the metabolic fluxes. If the nonlinear constraint is absent, the problem reduces to an optimization involving a bilinear system.

The dynamic optimization problem was solved by parameterizing the dynamic equations through the use of orthogonal collocation on finite elements (Cuthrell and Biegler, 1987). The time period (t₀-t_f) was divided into a finite number of intervals (finite elements). The fluxes, the metabolite levels, and the biomass concentration were parameterized at the roots of an orthogonal polynomial within each finite element. The details of the parameterization for a specific example are presented in the next section. Continuity of the metabolite and the biomass concentrations was imposed at the beginning of each of the finite elements. The time derivative of the variables was approximated as a linear combination of the value of the fluxes at each point, and the dynamic equations were transformed to algebraic equations. The nonlinear constraints were imposed at discrete points in the time interval considered. Thus the dynamic optimization was converted to an NLP problem. The resulting NLP was solved using the fmincon function in MATLAB (The MathWorks Inc., Natwick, MA).

Static optimization-based DFBA approach

In SOA, the time period was divided into N intervals. In the absence of the nonlinear constraints involving the fluxes, the optimization problem is reduced to an LP problem. The LP was solved at the beginning of each interval to obtain the fluxes at that time instant:

(4)

where T is the length of the time interval chosen.

The optimization problem was solved using CPLEX. The dynamic equations were integrated assuming that the fluxes were constant over the interval. The optimization problem was then formulated at the next time instant and solved. This procedure was repeated from t₀ to t_f. For the class of systems involving only bilinear terms with fluxes and the biomass concentration, it is possible to directly solve the dynamic equations and thereby eliminate the numerical integration.

	DIAUXIC GROWTH OF E. COLI ON GLUCOSE AND ACETATE

TOP ABSTRACT INTRODUCTION DYNAMIC FLUX BALANCE ANALYSIS DIAUXIC GROWTH OF E.... RESULTS AND DISCUSSION REFERENCES

The metabolic network considered for modeling the diauxic growth of E. coli is shown in Fig. 1. From a metabolic pathway analysis with glucose, acetate, and oxygen as the input and biomass and acetate as the output, a set of ~300 extreme pathways were identified (Schilling et al., 2000a). The biomass composition and the ratio of precursors required were obtained from the literature (Schilling et al., 2000b). From this set, four pathways were chosen based on the biomass yield and the known physiology of E. coli (Cronan and Laporte, 1996; Oh and Liao, 2000) to define a simplified metabolic network (see Figs. 2 and 3). The extreme pathways chosen represented both aerobic and anaerobic utilization of glucose and had the highest biomass yield from among the 300 pathways. The acetate utilization pathway was chosen to be consistent with experimental observations that the pckA gene coding for the PEP carboxykinase is expressed during growth on acetate (Oh and Liao, 2000). The simplified network was then used in all further studies presented in the paper.

View larger version (21K): [in this window] [in a new window]   FIGURE 1   Metabolic network of E. coli considered for the FBA. The network consisted of 54 metabolites and 85 reactions. Glycolysis, pentose phosphate pathway, TCA cycle with the glyoxylate bypass, anapleurotic reactions, and redox metabolism are included in the metabolic network.

View larger version (11K): [in this window] [in a new window]   FIGURE 2   Simplified metabolic network. The network identified after pathway analysis with glucose, acetate, and oxygen as the input and biomass as the output and selection based on biomass yield is presented above.

View larger version (30K): [in this window] [in a new window]   FIGURE 3   The metabolic pathways used to simplify the network v1 (top left), v2 (top right), v3 (bottom left), and v4 (bottom right). The details of the pathways in the simplified network are shown above. The active reactions are highlighted reactions in the pathways.

A dynamic model for the prediction of the time profiles for a batch bioreactor with glucose as the carbon source is presented in the equations,

(5)

where A^Glcxt, A^Ac, A^O₂ are the rows of the stoichiometric matrix associated with glucose, acetate, and oxygen, respectively, k_La is the mass transfer coefficient for oxygen and is assumed to be 7.5 hr¹ (Edwards et al., 2001).

The key variables in the mathematical model of the metabolic network are the glucose concentration, the acetate concentration, the biomass concentration, and the oxygen concentration in the gas phase. The oxygen concentration in the gas phase was assumed to be a constant (0.21 mM). A term for the oxygen transport from the gas phase (air at ambient temperature) was included in the model. The oxygen transport rate was assumed to be directly proportional to the difference in concentration. The oxygen uptake rate was constrained to allow a maximum possible flux of 15 mmol/gdw hr (Varma and Palsson, 1994b). Transport of acetate across the cell was assumed to be rapid (with respect to the metabolic flux); therefore, the internal and the external concentrations were assumed to be the same. The glucose uptake rate was bounded by Michaelis-Menten kinetics involving the glucose concentration (Wong et al., 1997). The DFBA formulation for the analysis of diauxic growth in E. coli is presented in the next subsection.

DFBA: DOA formulation

The DOA formalism of DFBA was used to analyze the diauxic growth of E. coli. The objective function for the DOA formalism is detailed in the equations,

Case 1:  Instantaneous objective

(6a)

Case 2:  Terminal time objective

(6b)

where N_s is the number of collocation points for the parameterization of the metabolite and biomass concentrations; Z^st  R^4N_s is the stacked vector containing the metabolite and biomass concentrations in time; K_m is the saturation constant (0.015 mM, Wong et al., 1997); z₀ is a vector consisting of the initial glucose, acetate, and oxygen concentrations; N_v is the number of collocation points of the fluxes; V^st  R^4N_v is the stacked vector containing the fluxes in time; _max is the rate of change of flux constraints imposed; C⁰ is the matrix containing the derivative weights; f(Z^st, V^st) is the function containing the derivative vector along with the continuity condition (determined from Eq. 5); and µ^sc is the growth rate determined from the initial and final biomass concentration measurements used in scaling the objective function.

For the DOA formalism, each time interval was divided into five finite elements, and the variables were parameterized at the roots of the fifth-order Legendre polynomial, resulting in 204 variables. The flux rates-of-change constraints were included in the optimization problem as linear constraints. The NLP was solved for two different objective functions involving the biomass concentration, and the results are presented in the next section. The first objective function (instantaneous objective (Eq. 6a) involved maximizing the scaled sum of the biomass concentration at the collocation points. As the biomass concentration increases 1000-fold during the course of the batch, the concentrations at different time points were scaled, so that all the time points are equally weighted. The second objective function (terminal time objective (Eq. 6b) maximized the biomass concentration at the final time.

DFBA: SOA formulation

For DFBA using SOA, the time of the batch (10 hrs) was divided into 10,000 intervals, and the optimization was formulated as described in the Eq. 4 and was solved using CPLEX. The number of variables in the optimization problem was four (corresponding to the number of the fluxes). The optimization was solved at the beginning of each interval, and the metabolite concentrations at the beginning of the next interval were found by direct integration.

The parameters used for the DFBA were the maximal oxygen and the glucose uptake rates (Varma and Palsson, 1994b), the mass transfer coefficient (Edwards et al., 2001), the substrate saturation constant (Wong et al., 1997), and the flux rate-of-change constraints. The only parameter that could not be identified based on the existing literature was the flux rate-of-change constraints. These parameters, however, can be estimated from biochemical parameters such as the transcription and translation rates and genomic information involving regulatory elements, microarray data, and proteomics (Tavazoie et al., 1999; Cohen et al., 2000; Kirkpatrick et al., 2001). Thus, in the case where the transcription and translation rates are known, the rate of change of flux constraints can be identified precisely. For the current study, a range of values for the rate of change of fluxes provided reasonable agreement between the model predictions and the experimentally observed time domain data. A single parameter set within the range was chosen for the present study.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION DYNAMIC FLUX BALANCE ANALYSIS DIAUXIC GROWTH OF E.... RESULTS AND DISCUSSION REFERENCES

The DFBA approaches were used to simulate batch growth of E. coli on glucose, where acetate is initially secreted and subsequently utilized. The data from a batch fermentation (Varma and Palsson, 1994b) is also plotted in all the figures.

Static optimization-based approach: Results

The results from the DFBA using the SOA are shown in Figs. 4 and 5. In Fig. 4, the flux rate-of-change constraints were relaxed for the purpose of comparison. The DFBA solution was used to identify the constraints governing cellular growth. It was determined that different constraints were active during different times in the batch culture. We defined distinct phases of the fermentation based on differences in the active constraint. It was observed that, up to 4.6 hr, the constraints on the oxygen and glucose uptake rates were limiting growth and were the active constraints. In the next phase of the fermentation (from 4.6 to 6.9 hr) the oxygen concentration in the fermentation environment approached zero, and the system was constrained by the transport of oxygen (governed by k_La term). At 6.9 hr, the glucose was nearly completely consumed, and, from this point until the glucose concentration reached zero, the system was limited by the glucose (Michaelis-Menten kinetics) and oxygen uptake rate constraints. When the glucose concentration was zero, the acetate utilization began, and the growth was characterized by the oxygen transport limitation, which was influenced by the k_La term. The growth in the final phase (on acetate) was linear, and not exponential as in the previous phase, due to the k_La constraint.

View larger version (28K): [in this window] [in a new window]   FIGURE 4   Model prediction using the SOA for DFBA in the absence of the rate of change of flux constraints. Interpretation of the constraints governing the growth of E. coli in the three phases is shown above. In the first phase, the constraints are the oxygen and the glucose uptake rate. The transport of oxygen along with the glucose uptake constrained the growth in the middle phase. Growth on acetate in the final phase was again constrained by oxygen transport. Glucose, acetate, and biomass concentrations from experimental data are plotted along with the model predictions (Varma and Palsson, 1994b).

View larger version (28K): [in this window] [in a new window]   FIGURE 5   Model prediction using SOA for DFBA in the presence of the rate of change of flux constraints. The constraints governing the growth are similar to the previous figure except for the region where the growth is constrained by the rate of change of flux constraints, and pathway 3 is active. Glucose, acetate, and biomass concentrations from experimental data are plotted along with the model predictions (Varma and Palsson, 1994b).

The flux rate-of-change constraints were also imposed on the metabolic network, and the simulations produced similar results (Fig. 5), with the exception of additional phases that were governed by the flux rate of change constraints. The flux rate-of-change constraints were active from 5.5 to 6.5 hr, where the flux from pathway 3 that produced both biomass and acetate was nonzero. This was due to the constraint on the flux rate of change of the pathway that produced acetate in the absence of oxygen (pathway 4).

Sensitivity to the oxygen uptake rate

View larger version (27K): [in this window] [in a new window]   FIGURE 6   Initial flux distribution as a function of the oxygen uptake rate for the SOA to DFBA. When the oxygen uptake rate is not sufficient to support aerobic growth (pathway 2), then the anaerobic pathway (v4) becomes active.

Sensitivity to the glucose uptake rate

View larger version (27K): [in this window] [in a new window]   FIGURE 7   Model prediction using SOA for DFBA in the presence of the rate of change of flux constraints for a glucose uptake rate of 11 mmol/gdwhr. Insufficient oxygen uptake rate due to the increased glucose uptake results in the shutting down of the acetate utilization pathway in the initial phase. Glucose, acetate, and biomass concentrations from experimental data are plotted along with the model predictions (Varma and Palsson, 1994b).

Sensitivity to the mass transfer coefficient (k_La)

1

1

View larger version (26K): [in this window] [in a new window]   FIGURE 8   Model prediction using SOA for DFBA in the presence of the rate of change of flux constraints for the case where kLa = 12.5 hr1. Final phase involving acetate utilization is constrained by the transport of oxygen. Increased oxygen availability results in higher rate of acetate utilization. Glucose, acetate, and biomass concentrations from experimental data are plotted along with the model predictions (Varma and Palsson, 1994b).

Dynamic optimization based approach: Results

The results from the DOA for DFBA are presented in Fig. 9. The rate of change of flux constraints were imposed at all time instants, unlike SOA (where the constraints were relaxed whenever the concentrations were close to zero). Therefore, for this case, the time evolution was marginally slower than the case of dynamic FBA using SOA for the same parameter set. However, when the flux rate of change constraints were relaxed, time evolution was rapid (Fig. 10). Sensitivity studies similar to the previous approach (SOA) were performed, and the results of the simulations for this approach (DOA) were similar. This is to be expected because these two approaches were formulated to produce the same results. The differences in the two approaches are related to the flexibility in problem formulation and the computational requirements (see Discussion).

View larger version (26K): [in this window] [in a new window]   FIGURE 9   Model prediction using DOA for DFBA in the presence of the rate of change of flux constraints. The results shown above are similar to the earlier results obtained using SOA. Glucose, acetate, and biomass concentrations from experimental data are plotted along with the model predictions (Varma and Palsson, 1994b).

View larger version (26K): [in this window] [in a new window]   FIGURE 10   Model prediction using DOA for DFBA in the absence of the rate of change of flux constraints. Glucose, acetate, and biomass concentrations from experimental data are plotted along with the model predictions (Varma and Palsson, 1994b).

Sensitivity to the objective function

View larger version (26K): [in this window] [in a new window]   FIGURE 11   Model prediction using DOA for DFBA where the objective is maximizing the end-point biomass concentration. The results obtained for this objective function do not match the experimental data. The biomass concentration achieved is higher than the previous case, and pathway 2 is active until the end of the batch. Glucose, acetate, and biomass concentrations from experimental data are plotted along with the model predictions (Varma and Palsson, 1994b).

Discussion

We have extended the classical FBA for analyzing the dynamic reprogramming of a metabolic network. In particular, we have examined the reprogramming of the metabolic network that occurs at different stages of diauxic growth of E. coli on glucose. Two approaches for DFBA were introduced, and the sensitivity to the different parameters was analyzed. The results were compared to the data presented in Varma and Palsson (1994b).

DFBA using SOA extended the FBA approach presented in Varma and Palsson (1994b) through the incorporation of the flux rate-of-change constraints. In this paper, the model for diauxic growth of E. coli considered the effect of oxygen transport, and the metabolic network studied was simplified using pathway analysis to obtain a compact representation. The scope of the results obtained for modeling the metabolic reprogramming during diauxic growth presented here were similar to those based on FBA. Cybernetic models have also been proposed for the study of diauxic growth (Ramakrishna et al., 1996; Narang et al., 1997). The fluxes in the cybernetic approach are obtained as a solution to an optimal resource allocation problem with an instantaneous objective function. Typically, only a subset of the network is considered in the optimization problem (Varner, 2000). The solution obtained is analytic, and one can represent the system with a dynamic model. However, the cybernetic approach requires kinetic information for all the reactions in the network. DFBA does not require kinetic information and considers the entire network, although the solution for the fluxes is not analytic and is obtained by solving an optimization problem.

DFBA using DOA allows the formulation of a dynamic objective function describing characteristics, such as, reduction of transition time between metabolic states (Torres, 1994) or end-point biomass optimization, into a rigorous mathematical framework. A dynamic objective function based on the desired goal could provide information useful in the design of genetically modified metabolic networks for metabolic engineering by taking into account the dynamic responses to fluctuations in the system. The static optimization-based DFBA would not allow such a dynamic formulation, because the optimization performed is at a specific time instant. However, in SOA, the number of variables that have to be solved is far fewer (in each optimization) in comparison, and the optimization problem is an LP problem as opposed to the NLP for DOA. As the size of the network increases, the number of variables and the number of constraints would increase proportionally in the NLP. Thus, SOA is scalable to larger metabolic networks.

DFBA provides a framework for modeling the dynamic responses of a metabolic network to various perturbations. In this paper, we have examined the applicability of this framework for modeling the diauxic growth in E. coli. The results from DFBA are qualitatively similar to the experimental observations. DFBA was able to predict the onset of acetate production and also the preference of E. coli for sequential utilization of acetate and glucose over the simultaneous utilization. The constraints governing the behavior were identified at various phases in the batch culture. It was found that, in the initial phase, the glucose and oxygen uptake rates were the active constraints. In the middle phase, the oxygen concentration is close to zero, and the mass transfer coefficient (k_La) and the maximum allowed rate of change of flux was found to constrain the system. Acetate utilization (last phase) was found to be constrained completely by the oxygen mass transfer coefficient.

The sensitivity to the various parameters was studied, and it was found that the dynamic model was most sensitive to k_La, whereas it was less sensitive to other parameters. The importance of the objective function was examined, and it was found that an instantaneous objective function was more representative of the experimental results than an end-point objective function. Another advantage of dFBA is that can incorporate kinetic expressions for reactions that are well-studied. This approach could also be used to identify regulatory phenomena and obtain insight into the functioning of the metabolic pathways. Changes in the regulatory structure that optimize the dynamics of a particular metabolic process could be obtained as a solution to a modified DFBA problem.

In conclusion, we have presented analysis tools for the quantitative study of the dynamic reprogramming of metabolic networks. These tools, along with experimental technologies such as microarrays, GeneChips, and proteomics, will help further understanding of the dynamic behavior of metabolic networks. Additionally, the DFBA approach can be used to provide strategies for the design of a network with a desired objective for metabolic engineering. Finally, the DFBA approach is an extension to classical FBA and has demonstrated great potential; however, further analysis is needed to improve the predictive capabilities in the biological sciences.