Analysis of Growth of Lactobacillus plantarum WCFS1 on a Complex Medium Using a Genome-scale Metabolic Model*

A genome-scale metabolic model of the lactic acid bacterium Lactobacillus plantarum WCFS1 was constructed based on genomic content and experimental data. The complete model includes 721 genes, 643 reactions, and 531 metabolites. Different stoichiometric modeling techniques were used for interpretation of complex fermentation data, as L. plantarum is adapted to nutrient-rich environments and only grows in media supplemented with vitamins and amino acids. (i) Based on experimental input and output fluxes, maximal ATP production was estimated and related to growth rate. (ii) Optimization of ATP production further identified amino acid catabolic pathways that were not previously associated with free-energy metabolism. (iii) Genome-scale elementary flux mode analysis identified 28 potential futile cycles. (iv) Flux variability analysis supplemented the elementary mode analysis in identifying parallel pathways, e.g. pathways with identical end products but different co-factor usage. Strongly increased flexibility in the metabolic network was observed when strict coupling between catabolic ATP production and anabolic consumption was relaxed. These results illustrate how a genome-scale metabolic model and associated constraint-based modeling techniques can be used to analyze the physiology of growth on a complex medium rather than a minimal salts medium. However, optimization of biomass formation using the Flux Balance Analysis approach, reported to successfully predict growth rate and by product formation in Escherichia coli and Saccharomyces cerevisiae, predicted too high biomass yields that were incompatible with the observed lactate production. The reason is that this approach assumes optimal efficiency of substrate to biomass conversion, and can therefore not predict the metabolically inefficient lactate formation.

Genome-scale metabolic models have become available for an increasing number of organisms (1). These models link genomic information to metabolic reaction networks and have gained considerable attention, not only within the context of metabolic engineering (2,3), but also from the bioinformatics field (4,5). An increasing repertoire of modeling techniques is available for exploring the capabilities of the metabolic network, and for understanding genotype-phenotype relationships (2,6). These techniques are often referred to as constraint-based modeling techniques.
So far, constraint-based modeling techniques have mainly been applied to microorganisms that grow on a minimal salt medium containing a single carbon source (1). However, in many biological niches, as well as in multicellular organisms, cells encounter more complex nutritional environments. In analogy, in many industrial and biotechnological applications of microorganisms and cell lines, complex media are used, either because the cells are auxotrophic for nutrients or such media are cheaper. Multiple inputs for the metabolic network complicate constraint-based modeling approaches considerably. Moreover, one may expect quite different metabolic behavior in such a complex environment compared with that in minimal salts media. The question arises whether and to what extent constraint-based modeling is useful in such biotechnologically relevant cases.
We used constraint-based modeling methods to explore the metabolism of Lactobacillus plantarum, a lactic acid bacterium (LAB) 2 used in a variety of industrial food fermentations, and marketed as a probiotic (6). L. plantarum grows only in rich media, such as decomposed plant materials and meat, and consequently has become auxotrophic for many vitamins and amino acids. Hence, the chemically defined medium used in our laboratory contains glucose, citrate, acetate, nucleosides and nucleobases, vitamins, and 18 amino acids (7). We have recently performed a metabolic reconstruction, based on a careful annotation of the genomic content (8), combined with single omission experiments to infer functional pathways (7). Here we report on the construction and analysis of a genomescale metabolic model of L. plantarum, and the comparison with experimental fermentation data.
We illustrate how different stoichiometric modeling approaches provide a detailed analysis of the metabolism in this lactic acid bacterium, and how it generates testable hypotheses. Unlike the apparent success to describe and predict metabolism in Escherichia coli and Saccharomyces cerevisiae, however, optimization of biomass formation using flux balance analysis failed to predict metabolism in L. plantarum. We will discuss this result in the light of microbial physiology.

EXPERIMENTAL PROCEDURES
Fermentation and Growth Conditions-L. plantarum WCFS1 was grown at 37°C in chemically defined medium (7) in a 10-ml tube, and used as inoculum of 500 ml of pH-controlled (pH 5.5) cultures in a 1.0-liter fermentor (Applikon Biotechnology BV, the Netherlands). The headspace was flushed with nitrogen at a stirring speed of 120 rpm. At A 600 ϭ 1, the medium pump was switched on at dilution rates 0.1-0.5 h Ϫ1 , and steady state was achieved within 5 volume changes (as judged by constant A 600 ). At steady state, 150-ml samples were taken and immediately spun down in a Megafuge 1.0R at 4°C (Heraus Instruments, Germany). Cell pellet was used to measure biomass composition as detailed in supplementary Materials I.
Supernatant was used to measure compound concentrations. Organic compounds were measured by high performance liquid chromatography as described previously (9). Free amino acids were determined by Ansynth BV (Roosendaal, The Netherlands). At the time of harvest, samples from the feed stock were also taken and analyzed for organic compounds and amino acids, as described above. Fluxes q i (in mmol h Ϫ1 gDW Ϫ1 ) were calculated as: Construction of a Genome-scale Model and Constraint-based Modeling-The genome-scale metabolic model was based on a metabolic reconstruction that was described in detail elsewhere (7). The model was developed within the Simpheny TM software environment (Genomatica Inc., San Diego, CA). Also the constraint-based modeling techniques, i.e. flux balance analyses (2,10) and flux variability analysis (11), were carried out within Simpheny TM . FBA and FVA are based on linear programming as described (10,11). FBA finds a flux distribution that optimizes an objective function; FVA finds the flux value range (called the "span") of each reaction by minimizing and maximizing the flux through that reaction, given a set of constraints. In Ref. 11 the constraints were set by a FBA solution; here we used experimentally derived flux values (of uptake rates and secretion rates) to constrain the solution space.
Blocked Reactions and Unbalanced Metabolites-To identify all blocked reactions that could not carry any flux (12), FVA was carried out: blocked reactions are those reactions that have a minimal and maximal flux value of zero. Note that the results of FVA and hence, the number of blocked reactions, depend on the constraints of the system. Unbalanced metabolites are those metabolites (row-vectors in the stoichiometry matrix) that have only zero entries after deletion of the blocked reactions (columns) in the stoichiometry matrix.
Elementary Flux Mode Analysis-For elementary flux mode (EFM) analysis, the reactions and metabolites from the Simpheny TM output files where parsed via Python scripts into a Meta-Tool (13) input file. Only ATP, ADP, phosphate, water and protons were set as external metabolites, i.e. as metabolites that can act as source or sink in an EFM. Because all reactions in the network are elementary balanced, only circulations and futile cycles can form EFMs under these conditions. Determination of Energy Parameters for Maintenance and Growth-To determine how much ATP is needed for biomass formation and maintenance, it is needed to estimate the total ATP production rate in the system at different growth rates. The slope of a plot of ATP production rate against growth rate gives the amount of ATP needed to form biomass Y ATP ; the intercept is the maintenance requirement m ATP , i.e. the ATP production rate at zero growth rate.
It is important to note that the biomass equation used to simulate growth in the final model is: biomass components ϩ y ATP ϩ y H 2 O 3 biomass ϩ y ADP ϩ y P i ϩ y H. Y ATP , the total amount of ATP needed for growth, is actually the sum of the ATP needed to form the biomass components, and the amount of ATP to assemble the biomass components into biomass (y). The amount of ATP needed to make biomass components (defined as x), is explicitly accounted for by the model reactions and will be calculated below. We first seek to find y, which can be done with a genome-scale model by setting y ϭ 0 in the biomass equation and instead introduce an ATP dissipation reaction, v ATP (ATP ϩ H 2 O 3 ADP ϩ P i ϩ H). We use FBA to find the maximal v ATP . The problem is defined as: maximize v ATP subject to, where S is the stoichiometry matrix, v is the flux vector, SD standard deviation, is the growth rate (with y ϭ 0), and D is the dilution rate of the chemostat. The result is the maximal rate of ATP production that the system can attain at the measured fluxes and at the growth rate. In steady state, under the assumption of full coupling, ATP production equals ATP consumption, so that v ATP ϭ m ATP ϩ y⅐. Hence, in a plot of v ATP against the growth rate, y is the slope and m ATP is the intercept (Fig. 1). To be able to compare total ATP production rates estimated by the traditional approach (see main text) and the genome-scale approach we compare: traditional, V total,trad ϭ v lactate ϩ 2v acetate and genome-scale, V total,GS ϭ (x ϩ y)⅐ ϩ m ATP . Where v lactate and v acetate are the rates of lactate and acetate production, respectively. To compute V total,GS we still need to estimate the amount of ATP required for precursor biosynthesis (x). For this, we can use the reduced cost of the biomass flux for ATP dissipation. The reduced cost r i of a flux v i is a sensitivity coefficient generated by the linear programming algorithm, indicating how much the objective function would change if the boundary constraint of that flux would be changed: r i ϭ dv ATP /dv i . Thus, for the biomass formation flux , its reduced cost for ATP dissipation indicates how much less ATP can be dissipated when biomass formation is increased. That reduced cost (given in Table  3) therefore represents x.
Flux Variability and Modeling of Uncoupling-A flux variability analysis is formulated: @v i ⑀v, minimize and maximize v i subject to Equation 2.
In the section of energy parameters, it was assumed in steady state that all ATP that is being produced by catabolism is completely used for growth. Hence, coupling is modeled as v ATP ϭ max(v ATP ), where again y ϭ 0 in the biomass equation. In the uncoupled case, the amount of ATP needed for maintenance and assembly is left open, i.e. v ATP was allowed to vary between 0 and max(v ATP ).

RESULTS
Genome-scale Metabolic Model: Features-We constructed a genome-scale model on the basis of the metabolic reconstruction reported previously (7). Such a model defines all geneprotein reaction associations, which can be complex (6,7). The reactions form a metabolic network that needs to be carefully curated for the purpose of constraint-based modeling (see Refs. 1 and 14 for reviews on this process). In a metabolic reconstruction that aims merely at putting genes in a metabolic context, such as pathways, uncertainties and omissions are not essential; for a quantitative model, however, essential aspects are (i) definition of a quantitative biomass equation, involving macromolecular content and composition (see supplementary Material I), (ii) analysis and closing of network gaps, i.e. reactions in essential metabolic routes to which no gene could be assigned, (iii) choosing specific stoichiometry of reactions, such as cofactor usage and transport mechanisms, (iv) determination of ATP coefficients quantifying growth-and non-growth associated ATP consumption processes not explicitly accounted for in the model (1, 6) (see below). Given the different perspectives of each network, it is difficult to directly compare the reconstruction of Ref. 7 with the current model. An indication of the extent of curation is the fact that the current model has 120 reactions that are not associated to genes (see below).
The final model consisted of 721 genes (23.5% of the total number of genes) and 643 balanced reactions of which 120 (18.5%) were not associated with a gene but inferred from biochemical and physiological data (7). The model had 531 balanced metabolites, 97 of which had an extracellular equivalent (see Table 1 for further details). When the genome is taken as the starting point in metabolic reconstruction, reactions may be included in the model that do not have a source and sink (e.g. isolated reactions). Such reactions cannot carry any flux at steady state. The number of reactions that participate in these so-called blocked reactions (12), and the degrees of freedom of the network depend on the available compounds that can be exchanged with the environment (see Table 1).
The complete model is available in supplementary Material IV. An important feature of the model is the balancing of protons. Proton gradients are essential for the overall energetics of bacteria (for review on this in LAB, see Ref. 15). Another feature of our model is the use of vitamins and cofactor pools in the biomass equation: new cells form the only sink for these com-pounds, and all reactions associated to the biosynthesis of these compounds would carry no flux without inclusion in the biomass equation. This feature was only recently introduced also in the S. cerevisae model, leading to improved predictability of lethality (16).
Maintenance and Growth-associated Energy Coefficients-Essential parameters in a stoichiometric model are the energy requirement for maintenance, and for growth (17). The growthassociated energy requirement lumps (partly unknown) energy requirements for biomass precursor biosynthesis and for assembly of these biomass components into biomass. ATP parameters determine how much biomass can be made from the available amount of free energy. Assuming that the catabolic ATP production rate is balanced by the anabolic ATP consumption rate under energy limitation, the ATP parameters can be estimated from ATP production rates, calculated from physiological data, at different growth rates (18). Therefore, L. plantarum was grown (anaerobically) in glucose-limited chemostats at different dilution rates, and product formation and substrate consumption rates were measured. Table 2 summarizes the results for the organic compounds; amino acid fluxes can be found in supplementary Material II. Because we found no differences in relative product formation and biomass composition over the dilution rates tested, data are presented as the average yield over the dilution range. Of all carbon taken up, 67 Ϯ 1% (g/g) was converted into lactate ( Table 2). For many amino acids, uptake rates were larger than required for biomass production (see supplementary Material II for details).
Traditionally, anaerobic ATP production is estimated from the lactic acid (1 ATP per lactate) and acetic acid (2 ATP per acetate) fluxes (19). When the ATP production was calculated in this way, and plotted against the growth rate, a linear relationship was found (black circles, Fig. 1). Such an approach, however, causes potential errors in the case of multiple input fluxes (i.e. citrate and amino acids). For example, catabolism of  (7)) for chemically defined medium, and additionally glutathione and menaquinone for the complete network. b Calculated as number of balanced reactions Ϫ (number of balanced metabolites Ϫ number of linear dependencies).
serine can contribute to the lactic acid flux without generating ATP. The genome-scale model can be used to circumvent this problem: the measured fluxes were set as constraints and FBA was used to calculate the maximal amount of ATP that could be generated. The resulting ATP production v ATP reflects the amount of ATP that the system requires (i.e. can generate) for growth and maintenance, additional to the ATP required for precursor biosynthesis (see "Experimental Procedures" for details). When the maximal ATP production rate v ATP is plotted against growth rate (white circles, Fig. 1), the intercept gives the maintenance coefficient (m ATP ϭ 0.36 mmol h Ϫ1 gDW Ϫ1 ), and the slope gives the amount of ATP required for assembly (y ϭ 27.4 mmol gDW Ϫ1 ). These numbers are similar to those for S. cerevisiae (20) and L. lactis (21). To compare the traditional approach for total ATP production estimation, based on lactate and acetate, with that using the genome-scale model, we calculated the total amount of ATP generated by the model under the experimental constraints. This is not trivial in the case of a complex medium with so many different inputs (see "Experimental Procedures"). When the traditional approach for total ATP production estimation was compared with that using the genome-scale model (gray circles, Fig. 1), virtually the same result was obtained. Apparently amino acid and citrate metabolism led to many small effects that compensated each other.
The genome-scale model approach did have an important additional benefit, however. Through its sensitivity analysis, the optimization algorithm provided an immediate insight into which fluxes contributed to the total ATP production, and which would negatively impact on them. In Table 3, exchange fluxes are shown with non-zero scaled reduced costs. The scaled reduced cost, R i , was calculated as:

TABLE 2 Amount of substrates consumed and products formed per g of biomass in continuous fermentations on chemically defined medium (CDM) with low (25 mM) and high (100 mM) glucose concentrations
Data presented are yield data averaged over the dilution rates indicated. For obtaining flux constraints for simulation, these data can be multiplied by the growth rate to obtain fluxes (in mmol h Ϫ1 gDW Ϫ1 ).    Where r i is the unscaled reduced cost dv ATP /dv i . Presented in this way, R i resembles a response coefficient as known from metabolic control analysis, which quantifies the relative effect of a change in a parameter to the variable of the system (22, 23). The scaled reduced cost thus gives a better indication of how important a step is for the objective value. In this case, very small fluxes may stoichiometrically appear to contribute a lot to ATP generation (such as nucleotide metabolism, Table 3), but a doubling of the corresponding flux will quantitatively have a small impact. Not unexpectedly, biomass had the largest negative effect on ATP dissipation flux. Other important contributors to energy metabolism were glucose, citrate, more acetate production (instead of lactate), and uptake of many amino acids.

Uptake or production CDM (25 mM glucose) CDM (100 mM glucose)
Catabolic pathways of amino acids were identified that were already known to contribute to ATP production through proton translocation, thereby reducing the amount of ATP needed for proton balancing (15). However, the analysis also indicated that catabolism, via transamination, of aromatic and branchedchain amino acids can generate ATP. This is very relevant because transamination and subsequent catabolism of these amino acids is a major contributor to flavor formation in fermented food products (24), but it was not previously implicated in ATP production. Closer inspection (detailed in supplementary Materials II) showed that transamination and degradation of these amino acids results in a proton-motive force-driven transhydrogenase reaction, converting NADH into NADPH. Despite the cost in proton-motive force, the overall transhydrogenase reaction has two beneficial effects: (i) NADPH is generated and relieves NADPH production via the pentose phosphate pathway that produces a lower ATP yield from glucose than does glycolysis; and (ii) NADH is reoxidized without a concomitant production of a reduced fermentation product such as lactate or ethanol and hence, relatively more acetate can be produced. Note that L. plantarum does not have a transhydrogenase enzyme of its own.
Analysis of the head space through purge-and-trap and subsequent gas chromatography-mass spectrometry analysis showed that catabolic products of the relevant amino acids are produced during fermentation (see supplementary Materials II), indicating that these catabolic pathways are indeed active. Thus, genome-scale analysis of fermentation data led to a new potential role for anaerobic degradation of some amino acids in L. plantarum.
Uncoupling, Flux Variability Analysis, and Genome-scale Detection of Futile Cycles-The estimation of ATP parameters was based on the assumption that all generated ATP is fully coupled for biosynthesis and growth. This may be reasonable for cells under energy limitation, but it is not always correct. Uncoupling between ATP production and anabolic ATP demand is commonly observed under conditions of energy excess (25), and is in fact the basis of many applications such as dough leavening and lactic acid production.
To test the impact of uncoupling on the metabolic capacities of the model, flux variability analysis (11) was used to calculate for each reaction in the model the range of flux values that were compatible with the experimental constraints (called the span). In this analysis, we compared two situations: (i) full coupling, where all ATP generated by the model is used for growth and maintenance; (ii) complete uncoupling, where the amount of ATP needed for growth and maintenance is completely free to vary (see "Experimental Procedures" for details). Fig. 2 shows the span of reactions, normalized to the maximal ATP production rate max(v ATP ). In the uncoupled case (blue in Fig. 2), much more internal flexibility in the reaction rates was observed compared with the coupled case (green in Fig. 2), as indicated by a higher span. Some reactions had an infinite span, corresponding to cycles with no net conversion reaction. We call these cycles circulations, defined in graph theory as paths in which all nodes have zero convergence, i.e. are completely balanced and cannot act as sink or source (26). These circulations cannot operate due to thermodynamic constraints (27).
Many other reactions had a normalized span of exactly 1.5, 1, 0.5, 0.2, and 0.167, indicating that these reactions were involved in futile cycles with 2/3, 1, 2, 5, and 6 ATP molecules hydrolyzed per cycle, respectively. To systematically identify all circulations and futile cycles from the genome-scale model, we used EFM analysis (28). By setting only ATP, ADP, phosphate, water, and protons as external metabolites, i.e. as the sources and sinks of the network, a combinatorial explosion (29) was prevented and 9 circulations and 28 futile cycles were identified (see supplementary Material III for details).
Reactions involved in a futile cycle all had an ATP-dependent span. However, there were also many reactions that were not involved in futile cycles, yet showed an ATP-dependent span (non-zero normalized span). These reactions form alternative pathways in the network that have different energetic consequences. An example is serine metabolism: if both catabolic and biosynthetic pathways were active at the same time, they would constitute a by-pass of glycolysis from 3-phosphoglycerate to pyruvate that does not yield ATP at the level of pyruvate kinase (Fig. 3). Input and output flux data from a chemostat at 0.11 h Ϫ1 and with a feed of 25 mM glucose were used to set the constraints. The green bars indicate the span of reactions under the assumption of full coupling between catabolic ATP production and anabolic ATP consumption (ATP demand for growth and maintenance was fixed at the maximal ATP production value). The blue bars indicate the span of each reaction when the ATP demand flux was allowed to vary, i.e. the assumption of full coupling was relaxed. The span was normalized to the maximal ATP production rate. Reactions with a normalized span Ͼ0.01 are shown, sorted by magnitude in the uncoupled case (blue).
Flux Balance Analysis Fails to Predict Homolactic Fermentation in L. plantarum-To test the use of the model further, we repeated the experiments but increased the concentration of glucose in the feed from 25 to 100 mM. Again we observed mainly lactate production and excess consumption of essentially the same set of amino acids as found with 25 mM glucose in the feed (Table 2 and supplementary Material II). A decrease in biomass yield on glucose was observed and could be attributed to increased maintenance requirements. The increase in maintenance is probably caused by lactate stress, the reactor concentration of which was 40.1 Ϯ 0.2 and 160.9 Ϯ 9.0 mmol liter Ϫ1 for 25 and 100 mM glucose, respectively (averages over dilution rates). This was inferred from a consistently higher v ATP estimated by FBA (on the basis of measured fluxes at 100 mM glucose, see "Experimental Procedures") compared with expected rates (i.e. m ATP ϩ y⅐ based on the 25 mM glucose experiment; difference of 1.23 Ϯ 0.31 mmol ATP h Ϫ1 gDW Ϫ1 over all dilution rates).
When FBA was applied to determine which fluxes limited biomass formation and which fluxes limited ATP production, the fluxes with non-zero reduced costs overlapped. This indicated that growth, unexpectedly, was still energy-limited even at the high glucose concentration. Indeed, addition of glucose led to enhanced biomass formation (results not shown). These analyses illustrate, again, the added value of a genome-scale model, and FBA, in interpreting fermentation data.
FBA, however, has also been used to predict growth rates and fermentation product formation based on evolutionary arguments (20,30,31). We therefore also tried to use FBA to calculate the optimal biomass and product formation rates, taking as constraints the nutrient feed rates (and the proper ATP parameters). Note that up until now we used measured fluxes of products as constraints. When the product formation fluxes were not constrained to observed fluxes, FBA predicted higher biomass formation rates by an average factor of 1.4 (for the 25 mM glucose experiment). Concomitantly FBA predicted exclusively mixed acid fermentation, i.e. acetate production with formate and ethanol as necessary by-products for NADH reoxidation. No optimal solution was compatible with lactate formation.

DISCUSSION
In this study we have presented a genome-scale metabolic model of L. plantarum. We have taken much effort to ensure high quality of the reconstructed metabolic network and the corresponding model. Presence and absence of functional pathways was tested experimentally as presented before (7). Many components of the biomass were measured, under the same experimental conditions under which parameters concerning the energetics of growth have been experimentally determined. Finally, we have measured as many input and output fluxes as possible, again, under the same experimental conditions. The work presented thus forms a coherent set of physiological data that is of high value on its own.
We demonstrated the added value of a genome-scale model in analyzing the data set. We illustrated this by estimating energy parameters for growth and maintenance in a complex medium. We also showed how optimization through linear programming can give clues of which components in the medium can limit certain objectives, such as ATP production or biomass formation. In this way, we could predict and experimentally validate that growth in a chemostat was still energy limited with a chemically defined medium with 100 mM glucose in the feed.
Furthermore, opportunities for discovering previously unknown metabolic capabilities of the reconstructed network have been illustrated by the degradation of amino acids, but also by the identification of futile cycles and parallel pathways. The relevance of identified transhydrogenase activity hidden in the valine catabolic pathway is uncertain, but one may predict that excess catabolism of amino acids through transaminase reactions may be stimulated by NADPH demand, e.g. only during growth. In fact, in the model, catabolism of the measured excess of amino acid uptake could completely replace the need for the NAPDH generating part of the pentose phosphate pathway. We could not find any study that specifically addressed the relation between flavor formation and growth rate.
We performed a genome-scale elementary flux mode analysis, by carefully defining the external metabolites, i.e. the sources and sinks in the network that need not be balanced. Earlier methods to overcome the combinatorial explosion that is associated with EFM analysis also adjusted the status of metabolites between internal and external, but more based on topological properties (29,32). The method can be extended to other processes, e.g. it will also work for identifying all EFMs that result in transhydrogenase activity, by setting only NAD(P) and NAD(P)H as external metabolites.
We applied flux variability analysis with experimentally derived flux constraints. The result of the FVA analysis shows which fluxes in the large network have been well determined by the measurements (small flux span), and which fluxes cannot be estimated (large flux span). Additional flux constraints, for example, obtained by 13 C labeling studies (33) could then be used to estimate those uncertain fluxes further. Here we used FVA to identify flexible parts in the network that were dependent on ATP availability (Fig. 2).
The futile cycles and parallel pathways that were identified by the methods discussed above are very relevant for an important physiological phenomenon: uncoupling between catabolic ATP production and anabolic ATP consumption, especially under energy excess. The identified flexibility in the metabolic network of L. plantarum for ATP production and consumption may help in explaining a largely unresolved question in the case of metabolic uncoupling: where does all the ATP go (18,25)? Isotopic labeling studies have failed to identify a clear primary futile cycle that can consume all the excess ATP generated by catabolism under glucose excess (34), although high fluxes through certain futile cycles have been observed under specific (often, paradoxically, glucose-limited) conditions, in particular involving pyruvate metabolism (35)(36)(37), and proton (38) or potassium leakage currents over the cell membrane (39). Our analyses identified many more futile cycles that could potentially contribute to ATP dissipation. We have tried to find indications whether and which futile cycles may be used specifically in L. plantarum under glucose excess, by means of microarray analysis, but failed to get a clear answer. 3 This indicates, not unexpectedly, that these cycles are regulated on the metabolic level, e.g. through allosteric activation or inhibition. Alternatively, ATP may simply not be generated by catabolic pathways. In E. coli the methyl glyoxal pathway has been described as a means to catabolize excess glucose without ATP generation (25), and we identified other alternatives for L. plantarum such as serine metabolism (Fig. 3).
These analyses, illustrating the use of a genome-scale model in data analysis, are important not only for fermentations involving LAB, but also for many other biotechnological and (bio)medical applications that either use complex growth or production media, or use media with excess sugars. The flexibilities and high number of futile cycles that we observed are by no means specific for L. plantarum or LAB: most reconstructed metabolic networks are characterized by high degrees of freedom (1).
Prediction of fermentative behavior based on optimization of biomass formation via the FBA approach failed for L. plantarum. Lactic acid production was not compatible with optimized growth, and the optimized growth yield was much higher than observed. Because FBA optimizes for biomass formation at a given glucose uptake rate, and therefore in essence predicts optimal yields, not rates, the FBA result can be rationalized from the higher ATP yield for mixed acid fermentation (3 per glucose) than for homolactic fermentation (2 per glucose). Thus, L. plantarum, even at energy limitations, uses a catabolic route that is less efficient in ATP production, and this cannot be predicted by FBA.
This failure of FBA to accurately predict fluxes and growth yields has been observed before. In L. lactis a similar problem occurred (21), which was treated by constraining the capacity of pyruvate-formate lyase, the committed step for mixed acid fermentation, to an experimentally determined value. Similarly, ethanol formation in S. cerevisiae was only predicted when the oxygen consumption rate was restricted to the observed value (20). Thus, energetically inefficient metabolic behavior, referred to as overflow metabolism and observed in many microorganisms (34,40,41), has been modeled within the constraint-based modeling paradigm as, indeed, overflow from the optimal pathway that has reached its maximal capacity.
It is important to realize that when fluxes are constrained by additional observations, the model becomes more and more descriptive, and less and less predictive. Moreover, the resulting fit to experimental results masks the fact that, from an evolutionary perspective, the underlying reasoning is imprecise. Constraining fluxes and optimizing for biomass yield assumes that "being efficient" is the main driving force for evolution, and that cells exhibiting overflow metabolism are still thriving for efficiency. A high efficiency, however, is a strategy that may result in good fitness, but it is not an objective on its own. There are many other strategies toward fitness.
Growth rate comes close to harnessing fitness, at least under prolonged cultivation (42) (although robustness is clearly also important and seems to come at the cost of growth efficiency and rate (43)). The thermodynamic principle that the (bio)chemical rate is driven by free-energy dissipation (44), underlies recent results that fast but inefficient metabolic strategies may be winning when microorganisms compete for common substrate (45,46). Alternative and more complex strategies for fitness are also relevant, also for microorganisms (47). For LAB the role of lactate as a pH-lowering and hence cytotoxic agent needs to be taken into account as a means to safeguard a niche in a nutrient-rich environment.
Irrespective of which strategy is actually being used, overflow metabolism should be viewed, at least in LAB, as a switch in strategy away from efficiency, not as "just" overflow from the optimal pathway that has reached its maximal capacity. The fact that pyruvate formate lyase in L. lactis is repressed by glucose and allosterically inhibited by glycolytic intermediates (48,49), strongly indicates such a switch in strategy. Also studies on adaptive evolution in E. coli have shown that some winning strategies showed inefficient (overflow) metabolism, depending on the substrate (50). Thus, one may expect efficiency to be the winning strategy only under certain nutrient-poor conditions, such as batch growth on acetate (30) or glycerol (51). These are the conditions where FBA, with biomass yield as an objective function, has been predictive. Under other conditions, such as studied here, FBA has limited predictive value.
In conclusion, we have illustrated and discussed the use of a genome-scale metabolic model for analysis of fermentative behavior of L. plantarum. This study can be seen as a model case for many biotechnological and biomedical situations where one is confronted with complex and nutrient-rich conditions. Under such conditions, the strategy of efficiency is probably less relevant for fitness than strategies aiming at robustness, fast growth rate, or preventing others to grow. Because constraint-based modeling techniques are in essence dealing with yields, extensions to the constraint-based modeling repertoire that can capture these other types of objectives would be much needed to enhance the predictive capabilities of genome-scale metabolic models.