Metabolic modeling of a mutualistic microbial community

The rate of production of methane in many environments depends upon mutualistic interactions between sulfate-reducing bacteria and methanogens. To enhance our understanding of these relationships, we took advantage of the fully sequenced genomes of Desulfovibrio vulgaris and Methanococcus maripaludis to produce and analyze the first multispecies stoichiometric metabolic model. Model results were compared to data on growth of the co-culture on lactate in the absence of sulfate. The model accurately predicted several ecologically relevant characteristics, including the flux of metabolites and the ratio of D. vulgaris to M. maripaludis cells during growth. In addition, the model and our data suggested that it was possible to eliminate formate as an interspecies electron shuttle, but hydrogen transfer was essential for syntrophic growth. Our work demonstrated that reconstructed metabolic networks and stoichiometric models can serve not only to predict metabolic fluxes and growth phenotypes of single organisms, but also to capture growth parameters and community composition of simple bacterial communities.


Introduction
By interacting with other species, organisms are able to occupy environmental niches that they otherwise could not. It is the complex web of interactions between species that defines the structure of communities and how those communities respond to environmental change. One form of interspecies interaction that may significantly influence community structure and stability is a mutualism, in which two or more species provide a net benefit to one another (Hay et al, 2004). There are numerous examples of mutualisms of different degrees of coupling in nature. For instance, mutualisms between tubeworms and methane-oxidizing or sulfide-oxidizing bacteria support life in hydrocarbon seeps in the ocean (Brooks et al, 1987;Cavanaugh et al, 1987); coral-zooxanthella mutualisms drive the growth of corals (Hay et al, 2004); a mutualistic interaction between Archaea and Bacteria is suggested to be responsible for the anaerobic oxidation of methane, an important component of global methane cycles (Boetius et al, 2000). Mutualisms may be especially important in microbial communities where multiple species are involved in degrading organic substrates (Schink, 2002;Schink and Stams, 2002;Stams et al, 2006). However, the full extent of the interaction in any given microbial mutualism is difficult to discern. Behavioral interactions between individual microbes cannot be observed directly as they can with macro-organisms. Instead, the manner in which products are exchanged between species must be inferred biochemically and genetically. In the post-genomic era, one possibility is to develop a mutualism model composed of strains that have sequenced genomes and are more experimentally tractable than natural communities. With this approach, it is possible to test diverse hypotheses by manipulating the metabolism of each species and employing advanced high-throughput technologies for genome-wide expression profiling, proteomics, and metabolic modeling. In this work, we explored the practicality of this approach by developing and evaluating a metabolic model for an ecologically relevant microbial mutualism, a syntrophic association between Desulfovibrio vulgaris Hildenborough and Methanococcus maripaludis S2.
A 'syntrophy' is one form of microbial mutualism that is commonly involved in the degradation of organic substrates by microbial communities. In syntrophic interactions, the transfer of metabolites between species is essential for growth (Schink, 1997(Schink, , 2002Schink and Stams, 2002). An example of such interactions is the widely distributed phenomenon of interspecies hydrogen transfer between sulfate reducers (or related syntrophs, e.g., Syntrophobacter spp) and methanogens, which use hydrogen to gain energy by reducing carbon dioxide into methane (Bryant et al, 1977;McInerney and Bryant, 1981;Traore et al, 1983;Pak and Bartha, 1998). The paradigm for the syntrophic association examined in this study is that in the absence of a suitable electron acceptor for the sulfate reducer (e.g., sulfate), methanogens create favorable thermodynamic conditions by scavenging hydrogen and keeping its partial pressure low, allowing the sulfate reducer to grow by 'fermenting' the carbon source (Schink, 1997). It is thought that this type of interaction is widely distributed and plays a pivotal role in many natural and engineered systems that are dependent upon microbial processes (Figure 1).
Numerous studies have examined the syntrophic growth of sulfate reducers and methanogens, yet it is still unclear exactly what is being exchanged in these interactions, how energy is conserved, and whether the metabolite(s) exchanged vary with changing ecological conditions. For example, it was suggested from theoretical considerations that formate may also act as electron carrier in such syntrophies (Boone et al, 1989). Although the same energy could be obtained from utilization of formate or hydrogen, these molecules have different chemical properties that may affect their relative functions across environments or across species with diverse physiological properties (Boone et al, 1989;Schink, 1997). Unfortunately, the carrier of electron exchange will always be at a low concentration, making it difficult to answer these questions with present analytical technology (Schink, 1997). To develop a predictive mechanistic framework for identifying key features of the extended metabolic network of syntrophic growth, we constructed a flux-balance model that captures the core metabolism of the two species. We compared model predictions for growth in monoculture and co-culture with experimental results.
Flux-balance models for several organisms have been used to define the range of phenotypes theoretically possible for a given genotype by simulating environmental conditions or turning off fluxes through particular pathways (Edwards and Palsson, 2000a;Klamt et al, 2002;Van Dien and Lidstrom, 2002;Borodina et al, 2005;Feist et al, 2006). We report here on the first construction of a flux-balance model for a twoorganism system and the testing of its performance by comparing its predictions to measurements of various growth parameters. In our model, lactate was supplied as the sole carbon and energy source for growth of the co-culture and for Desulfovibrio growing in monoculture on sulfate. The model consists of a series of linear equations, each expressing a relationship between metabolites, and the rate of flux through the reaction (Edwards et al, 2002). The reactions that make up the model were based on curated annotation of the two genome sequences. Because of the number of unknowns in the model, fluxes were calculated using linear optimization of a specific reaction, usually biomass production (Price et al, 2004).
We limited the size of the model for each organism to reactions involved in central metabolism for initial testing of a two-species model. First, we performed simulations with the sub-models for each species to test the feasibility of various physiological scenarios. We examined the accuracy of these models by comparing simulations against known growth characteristics of that organism in monoculture. These The digestive system of ruminants and other animals, fresh water sediments, and anaerobic digesters for the treatment of agricultural waste are representative environments where syntrophy-based methanogenesis is vigorously active, likely as a result of syntrophic associations between methanogenic archaea and hydrogen-producing microorganisms such as D. vulgaris. (B) Schematic representation of an anaerobic trophic cascade. Hydrogen produced in the fermentation of biopolymers and fatty acids serve as the electron donor for the secondary fermenters and subsequently for methanogenic reduction of carbon dioxide. (C) FISH microphotograph of a syntrophic association between M. maripaludis (red) and D. vulgaris (green), the model syntrophy for this study. (D) Diagram representing the main metabolic interactions between D. vulgaris and M. maripaludis in methanogenic co-cultures. Hydrogen, carbon dioxide, and acetate are known substrates for M. maripaludis; the transfer of formate from D. vulgaris to M. maripaludis is proposed based on our data. sub-models were then combined to generate a model of syntrophic growth and simulation results were compared to experimental data. The two-species model accurately predicted several features of co-culture growth, including the ratio of the two types in culture. Finally, we tested the performance of a genetically perturbed syntrophic system, both computationally and experimentally, to gain insight into the prevalence of formate transfer in this association. These experiments showed that formate transfer is feasible, but hydrogen must be the dominant electron carrier during syntrophic growth of this species pair.

Results
Construction and evaluation of flux-balance model for syntrophic association A major goal of this study was to identify the core reactions of carbon and energy metabolism of this organism pair that account for the experimentally observed properties of syntrophic growth. We therefore constructed flux-balance models of central metabolism of D. vulgaris and M. maripaludis and combined the two to model syntrophic growth. For some of the reactions in central metabolism, there was more than one putative enzyme encoded in the genome and these enzymes used different cofactors. In addition, the stoichiometry of some reactions was unknown and could affect biomass predictions. Thus, we first explored the relevance of alternate metabolic pathways and stoichiometries by comparing model predictions for each monoculture to experimental data.

Analysis of D. vulgaris model
There are several reactions in the central metabolism of D. vulgaris for which more than one enzyme could be responsible. These include conversion of D-and L-lactate to pyruvate and oxidation of pyruvate to acetyl-CoA. Conversion of lactate to pyruvate could be catalyzed by a NADindependent (dye-linked) lactate dehydrogenase or a NADdependent dehydrogenase. Activity of a NAD-independent but not NADH-dependent dehydrogenase has been supported experimentally for several strains of Desulfovibrio, including D. vulgaris, D. vulgaris Marburg, and D. desulfuricans (Ogata et al, 1981;Stams and Hansen, 1982;Pankhania et al, 1988). These activities were measured with artificial electron acceptors (dyes). The natural intracellular electron acceptor for these dehydrogenases is still unknown, but is most likely a ferredoxin; therefore, this enzyme is referred to as ferredoxinlinked lactate dehydrogenase. In addition to the evidence for a NAD-independent, ferredoxin-linked lactate dehydrogenase, there are also several open reading frames in the D. vulgaris genome, which may encode NAD-dependent lactate dehydrogenases. These include DVU0600, which was annotated as a NAD þ -dependent L-lactate dehydrogenase by TIGR (COG39, EC 1.1.1.27), and DVU1412, which was annotated as a NAD þdependent D-isomer-specific 2-hydroxyacid dehydrogenase family protein by TIGR (COG1052, EC 1.1.1.29). To assess the potential importance of these putative reactions to the function of our model, we compared simulations of models with and without a NAD þ -dependent lactate dehydrogenase reaction in the presence of ferredoxin-linked dehydrogenase ( Table I). The presence or absence of the NAD þ -dependent dehydrogenase did not significantly affect the predicted biomass yield in simulations of D. vulgaris in syntrophic association or in monoculture with sulfate. In addition, by setting one of these fluxes to 0, we simulated the metabolic behavior of 'a deletion mutant' strain. The model predicted that the ferredoxin-linked enzyme is essential, whereas the NAD þ -linked enzyme is not. The ferredoxin-linked enzyme is necessary in the model because it is one of two sources of reduced ferredoxin, which is essential for hydrogen production and sulfate reduction. NADH supplies, on the other hand, can be replenished by several alternative D. vulgaris reactions.
The second reaction, which could be carried out by two different enzymes, is conversion of pyruvate to acetyl-CoA. The genome encodes a pyruvate oxidoreductase and a pyruvate-formate lyase. These two enzymes convert pyruvate to acetyl-CoA by different mechanisms. The pyruvate oxidoreductase produces reduced ferredoxin, whereas the pyruvateformate lyase produces formate that could be excreted by D. vulgaris for utilization by M. maripaludis when growing in co-culture. Alternatively, when growing in monoculture, it could be oxidized to CO 2 by reducing NAD þ or donating electrons directly to the electron transport chain. Simulations with these two enzymes showed that either enzyme (but not both) could be removed from the model with no effect on the predicted maximum biomass yields (Table I).
Although only cytochrome c-dependent formate dehydrogenase activities in D. vulgaris have been reported in the literature (Sebban et al, 1995), the genome sequence indicates the presence of at least three two-subunit formate dehydrogenases (DVU0587 and DVU0588; DVU2481 and DVU2482; DVU2811 and DVU2812). All three are annotated as anaerobic selenocysteine-containing formate dehydrogenases (EC 1.2.1.2, COG243). We cannot rule out that one of the D. vulgaris formate dehydrogenases could potentially reduce NAD þ . As their individual physiological functions remain unresolved, both mechanisms (as two separate reactions) were initially included in our model. However, model simulations revealed that deleting one or the other reaction (Table I) did not significantly affect biomass predictions. Thus, only the cytochrome c-dependent formate dehydrogenase was left in the model.
It is most likely that D. vulgaris uses only the interconversion (non-oxidative) branch of the pentose-phosphate pathway. It does not appear to have the genes for glucose-6phosphate dehydrogenase and 6-P-gluconate dehydrogenase. These two enzymes are key to utilization of this pathway for growth in some organisms. Biomass predictions were not affected by removing only one of these enzymes from the model. Removing both of these enzymes from the model only decreased the predicted biomass by approximately 0.6% (Table I), indicating that these steps of the pentose-phosphate pathway were not essential, but their presence may slightly improve growth of D. vulgaris on lactate.
Another interesting feature of the metabolism of this organism, predicted from the genome sequence, is the incomplete (branched) TCA cycle. The absence of genes for 2-oxoglutarate oxidoreductase (the genes for this enzyme were not found in the genome) would prevent flux through 2-oxoglutarate, breaking the cycle. Thus, D. vulgaris central metabolism might resemble that of anaerobically growing Escherichia coli and other fermenting bacteria and also some autotrophs (e.g., Beggiatoa; Hagen and Nelson, 1996) and methylotrophs (e.g., Methylobacterium extorquens; Van Dien and Lidstrom, 2002). Our model analysis showed that the system is insensitive to the addition of this enzyme.
In addition to the above issues regarding the composition of enzymes in D. vulgaris, there is also some uncertainty about the stoichiometry of reactions involved in electron transportdriven phosphorylation. More specifically, it is unclear which proton pumps function in D. vulgaris and how many protons are translocated across the membrane per electron transferred via the electron transport chain to sulfate or other electron acceptors. The ratio of H þ translocated per ATP synthesized is also unknown for D. vulgaris. These values have been determined experimentally for several species (Maloney, 1977(Maloney, , 1983Otto et al, 1984;Cook and Russell, 1994;Sturr et al, 1994;Olsson et al, 2003) and vary from 1 to 13. In this model, we assumed a ratio of 2H þ per ATP synthesized (Reaction 9). This stoichiometry is supported by good agreement between the model and published growth yields for D. vulgaris (Table II). When the model is constrained by lactate and sulfate uptake, the amount of hydrogen produced is overpredicted by several orders of magnitude; however, this agreement could not be improved through changes in the model energetic parameters. It is likely that in vivo hydrogen production is reduced because of regulatory or other effects on hydrogenases, or thermodynamic constraints that are not incorporated in the model. If the H 2 evolution rate is fixed at the experimental value, the model gives a feasible solution with acetate, sulfate, and biomass flux rates that are similar to experimentally determined levels.

Analysis of the M. maripaludis model
The energetics of M. maripaludis metabolism are not very well characterized, because it is not known exactly how many protons are pumped by each of the several potential ion pumps expressed by this methanogen. However, a redox balance can be used to determine if the model predictions are energetically feasible (Blanch and Clark, 1996). To perform this type of analysis, we rely on the fact that all the redox equivalents contained in the growth substrate must be transferred either to biomass, products, or an external electron acceptor. In other words, available reducing equivalents in all products must be equal to that of all reactants. M. maripaludis translocates protons or proton equivalents (in the form of Na þ ) in two steps of methanogenesis: the methyl transfer to coenzyme M (Reaction 11 in Supplementary information, Appendix B) and the last step to produce methane itself (Reaction 13 in Supplementary information, Appendix B). In addition, the ferredoxin-linked hydrogenase imports some amount of protons to drive the energetically unfavorable first step of methanogenesis and to convert acetyl-CoA to pyruvate (Deppenmeier, 2004;Porat et al, 2006), thus decreasing the   Biomass production rates are in mg DCW/l h.
Metabolic model of microbial mutualism S Stolyar et al net amount of excreted protons by two. We initially assumed that six protons are translocated during methanogenesis, resulting in a net of four protons per methane produced. Under these conditions, electron balance was not satisfied by the fluxes predicted by the model. Specifically, 16% of the electron potential contained in the M. maripaludis substrates (H 2 , CO 2 , and acetate) was not present in the products (biomass and methane). This outcome indicated that not enough reduced product (i.e., methane) was being produced in the predicted flux distribution. If fewer protons are pumped per carbon atom directed to oxidation pathways, more CH 4 (greater CO 2 -CH 4 flux) must be formed to support growth. By setting the number of protons translocated in the two final steps of methanogenesis equal to four (we assumed two in each step, two net protons per methane produced), the redox balance was satisfied to 2.8%. Hence, in all further simulations, Reactions 11 and 13 were assumed to pump two protons each. This value is equal to that reported for other methanogens (Deppenmeier, 2004).
To assess whether the M. maripaludis model is capable of accurately predicting the relationship between growth rate and uptake of substrates, we measured hydrogen, acetate, and carbon dioxide uptake, methane and protein production over time for three independent cultures of M. maripaludis growing on hydrogen, acetate, and carbon dioxide (see details in Supplementary information). We then performed a series of simulations and compared them to our experimental results with the aim of evaluating the model and the extent to which biomass carbon might originate from acetate versus carbon dioxide. M. maripaludis is capable of obtaining carbon autotrophically from carbon dioxide, but may obtain 60% of its carbon from acetate as well (Shieh and Whitman, 1987). During syntrophic growth, it is possible that M. maripaludis may also use the acetate produced by D. vulgaris as a carbon source. To explore how this variable affected biomass predictions of the model, we ran several simulations with the percent biomass carbon derived from carbon dioxide set at levels varying from 25 to 0%. We also ran a simulation with this feature unconstrained so that Matlab would select the configuration resulting in optimum biomass flux. As inputs for these simulations, we used the overall flux of hydrogen and acetate that we measured over 25 h of growth of M. maripaludis batch cultures on hydrogen, carbon dioxide, and acetate (see Supplementary information). Predicted CO 2 and CH 4 flux rates were similar in all simulations regardless of the percentage of biomass carbon derived from acetate versus carbon dioxide (Table III). Predicted biomass fluxes varied from 0.4 mg to 0.85 mg DCW (dry cell weight)/h and they tended to be lower when the model was constrained such that most or all biomass carbon originated from acetate and growth was therefore limited by carbon availability. All of these predicted biomass fluxes were higher than what was observed experimentally, suggesting that M. maripaludis may not fix carbon autotrophically from CO 2 in our culture conditions. Alternatively, other aspects of M. maripaludis biology may not be correctly represented in the model. For example, we assumed that 3H þ were translocated per ATP produced, but if four or more protons are required, there would be less energy available for growth from the same substrate. Making this change to the model decreased the biomass flux significantly, but not enough to match the measured flux rates. Another possibility is that the utilization of ATP is not correctly modeled. If M. maripaludis uses a significant amount of ATP to maintain osmolarity in this relatively low-salt environment or for other processes (see below), then the model may overpredict the amount of biomass that could feasibly be produced.

Analysis of the model of syntrophic association
The D. vulgaris and M. maripaludis flux-balance models were combined to form one model describing growth and metabolite accumulation when the organisms were growing syntrophically in co-culture. To model the interaction between the two species, we conceived a system of three 'compartments'. The first two 'compartments' each contained the metabolite fluxes for one of the single-species models analyzed above. These species 'compartments' could each represent the action of one cell, or the combined flux of many cells of the same species. To model the interaction between the two species, a third 'compartment' was added to the model, through which metabolites could be transferred between organisms. Exchange fluxes were added to the model in this compartment. With this modification, species could take up metabolites excreted by the other species. For syntrophic growth, the model parameters were set to reflect that there was no sulfate or hydrogen available in the culture medium. Metabolic fluxes are frequently scaled to biomass in analyses of flux-balance models. In a two-species model, however, there are two types of biomass present in not necessarily equal amounts. To avoid unnecessary confusion caused by scaling the fluxes in each species differently, we calculated our fluxes in mmol/l h as others have calculated previously (Pramanik and Keasling, 1998). Thus, in these analyses, we are using the model to predict community-level fluxes rather than the fluxes of individual cells. Initially, we explored the theoretical metabolic capability of the system by running a simulation with the minimal amount of constraints applied (Figure 2). The lactate uptake rate was set equal to 10 mmol/h, so all predicted fluxes are implicitly defined per 10 mmol lactate. Ammonia was assumed to be present in excess. Hydrogen and formate exchange fluxes were set to zero, so that all of the formate and hydrogen secreted by D. vulgaris had to be utilized by M. maripaludis. Acetate, carbon dioxide, and methane were allowed to accumulate. Finally, other potential fermentation products of D. vulgaris (glycerol and ethanol) were set to zero because no glycerol and only micromolar quantities of ethanol were detected in the medium of either the syntrophic association or monocultures of D. vulgaris (see Supplementary information for more details). D. vulgaris utilizes the primary carbon source supplied (lactate) and it is therefore reasonable to consider D. vulgaris biomass as the primary variable that is optimized ecologically. In addition, D. vulgaris dominates numerically in co-cultures by a ratio of at least 2:1 (Supplementary Figure 1B). Thus, the objective function was chosen as a combination of biomass from both organisms, with that of D. vulgaris heavily weighted. This function maximizes D. vulgaris growth as the primary objective. Secondly, of all equivalent flux distributions satisfying this objective, the one producing the most M. maripaludis was selected. M. maripaludis is considered the secondary objective because the amount of materials supplied by D. vulgaris (formate, H 2 , and acetate) constrains the optimal growth rate of M. maripaludis. Because growth of the two species is coupled by a chemical reaction, the ratio of biomasses of the two species is constrained by the stoichiometry of the reaction, and therefore should not vary widely. We wanted to explore whether this feature of syntrophic growth was captured by the model, so we compared results of three simulations that differed only in how the biomass from each species was weighted in the objective function (D. vulgaris to M. maripaludis biomass ratios 10:1, 1:1, and 1:10). In all three simulations, the predicted D. vulgaris biomass was the same, but M. maripaludis biomass was slightly higher when the objective function was heavily weighted toward M. maripaludis (Table IV). Thus, the biased function maximizing biomass of D. vulgaris was chosen for most other simulations.
The results of these initial model simulations predicted that D. vulgaris growing optimally converts the majority of the carbon contained in the substrate lactate into acetate, and that only 4.8% of the carbon is directed to biomass (Figure 2). The remaining carbon is lost as CO 2 or formate. To maintain redox balance in the absence of an external electron acceptor, evolution of 19 mmol of a reduced compound, either formate or hydrogen, is predicted. The model shows that when acetate is available, it is the preferred source of biomass carbon for M. maripaludis over the CO dehydrogenase pathway. This prediction is consistent with published data and our experimental results showed that M. maripaludis consumes acetate and presumably uses it as a carbon source (Supplementary Figure 2; Shieh and Whitman, 1987).

Relationship between the syntrophic growth model and experimental data
To assess the validity and applicability of the co-culture model, we tested how well it could predict metabolic activities across stages of syntrophic growth using only experimentally determined lactate and hydrogen flux rates. At each of several time intervals, a model simulation was run in which the experimentally determined flux rates for lactate and hydrogen were used as constraints. The predicted and experimentally determined fluxes for acetate, methane, carbon dioxide, and biomass are shown in Figure 3. The lactate and hydrogen flux measurements used as constraints in these simulations each had variability associated with them; therefore, the average flux at each interval may not represent the true flux of lactate or hydrogen. To make comparisons between the data and simulation results easier to interpret in light of possible measurement errors, we performed three simulations of the model at each time interval. In the first simulation, we used the average lactate and hydrogen flux rates. In the second and third simulations, we constrained the model with the upper or lower bounds of the 95% confidence interval for lactate and hydrogen flux rates, as determined from the variation in the measurements. These simulation results were used to form the error bars on model results expressed in Figure 3.
The flux-balance model most accurately predicted metabolite fluxes when both species were in active growth phase. This situation is most likely to occur after 50 h of growth, when the hydrogen concentration has peaked, and before 123 h of growth, when all of the lactate has been consumed. Several flux calculations during this time interval are presented in Figure 3. In most of these cases, predicted acetate and methane fluxes were very similar to measured fluxes ( Figure 3A and B). As predicted by the model, carbon dioxide and biomass fluxes increased over time as the lactate flux increased; however, the measured fluxes were lower at every time interval than predicted by the model (Figure 3C and D). Between 62 and 76 h, which corresponds to the time just after hydrogen concentration has peaked, and the last time interval when the ratio of D. vulgaris to M. maripaludis RNA is at its highest before declining (see Supplementary Figure 1C), the flux of carbon dioxide is zero. This result indicates that for a short time interval, the rate of carbon dioxide production is equal to the rate of its consumption.
To further explore the predicted behavior of the co-culture as opposed to the individual organisms growing in monoculture, we looked at key ratios in the internal flux distributions. Flux ratios are a better means for comparison across multiple conditions than are the individual fluxes, because they remove bias due solely to differences in overall carbon throughput. The ratios chosen for study are those giving the best  information about the metabolic state. The ratio of formate secretion to hydrogen secretion by D. vulgaris, as discussed in more detail in the next section, gives the relative contribution of formate to redox transfer between organisms. The ratio of pyruvate-formate lyase to pyruvate dehydrogenase indicates the relative importance of the two routes for pyruvate assimilation, whereas the ratio of pyruvate carboxylase to malic enzyme shows the contribution of the two anapleurotic routes. These flux ratios were calculated for all the time intervals of the co-culture experiment, and compared to those for the D. vulgaris monoculture experiment (Table V). The primary difference is that in the monoculture, less acetate is secreted per mole of lactate uptake, suggesting that the presence of sulfate allows a more efficient metabolism, so more of the input carbon can be converted to biomass components. In addition, the contribution of pyruvate-formate lyase to pyruvate metabolism is consistently lower in the monoculture than in co-culture. M. maripaludis has fewer branched pathways, so the metabolism is well characterized just by the external fluxes discussed previously: acetate compared to CO 2 uptake, formate versus hydrogen uptake, and CH 4 production. The model was used to explore how much non-productive ATP hydrolysis (representing futile cycles and maintenance energy) is required to make the model correctly predict the amount of biomass produced. Optimization for biomass represents an idealized metabolism with no maintenance energy, futile cycles, or ATP leaks by non-proton pumping ATPases and may therefore predict unrealistic biomass accumulation rates. Net hydrolysis of ATP by such mechanisms is common to all organisms, and serves such cellular  Figure 3 Metabolite flux rates for syntrophic growth in batch culture predicted by simulations in comparison to experimental data. Predicted (black) and experimentally determined (gray) flux rates for (A) acetate, (B) methane, (C) carbon dioxide, (D) total biomass (solid line, squares), D. vulgaris biomass (long dashes, stars), and M. maripaludis biomass (short dashes, triangles). Error bars for experimental values indicate the bounds of 95% confidence intervals calculated from four independent measurements of each sample. For simulation data, error bars were calculated by running simulations on the upper and lower bounds of the 95% confidence intervals for measured lactate and hydrogen flux rates. functions as maintaining pH gradients, intracellular pH homeostasis, proper osmolarity, and DNA repair. Futile cycling also exists when two or more enzymes catalyze a set of reactions with no net product other than the conversion of ATP to AMP, for example, pyruvate kinase and PEP synthase (Russell and Cook, 1995). Starting with the average measured lactate and hydrogen fluxes as before, the ATP hydrolysis fluxes for both organisms were manipulated so that the results matched the total biomass ( Figure 3D) as well as the ratio between D. vulgaris and M. maripaludis biomass (Supplementary Figure 1B) at each time point. The resulting values for ATP hydrolysis fluxes are listed in Table VI. These fluxes were then divided by the biomass concentration of each cell type, estimated from the experimental data, to obtain the standard units for non-growth-associated maintenance energy of mmol/ (h mg DCW) or equivalently mmol/(h g DCW). For the most part, these values are on the order of 1 mmol/(h g DCW), which agrees well with values for other organisms: 7.6 for E. coli (Edwards and Palsson, 2000b), 1.0 for Saccharomyces cerevisiae , 4.8 for Methylobacterium extorquens (Van Dien and Lidstrom, unpublished result), and 0.45 for Geobacter sulfurreducens (Mahadevan et al, 2006). This exercise was not intended to be rigorously quantitative, owing to the approximations made in biomass measurement, but rather to demonstrate that addition of a small amount of ATP drain to the model could improve accuracy of the predictions. Furthermore, it shows that small changes in maintenance requirements throughout the lifetime of the culture-for example, owing to changing environmental conditions or growth phase of the organism-can explain the changes in abundance ratio of the two organisms observed during the co-culture experiment.

Analysis of interspecies electron transfer
To gain a deeper understanding of the model and the interaction between D. vulgaris and M. maripaludis, we explored the theoretical relative importance of formate and hydrogen as interspecies electron carriers. Simulations were compared between models with either of the two different pyruvate-metabolizing enzymes that were tested in the D. vulgaris monoculture model: pyruvate-formate lyase and pyruvate oxidoreductase. In the model, if pyruvate is converted with pyruvate-formate lyase, then formate is produced and can be used as an interspecies electron carrier, but formate is not produced when pyruvate oxidoreductase is the only enzyme oxidizing pyruvate. When the model was constrained such that pyruvate was metabolized exclusively by pyruvate-formate lyase, 83% of the formate was excreted and subsequently utilized by M. maripaludis, the remainder being oxidized to CO 2 . The ratio between formate and hydrogen exchange fluxes could be varied from zero (no formate) to almost unity with no effect on predicted growth yields; however, under no case examined could any growth be supported with more formate than hydrogen. In the model simulations, when the lactate consumption rate was set to 10 mmol/h, the minimum predicted hydrogen transfer rate was 9.10 mmol, which corresponded to a formate transfer rate of 8.96 mmol/h. In the pyruvate-formate lyase mutant, formate transfer was zero and hydrogen transfer was at a maximum, 18.03 mmol/h. These results suggest that hydrogen is the dominant electron carrier, but that some proportion of the total electron equivalents produced by D. vulgaris and available to M. maripaludis could be in the form of formate. If this relationship is correct, then the total number of electrons available for M. maripaludis consumption would be lower for a mutant that was unable to consume formate than for a M. maripaludis strain capable of consuming both molecules. The total biomass and growth rate of such a mutant and its growth in co-culture may be lower. We explored this possibility by establishing a syntrophic association between D. vulgaris and the M. maripaludis mutant, MM 709 (Wood et al, 2003), which is deleted in the two annotated formate dehydrogenases. The mutant is incapable of growth on formate and shows reduced growth on H 2 :CO 2 relative to the wild type (Wood et al, 2003). A co-culture established between this mutant and D. vulgaris was capable of stable growth, demonstrating that formate transfer is not required for syntrophic growth. However, the growth rate, biomass yield, and lactate degradation rate were lower in this pairing than observed for the wild type ( Figure 4). Formate was not detected in the growth medium.

Discussion
Syntrophic interactions among anaerobic microorganisms are important components of food webs that are responsible for degrading organic materials in a variety of environments; yet the physiology of such interactions is not well understood. To develop a framework for studying the complex metabolic network of microbial food webs, we extended modeling to the level of multiple metabolic networks of interacting microbial species. This was accomplished by separating metabolites into compartments, between which there is not necessarily free and unlimited transport. Not only must there be no net accumulation or depletion of metabolites in the system, as required in general for constraint-based simulations, but the input and output fluxes of each metabolite within a particular compartment must also balance. The existence of multiple intracellular compartments has already been treated in eukaryotic cell models Forster et al, 2003). For example, in yeast there are many reactions restricted to the mitochondria whose products must be Rate of ATP hydrolysis normalized by the amount of biomass in each organism present at the sample time. Biomass concentration was estimated by culture OD (Supplementary Figure 1A) and the cell ratio of the two organisms (Supplementary Figure 1B). explicitly transported to the cytosol for further catabolism. However, this type of compartmentalization is fundamentally different from that addressed in our work, in that the compartments of existing eukaryotic models are directly connected. In a multi-organism model, the compartments represented by different microbial species are separated spatially by the extracellular medium. Consequently, the presence or absence of a single transporter in one species may greatly affect the behavior of other species in the system. Some aspects of syntrophic growth are difficult to study using classical biochemical and genetic techniques, especially the mechanism of electron transfer between species, and we anticipate that multispecies models will help evaluate alternative mechanisms and determine which are feasible. Our studies demonstrated that the single organism models successfully predicted biomass formation in relationship to metabolite production and substrate consumption for each species growing in isolation. The fusion of these models also provided quantitative predictions that helped in understanding physiology during syntrophic growth. Even though the activity of each species varied across stages of lactate degradation and methane production (see Supplementary information), the predicted fluxes of these and other metabolites were similar to those we measured at most time intervals. The model also predicted the relative abundance of the two species. In our culture conditions, D. vulgaris cells were consistently more abundant than M. maripaludis cells (see Supplementary information) by a factor of 2 or more. The fluxbalance model predicted higher D. vulgaris biomass than M. maripaludis biomass, even when it was set to primarily optimize M. maripaludis biomass. The two-species model was then used to explore the significance of two alternative mechanisms of electron transfer between species during syntrophic growth-formate and hydrogen.
Previous studies have implicated both formate and hydrogen as agents of electron transfer by comparing hydrogenase and formate dehydrogenase activities in co-culture versus axenic cultures, or by testing the ability of hydrogen versus formate consumers to cooperate with a particular syntroph (Schink and Stams, 2002). However, the prevalence of either mechanism in syntrophic interactions is unknown (Schink, 1997). Thus, a key objective of the integrated model was to explore the possible metabolic constraints on the use of either hydrogen or formate as an agent of electron transfer. We first confirmed that both organisms had the necessary gene content for production and utilization of both hydrogen and formate. We then tested the relevance of these genes by assessing their function in the context of the network of metabolic reactions that composed the central metabolism of the two organisms. The resulting model predicted that both molecules could be used as electron carriers simultaneously, but that only hydrogen was essential. Our experimental results confirm this model prediction because syntrophic growth was possible even when M. maripaludis was genetically incapable of utilizing formate. Our flux-balance model results also suggested that hydrogen transfer should be dominant over formate transfer during degradation of lactate because D. vulgaris cannot produce enough formate for the methanogen and at the same time satisfy the internal redox balance. A dominant role of hydrogen transfer with possible augmentation of electron flux through formate transfer might also be inferred from the mild deficiency observed for syntrophic growth of the co-culture containing a formate dehydrogenasedeficient M. maripaludis strain. However, formate was not detected in the culture medium and both co-cultures produced the same amount of methane, suggesting that the availability of electrons for M. maripaludis was roughly similar. Thus, it is unclear whether this growth deficiency is a consequence of other pleiotropic effects of formate dehydrogenase genes on M. maripaludis or D. vulgaris growth, or both. One potential problem with using metabolic models to understand complex communities is that the network of metabolic reactions within and among species is very complex. Our development of a predictive two-species model is especially interesting, given that the model itself was so simple. Most flux-balance models published to date are comprised of hundreds of reactions, but our combined twospecies model (core model) contained mainly reactions involved in central metabolism and therefore had only 170 reactions and 147 metabolites. Core models of this type are the first step toward a larger modeling effort, involving construction of genome-scale models. Modeling at the genome-scale has the advantage that the full solution space can be explored, leading to potential phenotypes that cannot be predicted with the core model. Not only does genome-scale modeling involve the inclusion of full biosynthetic pathways, but it also expands the range of substrates and potential secreted products. The size of the solution space would also be expanded due to catabolic and interconversion pathways (e.g., amino acids).
Although expansion of the model may improve the fit, it also presents additional difficulties, such as the tendency to overestimate the impact of pathways that may only be expressed at low levels. These problems could be addressed with an even further elaboration on the constraint-based modeling framework, which is the inclusion of regulatory rules to moderate certain flux constraints . However, application of these advanced modeling techniques requires a large amount of effort, both in characterization of the organisms and in constructing the model itself. The core models used here represent an initial step toward using models to understand microbial interactions. They also demonstrate that genome-scale models are not necessarily essential for understanding basic physiological phenomena. This general approach was already proven useful for single organism models Van Dien and Lidstrom, 2002), and now our modeling results demonstrate applicability of simplified metabolic models for predicting structures of microbial communities. With this validation complete, future work will focus on development and use of genome-scale models.
To our knowledge, this is the first flux-balance model of Desulfovibrio and an association involving two species. This advance in modeling has implications for research on the ecology of microbial communities in natural and engineered systems. The basic composition of microbial communities is to a large extent based upon the successive transfer of metabolites from one species to another, which is capable of metabolizing it further to obtain energy or use it as a building block for biosynthesis. Metabolic models may therefore be useful for predicting how the availability of different substrates is likely to affect the abundance of species in microbial food chains such as those that occur in anaerobic digesters ( Figure 1). Thus, our development of a relatively simple model that captures the key features of syntrophic growth suggests that this technique should be useful for studying increasingly complex microbial communities.

Materials and methods
The D. vulgaris model used here was based on genome annotation (Heidelberg et al, 2004) and the Microbes online database (http:// www.microbesonline.org/). As the focus of our study is on fundamental reactions of carbon and energy metabolism in the two organisms, the complete model was simplified to include only central metabolism and amino-acid biosynthesis pathways. All other biosynthetic reactions were condensed into single equations for each major biomass component: DNA, RNA, protein, carbohydrate (assumed to be glycogen), phospholipid, peptidoglycan, and lipopolysaccharide. A reaction was then written to combine these components in the proper ratio for biomass synthesis. We assumed that biomass composition of D. vulgaris is the same as in E. coli (Pramanik and Keasling, 1997). Using this composition, 1 g DCW of biomass is equal to 43.4 mmol of carbon atoms. This value was used to calculate biomass yield on a carbon atom basis. The resulting model contained 86 reactions and 75 intracellular metabolites (Appendix A).
A model of M. maripaludis central metabolism was constructed on the basis of putative gene functions from the genome annotation (Hendrickson et al, 2004), supplemented by current physiological knowledge of the organism. Three pathways were included for nitrogen assimilation: NH 4 þ uptake from the medium, N 2 fixation, and utilization of externally supplied alanine. Carbon uptake was allowed to occur from CO 2 via the CO dehydrogenase/acetyl-CoA synthase pathway, or from acetate via acetyl-CoA synthetase. Formate could also be oxidized to CO 2 . Although there are gaps in some of the amino-acid biosynthesis pathways (Hendrickson et al, 2004), the complete pathways were included because this organism is not auxotrophic for any amino acids. As with the D. vulgaris model, all other biosynthetic pathways were condensed. Because there have been no reports of detailed biomass composition measurements in Archaea, in constructing the model it was assumed that DNA, RNA, phospholipids, and protein combine in the same proportion as in E. coli (Pramanik and Keasling, 1997). The methanogen does not have an outer membrane. Instead it possesses an S layer; therefore, peptidoglycan and lipopolysaccharide were not included. Although this biomass composition is a crude estimate, it has been demonstrated that changing macromolecular composition does not have a large effect on central carbon or energy metabolism (Edwards and Palsson, 2000b). The final M. maripaludis model contained 82 reactions and 72 intracellular metabolites (Appendix B).
For the simulation of syntrophic growth, the two metabolisms were treated as separate compartments within the same model, similar to the way by which different organelles are treated in eukaryotic cell models Forster et al, 2003). A third 'compartment' was added to represent the culture medium, through which certain metabolites could be transferred between the two organisms. Transport fluxes were added for metabolites that could be taken up by or secreted from each organism. Finally, exchange fluxes were added to represent accumulation or depletion of metabolites from the bulk medium. These exchange fluxes were either set equal to zero (indicating no net accumulation), set to the experimental measurement, or allowed to be flexible. The particular conditions for each simulation are described in the Results section. The final co-culture model contained 170 reactions and 147 metabolites. Note that the same compounds in different compartments are treated as separate metabolites and are connected by transport fluxes.

Simulations
Flux-balance analysis involves writing a mass balance on each metabolite in the network, applying a steady-state assumption (Savinell and Palsson, 1992;Vallino and Stephanopoulos, 1993), and solving the resulting system of equations. As there are more unknown reaction rates (metabolic fluxes) than metabolites, the system of equations is underdetermined and thus does not have a unique solution. A common method to obtain a unique flux distribution in such a case is to define an objective function, such as biomass production, and use linear programming to find the optimal solution with respect to this objective function (Varma and Palsson, 1994;Pramanik and Keasling, 1997). Simulations were performed using the Fluxanalyzer program (Stelling et al, 2002;Klamt et al, 2003), which operates within the Matlab environment (The Mathworks, Natick, MA). The objective considered, unless stated otherwise, was maximization of the following function: 1.0 Â D. vulgaris biomass þ 0.1 Â M. maripaludis biomass. This resulted in the maximization of D. vulgaris growth as the primary objective, with M. maripaludis as a secondary objective. There are possibly multiple solutions maximizing D. vulgaris growth, and this function assures that the solution that also maximizes M. maripaludis growth is chosen.

Construction and characterization of syntrophic interaction
Strains and culture conditions D. vulgaris Hildenborough (NCIMB 8303) was obtained from J Wall (Columbia, MO). Wild-type Methanococus maripaludis S2 (Whitman et al, 1986) and a M. maripaludis mutant strain MM 709 lacking both annotated formate dehydrogenases (Wood et al, 2003), were used in this study.
M. maripaludis S2 was routinely cultured in mineral McC medium amended with 10 mM acetate (Jones et al, 1983) at 371C in tubes or bottles with B80% headspace consisting of 80% H 2 :20% CO 2 at 30 psi and constant shaking at 200 r.p.m.
Although M. maripaludis and D. vulgaris should theoretically be able to grow syntrophically in the absence of sulfate, they are typically cultivated under different concentrations of sodium and magnesium. Before establishing a syntrophic culture, it was necessary to find concentrations for these two ions that were permissible for both organisms. To develop a medium suitable for syntrophic growth with D. vulgaris, the effect of sodium and magnesium concentration on M. maripaludis was investigated by modifying the sodium and magnesium chloride concentrations in McC as follows: 17 mM NaCl and 2 mM MgCl 2 , 50 mM NaCl and 4 mM MgCl 2 , 120 mM NaCl and 5.9 mM MgCl 2 . The methanogen demonstrated comparable growth rate at all ion concentrations tested. Growth of D. vulgaris was not affected at magnesium sulfate concentrations as high as 40 g/l, but sodium chloride slowed growth at concentrations as low as 2 g/l (the extra amount added to the media). On the basis of these results, we have modified McC media to accommodate syntrophic growth. This B3 medium (see the composition below) without sulfate added was used for most experiments with co-cultures.
D. vulgaris and M. maripaludis were cultured syntrophically at atmospheric pressure in B3m or B3n medium without Na 2 SO 4 at 301C without shaking in an 80% N 2 :20% CO 2 headspace that composed B30% of the culture vessel (25 ml Balsch tube or 250 ml bottles). As necessary, pressure in the headspace was monitored using a needle pressure gauge. Syntrophic associations were initiated by mixing equal volumes of stationary phase D. vulgaris and M. maripaludis cultures, removing excess culture medium by centrifugation and resuspending in B3n medium lacking sodium sulfate. Samples of established syntrophic associations and D. vulgaris and M. maripaludis samples were stored at À801C in 10% glycerol.
Microbial contamination was monitored by microscopic observation and by plating co-cultures on TSA (Difco) and R2A (Difco) agar and incubating plates at 301C in aerobic and anoxic conditions.

Experimental determination of growth and metabolic activity
To characterize growth of the syntrophic association between D. vulgaris and M. maripaludis, four 240 ml bottles containing 177 ml B3n media and an 80% N 2 :20% CO 2 headspace were inoculated with 10 ml of a stationary phase culture of the syntrophic association. This inoculum had been previously grown in B3m medium in batch culture by serial transfer through eight transfers at a 20-fold dilution rate. Thus, before the start of our experiment, both species had approximately 35 generations to acclimate to the syntrophic growth conditions. Immediately after inoculation and at several time intervals for the next 176 h, samples were taken from each bottle for measurement of total biomass, ratio of cell types, 16S rRNA ratio quantification, and lactate, acetate, methane, carbon dioxide, and hydrogen concentration. Pressure was measured at each time interval. The methods used for these assays are described below.
To characterize growth of M. maripaludis in monoculture, 0.5 ml of an overnight culture of M. maripaludis was inoculated into three 240 ml bottles with an 80% H 2 :20% CO 2 headspace at 30 psi pressure. Each bottle contained 60 ml of B3n media amended with 10 mM Na acetate, 2 g yeast extract per liter, and higher concentrations of NaCl (4 g/l), cysteine (3 mM), and sodium sulfide (0.5 g/l). At several time intervals for 25 h, pressure of the headspace was measured and samples were taken for measurement of biomass, lactate, acetate, methane, carbon dioxide, and hydrogen.

Analytical methods
The concentration of organic acids (lactate, pyruvate, acetate, formate, and fumarate) and inorganic ions (sulfate, phosphate) in culture fluids were determined using a Dionex 500 system equipped with an AS11HC column. In some cases, the concentration of organic acids was also measured on an HPLC equipped with an HPX 78 (Bio-Rad) column. Hydrogen concentrations were determined with a RGD2 Reduction Gas Detector (Trace Analytical) with 60/80 MOLE SIEVE 5A column (72 Â1/8 inch) with N 2 as carrier gas. The concentration of methane and carbon dioxide was measured on a GC equipped with a TCD and 80/100 HAYESEP Q column (72 Â1/8 inch) with helium as carrier gas.
The concentrations of ethanol and glycerol in the syntrophic culture supernatant were measured enzymatically. Ethanol was measured using alcohol dehydrogenase (Sigma A3263) and aldehyde dehydrogenase (Sigma A6338), and glycerol using the enzyme glycerol dehydrogenase, both using the protocols suggested by Sigma.

Biomass, protein, and RNA measurements
For protein measurements, cell pellets were resuspended in 250 mM NaOH and frozen at À201C. Protein was measured using Bradford reagent (Pierce) according to the manufacturer's recommendations. It was assumed that the composition of D. vulgaris and M. maripaludis cells was 60% protein in all culture conditions tested. Thus, the experimentally determined protein concentration was multiplied by 2 to estimate total biomass in most cases. Alternatively, the biomass of the syntrophic culture was estimated in some cases from the optical density using a conversion factor of 1 OD 600 ¼0.385 g/l dry weight. This conversion factor was determined by drying biomass of a syntrophic culture on filters overnight at 1051C.
At each time interval in syntrophic growth experiments, a biomass sample was frozen at À801C. RNA was extracted from these samples using the Master Pure kit (Epicentre Biotechnology) according to the manufacturer's directions. The RNA in these samples was subsequently separated and analyzed using capillary electrophoresis in a Bioanalyzer (Agilent) according to the manufacturer's instructions. The 16S rRNA fragments of D. vulgaris and M. maripaludis migrate at slightly different rates and produce separate peaks in this electrophor-esis method. The ratio of the peak areas of the two species was calculated to determine the 16S rRNA ratio.

Microscopic observation of syntrophic association
During cultivation, samples of the syntrophic association were viewed under the light microscope, Olympus BH2, to check for purity. During growth characterization studies, 100 ml of the culture was immediately fixed in a solution composed of 225 ml phosphate-buffered saline and 75 ml of 37% formaldehyde, and stored at 41C for at least 24 h. These samples were subsequently stained for 24 h by addition of 50 ml of a 500 mg/ml solution of dichlorotriazinylaminofluorescein, a fluorescent protein stain. Diluted samples were filtered onto a black nitrocellulose membrane with 0.22 mm pore size, mounted onto a microscope slide and viewed under an epifluorescence microscope LSM5 PASCAL (Zeiss). The ratio of D. vulgaris-shaped cells (vibriod) and M. maripaludis-shaped cells (cocci) was counted on 10 randomly selected fields for each slide.

Supplementary information
Supplementary information is available at the Molecular Systems Biology website (www.nature.com/msb).