Thermodynamics and H2 Transfer in a Methanogenic, Syntrophic Community

Microorganisms in nature do not exist in isolation but rather interact with other species in their environment. Some microbes interact via syntrophic associations, in which the metabolic by-products of one species serve as nutrients for another. These associations sustain a variety of natural communities, including those involved in methanogenesis. In anaerobic syntrophic communities, energy is transferred from one species to another, either through direct contact and exchange of electrons, or through small molecule diffusion. Thermodynamics plays an important role in governing these interactions, as the oxidation reactions carried out by the first community member are only possible because degradation products are consumed by the second community member. This work presents the development and analysis of genome-scale network reconstructions of the bacterium Syntrophobacter fumaroxidans and the methanogenic archaeon Methanospirillum hungatei. The models were used to verify proposed mechanisms of ATP production within each species. We then identified additional constraints and the cellular objective function required to match experimental observations. The thermodynamic S. fumaroxidans model could not explain why S. fumaroxidans does not produce H2 in monoculture, indicating that current methods might not adequately estimate the thermodynamics, or that other cellular processes (e.g., regulation) play a role. We also developed a thermodynamic coculture model of the association between the organisms. The coculture model correctly predicted the exchange of both H2 and formate between the two species and suggested conditions under which H2 and formate produced by S. fumaroxidans would be fully consumed by M. hungatei.


Introduction
Microorganisms in nature engage in a variety of interactions with other species in their environment. Syntrophy is one such type of inter-species interaction in which one species lives off the metabolic by-products of another [1][2][3]. Synthetic methanogenic communities [4] are typically tightly constrained by thermodynamics, as the oxidation reactions carried out by the first community member are thermodynamically unfavorable unless the degradation products are maintained at low levels by the second community member [5].
Genome-scale models can also be used to study the relationship between thermodynamics and metabolism, by ensuring that network predictions are consistent with thermodynamic principles [29][30][31][32][33][34]. In this study, we used thermodynamics-based metabolic flux analysis (TMFA) to develop a thermodynamic, coculture model of the syntrophic association between the anaerobic bacterium Syntrophobacter fumaroxidans and the methanogenic archaeon Methanospirillum hungatei. In association with M. hungatei, S. fumaroxidans converts propionate to acetate, CO 2 , and H 2 [35][36][37]. CO 2 and H 2 can be interconverted to formate [38][39][40], with H 2 and formate serving as the electron carriers between the two species. H 2 and formate production are only observed during syntrophic growth. Using a thermodynamic, constraintbased model, we set out to test the proposed hypothesis that this behavior is governed by thermodynamics [3,5,6,8,9].
We developed genome-scale metabolic reconstructions of both microorganisms, and verified proposed mechanisms of ATP production within each individual species. Additional constraints and a cellular objective function were identified to predict the proper flux through experimentally characterized carbon and electron transport pathways during monoculture and syntrophic (i.e., coculture) growth. Our analysis revealed that thermodynamic constraints alone are insufficient to explain why S. fumaroxidans does not produce H 2 in monoculture.
We also extended TMFA to model the syntrophic association between the two microorganisms. The association is modeled as a continuous coculture system with constraint-based models for each microbe and a mass balance around the reactor. Similar to cFBA [28], the coculture model accounted for the biomass concentrations of each species. We predicted the behavior of this syntrophic association under a variety of dilution rates, and identified regimes of behavior consistent with experimental observations.

Results
Metabolic reconstructions of S. fumaroxidans (iSfu648) and M. hungatei (iMhu428) were tested and parameterized using experimental data from growth on single substrates in monoculture and coculture. The iSfu648 and iMhu428 models were used to verify proposed mechanisms of energy conservation within each species. A coculture model was then developed to identify conditions where H 2 and formate, produced by S. fumaroxidans, is fully metabolized by M. hungatei.

Testing and Parameterizing the iMhu428 Metabolic Model
The iMhu428 reconstruction of M. hungatei was built from the iMB745 reconstruction of M. acetivorans [41]. A preliminary draft reconstruction was built from iMB745 using the RAVEN Toolbox [42] and the KEGG SSDB [43]; however, M. hungatei orthologs were found for only 428 of the 745 genes in iMB745. To avoid extensive gapfilling, reactions from the iMB745 were copied into the M. hungatei reconstruction, with modifications to reflect key metabolic features of M. hungatei (see S1 Text). As a consequence, the iMhu428 reconstruction is a draft reconstruction requiring further evaluation. A thermodynamic model for iMhu428 was built and TMFA was used to predict ATP generating mechanisms in minimal media monoculture conditions (see S1 Table in S1 Dataset for constraints used).
Experimental evidence suggests that M. hungatei is able to generate 0.5 mole ATP per mole of CO 2 converted to CH 4 [7], via the metabolic route shown in Fig 1. In order for this route to be thermodynamically feasible, the Δ r G '0 of one reaction (FMFTSPFT, formylmethanofurantetrahydromethanopterin formyltransferase) had to be allowed to vary within a 99% confidence interval of its estimated standard transformed Gibbs free energy of reaction (D r G 0 0 est ) (rather than the 95% interval used for all other reactions) in order to carry flux in the proper direction. The reactions for carbon source utilization in M. hungatei are well-characterized [7,[44][45][46][47], but uncertainty remains about the stoichiometry of small ion transport [7]. Na + transport stoichiometries associated with tetrahydromethanopterin S-methyltransferase (MTSPCMMT_CM5HBCMT, E.C. 2.1.1.86) and a Na + /H + antiporter (NAT3_1) were selected to give an ATP yield matching the experimental estimates: two Na + ions exported by MTSPCMMT_CM5HBCMT, and a one Na + per H + transported by NAT3_1. However, different stoichiometries for these reactions are also thermodynamically possible (see Discussion).
Experimental measurements of growth rates, yields, and maintenance costs were also used to identify substrate uptake rates (SUR) for CO 2 and formate, and the growth-(GAM) and non-growth-associated (NGAM) ATP maintenance requirements for M. hungatei, as described in S1 Text. NGAM represents the amount of energy spent to maintain the cell (i.e., maintenance energy), while GAM represents energy spent on growth-related functions (e.g., protein synthesis). For the iMhu428 model, the NGAM was estimated to be 0.6 mmol ATP/ gDW/day, GAM was estimated to be 47 mmol ATP/gDW, SUR CO2 was estimated to be 75.7 mmol/gDW/day, and SUR formate was estimated to be 955 mmol/gDW/day.

Testing and Parameterizing the iSfu648 Metabolic Model
The iSfu648 reconstruction of S. fumaroxidans was built from the KEGG database using the RAVEN Toolbox [42]. The resulting draft reconstruction was manually refined (see S1 Text), with particular attention paid to ATP production mechanisms. A number of studies have identified gene clusters encoding a variety of hydrogenases, dehydrogenases, and other electron transport enzymes [8,9,[48][49][50][51], whose expression levels vary across growth conditions [51]. All told, 17 enzymes which catalyze 12 different electron transport reactions have been identified (S3 Table in S1 Dataset). In many cases, the proposed reactions catalyzed by these enzymes differ between studies; a brief description of each reaction and justification for each annotation is given in S1 Text. The draft reconstruction was updated to be consistent with the reported carbon utilization and electron transport reactions, and the resulting stoichiometric model was converted to a thermodynamic model.

ATP Production Mechanisms in S. fumaroxidans
Experimental studies have elucidated five growth modes for S. fumaroxidans: four in monoculture and one in coculture with M. hungatei (S1 Text) [36,48,52]. This work examines the three most commonly studied growth modes (Table 1): monoculture growth on fumarate, monoculture growth on fumarate plus propionate, and coculture growth on propionate.
A variety of experimental findings were synthesized to develop theoretical flux distributions for these three growth modes (Fig 2 and S1 Text). These experimental findings suggested additional regulatory and flux-coupling constraints for the iSfu648 model, such as coupling between fumarate reductase and the cytosolic hydrogenase due to co-localization in the membrane (see S2 Table in S1 Dataset and S1 Text for the full set of constraints and their justification).
During monoculture growth on fumarate alone (Fig 2A), one mole of fumarate gets fully oxidized to CO 2 , while six moles of fumarate get reduced to succinate [48,52]: The oxidation of one fumarate to CO 2 generates one ATP and five reducing equivalents (three NADH and two pairs of reduced ferredoxin) [48,52], while the reduction of additional fumarate to succinate by fumarate reductase (FRD) consumes reducing equivalents (menaquinol) [48]. Electrons are transferred from NADH and reduced ferredoxin to menaquinone through the combined action of the Rnf complex (RNF), the ferredoxin-oxidizing hydrogenase (frH 2 ase), the cytosolic hydrogenase (cytH 2 ase) and formate hydrogen lyase (FHL). The reduction of fumarate to succinate also generates the proton motive force (PMF) responsible for driving the RNF reaction and producing ATP.
During monoculture growth on fumarate plus propionate (Fig 2B), one mole of propionate gets oxidized to succinate, while one mole of fumarate gets oxidized to acetate and CO 2 . Two additional moles of fumarate get reduced to succinate [36,48]: The oxidation of fumarate to acetate and CO 2 produces one NADH and one pair of reduced ferredoxin, while the reduction of fumarate to succinate by FRD consumes menaquinol. Electrons are transferred from NADH and reduced ferredoxin to menaquinone through the combined action of the confurcating hydrogenase (cH 2 ase) and cytH 2 ase. Oxidation of propionate to succinate produces one ATP, while FRD generates the PMF necessary for additional ATP production.
During coculture growth on propionate (Fig 2C), propionate gets oxidized to acetate and CO 2 via the methylmalonyl-CoA pathway [48,52]: ATP is generated during the oxidation of propionate to succinate, and this ATP establishes the PMF necessary to drive the endergonic oxidation of succinate to fumarate (SDH), producing menaquinone. cytH 2 ase then transfers electrons from menaquinol to two protons, generating H 2 . The oxidation of fumarate to acetate and CO 2 produces one NADH and one pair of reduced ferredoxin, and cH 2 ase couples NADH and ferredoxin re-oxidation with H 2 production. Unlike in the monoculture growth modes, the H 2 is not consumed intracellularly and must diffuse outside the cell. It has been proposed that the net production of H 2 by S. fumaroxidans is only thermodynamically favorable at the low H 2 concentrations maintained by methanogens, thereby explaining why S. fumaroxidans only produces H 2 during coculture growth. S. fumaroxidans exhibits considerable flexibility in its ATP production mechanisms during coculture growth ( [8,9,48,51]), and can produce formate instead of CO 2 (Fig 3) yielding an overall transformation of: In one mechanism (Fig 3A), activity of the cytosolic formate hydrogenase (cytFDH) substitutes for the activity of cytH 2 ase. In a second mechanism (Fig 3B), the confurcating formate dehydrogenase (cFDH) substitutes for cH 2 ase. Here, cFDH couples NADH and ferredoxin reoxidation with the conversion of CO 2 (from propionate oxidation) to formate.

Model Predictions of ATP Production by S. fumaroxidans
Experimental evidence and conceptual models of S. fumaroxidans energy metabolism suggest that the carbon and electron transfer pathways shown in Fig 2 provide the sole source for ATP production in S. fumaroxidans, either by substrate-level phosphorylation or through establishment of a proton gradient used by ATP synthase [5]. To test the computational model's predictions, TMFA was used to maximize ATP production under each of the three growth modes.
When flux was restricted to a reduced network containing all the reactions shown in Fig 2 (listed in S4 Table in S1 Dataset), the iSfu648 model correctly predicted the flux distributions shown in Fig 2. However, when flux was allowed throughout the entire network, additional flux distributions with higher ATP yields were identified. Additional reaction direction constraints were developed to ensure model-predicted flux distributions matched experimental observations (see S2 Table in S1 Dataset and S1 Text for details). However, the resulting flux distributions are not fully consistent with the hypothesis that S. fumaroxidans has adapted to maximize its energy yield.

Evaluating Thermodynamics of H 2 Production
Experimental studies of S. fumaroxidans have shown that H 2 is not produced during growth in monoculture [36,52,53], and it is widely thought that H 2 production is only thermodynamically favorable at low partial pressures [3,6,8,9]. In particular, methanogens in syntrophic communities enable sustained H 2 production by consuming H 2 and keeping its partial pressure low [3,5,6,8,9]. Indeed, when H 2 production was observed in monoculture, H 2 production ceased at a partial pressure of approximately 10 Pa [53].
However, when maximizing H 2 production under monoculture conditions, simulations reveal that H 2 production remains thermodynamically feasible. For example, during monoculture growth on fumarate, the iSfu648 model predicts that H 2 can be produced via the following In this scenario, H 2 molecules produced by the ferredoxin-oxidizing hydrogenase are exported outside the cell, instead of serving as substrates for the cytosolic hydrogenase. As a result, no PMF is generated by fumarate reductase, and the net ATP yield is zero. Thus, while H 2 production remains thermodynamically possible, H 2 production is only associated with sub-optimal mechanisms of ATP generation. This suggests that thermodynamic considerations alone may not explain the absence of H 2 production during monoculture growth, but that the observed flux distribution may instead be driven by demands for energy generation.
While H 2 production was not predicted for monoculture conditions when ATP production was maximized, H 2 production was initially predicted when growth was instead maximized. To eliminate monoculture H 2 production in the TMFA model, we first sought to constrain ratios of metabolite concentrations with an approach similar to that used to correct TMFA growth predictions [34] (see S1 Text for details). While metabolite ratio constraints could be identified to prevent some H 2 production mechanisms, H 2 production during monoculture growth could not be completely eliminated. If thermodynamics prevents H 2 production in monoculture, then the current thermodynamic model may contain too much uncertainty in its Gibbs free energy estimates (see Discussion). Regulatory effects could also potentially prevent H 2 production in monoculture conditions. To correct the model, all subsequent monoculture simulations were performed by preventing H 2 production.

Parameterization of the iSfu648 Metabolic Model
Model parameters were estimated after reaction direction constraints were added to the iSfu648 model (to be consistent with reported ATP generation and H 2 production mechanisms). Experimental measurements of growth rates, yields, and maintenance costs were used to identify the SURs, GAM, and NGAM parameters for S. fumaroxidans. These parameters were estimated using data from monoculture growth on fumarate alone and coculture growth on propionate alone (see S1 Text). For the iSfu648 model, the following parameters resulted in the best fit of the model to the experimental data: NGAM = 3.36 mmol ATP/gDW/day, GAM = 22.8 mmol ATP/gDW, SUR propionate = 37.7 mmol/gDW/day, and SUR fumarate = 27.6 mmol/gDW/day.
Using these parameter values, the in silico growth rates under each growth condition were predicted (Table 1). Not surprisingly, the predicted growth rates for fumarate alone and propionate alone conditions agree with experimental observations (since these were used to estimate the parameter values). However, the model significantly under-predicts the measured growth rate during monoculture growth on fumarate plus propionate (0.55 days -1 predicted, 0.73 days -1 observed). This discrepancy could be caused by differences in uptake rates or maintenance costs in the fumarate plus propionate condition compared to the conditions with propionate alone or fumarate alone.

Uptake and Secretion at Maximal Growth
When maximizing biomass production on the entire network, the iSfu648 model predicted a wide range of product secretion rates. When the enzyme cost (i.e., total flux) was minimized at the maximum growth (pTMFA [54], see Methods) the model-predicted product yields closely matched reported values for two of the three growth modes (Table 1)-monoculture growth on fumarate alone and coculture growth on propionate alone. These results indicate that the majority of carbon is diverted to fermentation products, consistent with the expectation that high fluxes through the low-energy fermentation pathways are needed to meet cellular energy demands.
However, for monoculture growth on fumarate plus propionate, the model failed to predict that fumarate and propionate should be consumed at the observed ratio of approximately three fumarate per propionate [36]. Instead, the model predicted both substrates would be consumed at their maximum SURs, resulting in a ratio of 0.73 fumarate per propionate. The experimentally observed 3:1 ratio is thought to arise due to coupling within the metabolic network ( Fig  2B), since oxidation of one fumarate produces one CO 2 (used to oxidize one propionate) and two pairs of electrons (used to reduce two fumarate). While this coupling arises naturally on the reduced network, the full metabolic network enables alternative coupling mechanisms (not shown) that permit other fumarate to propionate ratios. Since propionate oxidation generates carbon precursors and ATP for biomass, maximizing biomass production results in the model under-predicting the fumarate to propionate uptake ratio. Increasing the ratio of fumarate to propionate uptake (to 3:1) decreases the predicted propionate SUR and growth rate (S1 Fig in S1 Text), implying the experimental ratio is sub-optimal with respect to growth maximization. Instead of constraining SURs (since values were not reported in the literature), we incorporated a fumarate to propionate SUR ratio constraint for this condition.
While the predicted product yields closely matched experimental observations (after imposing the SUR ratio constraint), the pTMFA-predicted intracellular flux distribution during monoculture growth on fumarate and propionate substantially deviated from that shown in Fig 2B. Further constraints on reaction directions were required so that propionate and fumarate were metabolized in the model via the pathways shown in Fig 2B (results not shown). Taken together, the need for constraints on reaction directions and fumarate to propionate uptake ratio suggests that neither maximization of biomass nor minimization of enzyme cost are sufficient to explain the fluxes of S. fumaroxidans growing in this monoculture condition.
Cocultures of M. hungatei and S. fumaroxidans have been grown in both batch and continuous (chemostat) systems. During batch growth, H 2 pressure rose during the lag phase and became constant during exponential growth [53]. Continuous cultures also exhibited constant H 2 partial pressure [53]. However, to the best of our knowledge, measurements for the relative ratios of M. hungatei to S. fumaroxidans at constant H 2 pressure (where H 2 consumption and production rates are balanced) have not been reported. Instead, an overall reaction for the coculture of propionate ! acetate þ 0:25 is frequently discussed [5,53], which can occur at a ratio of three M. hungatei to four S. fumaroxidans. Since initial conditions for batch experiments were not reported, a continuous culture model was constructed and first evaluated using a 3:4 relative biomass ratio (M. hungatei: S. fumaroxidans). The continuous coculture model included constraint-based models for each microbe and mass balances around the reactor (S2 Fig in S1 Text), and accounted for the biomass concentrations of each species. Both species were constrained to grow at the dilution rate, and the model minimized the species-weighted total flux through the two metabolic networks (pTMFA). Propionate was the only substrate in the reactor feed (i.e., it had a net flux into the reactor), ensuring that all carbon and electrons used by M. hungatei were produced by S. fumaroxidans.
Predicted yields around individual species were first evaluated (Fig 4A). At low dilution rates, S. fumaroxidans was predicted to convert propionate to acetate, H 2 , and CO 2 /formate and M. hungatei was predicted to convert CO 2 /formate to CH 4 (Fig 4B and 4C, which show alternate solutions with CO 2 or formate being exchanged). As the reactor dilution rate increased, the predicted species yields of H 2 and CO 2 /formate (S. fumaroxidans) and CH 4 (M. hungatei) decreased slightly (Fig 4B and 4C). In addition, species' uptake rates of propionate and CO 2 /formate increased with dilution rate, as biochemical transformation of these substrates provides the energy needed for cellular growth and maintenance.
Species' uptake rates, secretion rates, and relative biomass ratio affects overall bioreactor yields, and these bioreactor yields were subsequently investigated (Fig 4D). At a 3:4 relative biomass ratio (M. hungatei: S. fumaroxidans), M. hungatei did not fully utilize all of the H 2 and CO 2 /formate produced by S. fumaroxidans, even at high dilution rates (Fig 4E and 4F). The net H 2 production by the community indicates that S. fumaroxidans produces more H 2 than M. hungatei's needs and suggests the community can maintain higher M. hungatei to S. fumaroxidans ratios or that S. fumaroxidans could support faster growth of M. hungatei. At a dilution rate of 0.05 days -1 , S. fumaroxidans produces H 2 in excess of M. hungatei's energy needs until the relative M. hungatei to S. fumaroxidans biomass ratio reaches approximately 1.6:1 (Fig 5). These simulations suggest that invariant external H 2 concentration requires high ratios of M. hungatei to S. fumaroxidans (higher than has been proposed in the literature based on overall reaction stoichiometries), M. hungatei growing at faster rates than S. fumaroxidans (e.g., in batch culture), or a combination of the two.
Furthermore, studies have shown that in coculture, S. fumaroxidans passes electrons to M. hungatei via formate, as well as H 2 [38,39,56]. Coculture simulations predicted that formate  could be exchanged in lieu of CO 2 (Fig 4B and 4C), without affecting the predicted bioreactor yields or species-weighted total flux (pTMFA objectives). When formate is exchanged the formate dehydrogenases of S. fumaroxidans and M. hungatei facilitate the interconversion of formate to CO 2 and H 2 . Finally, steady-state metabolite concentrations in the coculture were predicted using thermodynamic variability analysis [33,34] (S1 Text). Similar to a previous study of E. coli [34], the majority of steady-state metabolite concentrations were not constrained by thermodynamics (i.e., the concentration ranges were the global concentration bounds of 0.01mM and 20 mM). However, hypotheses for some extracellular (Table 2) and intracellular (S5 Table in S1 Dataset) metabolite concentrations could be made. The model predicts that propionate in the media must be greater than 0.36 mM for coculture growth to occur, and that acetate and CO 2 concentrations must be less than 4mM and 0.65 mM, respectively. The model also predicts critical concentrations for H 2 (0.0032 mM) and formate (0.0020 mM, when formate is being exchanged), which must be maintained in order for methanogenesis to occur. All of these predictions were insensitive to both the dilution rate and biomass ratio (M. hungatei: S. fumaroxidans).

Discussion
The iMhu428 and iSfu648 thermodynamic models were successfully used to verify proposed carbon and ATP production pathways in M. hungatei and S. fumaroxidans; however, some more efficient pathways (based on ATP or biomass yields) were found in some cases. The computational results highlighted topics that necessitate further discussion relating to (1) the stoichiometry of transport reactions, (2) the thermodynamics of H 2 production, and (3) the carbon and electron shuttles between species in coculture.

Ion Transport in the iMhu428 Metabolic Model
The majority of the iMhu428 model content comes from the iMB745 reconstruction of M. acetivorans and still needs to be verified. Changes to model content could affect the conclusions drawn about M. hungatei behavior. Despite containing a complete methanogenesis pathway, the iMhu428 model was unable to identify the H + /Na + transport stoichiometry of the energy-converting (Eha-or Ehb-type) hydrogenase (EHA, 1.12.7.2), which is thought to pump H + /Na + while reducing ferredoxin [7]. The heterodisulfide reductase (HDR, 1.8.98.1) can also reduce ferredoxin, and the iMhu428 model predicted HDR to be the only ferredoxin-reducing reaction required for methanogenesis. This observation is consistent with the observation that the expression of Eha/Ehb is considerably lower than that of HDR [57].
Additionally, different stoichiometries for other ion transport reactions important to methanogenesis remain thermodynamically possible. For example, the group contribution method predicted that tetrahydromethanopterin S-methyltransferase (MTSPCMMT_CM5HBCMT, E.C. 2.1.1.86) could drive transport of up to 4 Na + ions under standard conditions, instead of the 2 Na + ions used in the iMhu428 reconstruction. Furthermore, some studies suggest the archaeal A 1 A 0 ATP synthase is coupled to Na + instead of H + translocation [7,58]. When the iMhu428 model was modified to reflect this coupling, the model predicted the Na + /H + antiporter was no longer active, as Na + ions from MTSPCMMT_CM5HBCMT were directly used for ATP synthesis. Thus, while the modeled methanogenesis pathway is consistent with available data, it is not the only possibility.

H 2 Production in the iSfu648 Metabolic Model
Analysis of the iSfu648 model revealed that H 2 production is thermodynamically feasible in monoculture, implying there may be other biological reasons why H 2 production is not normally observed under this condition, or that tighter estimates of thermodynamic parameters are needed. The iSfu648 thermodynamic model also has some important limitations. In particular, the iSfu648 model does not contain enough thermodynamic information to predict the directions of important electron transport reactions that involve ferredoxin, including the confurcating hydrogenase and formate dehydrogenase, the ferredoxin-oxidizing hydrogenase and formate dehydrogenase, and the RNF-type oxidoreductase (S3 Table in S1 Dataset). This is because the group contribution method is unable to estimate the standard transformed Gibbs free energy of formation (Δ f G '0 ) of ferredoxin, resulting in no Δ r G '0 estimates for these reactions. Fortunately, new quantum chemical approaches for estimating the thermodynamics of metabolism [59] may potentially provide additional Δ f G '0 estimates. This work also raises important questions about the appropriate mathematical basis for representing thermodynamic constraints. Previous studies used the ΔG '0 of groups directly when modeling thermodynamics [34], and found that introducing uncertainty into a thermodynamic model of E. coli made the model computationally difficult to solve. In this work, using either the ΔG '0 of molecules or groups to model thermodynamics proved computationally difficult (results not shown). Instead, only using Δ r G '0 as the basis for thermodynamic calculations enabled uncertainties in free energy estimates to be handled without any computational difficulties. However, using Δ r G '0 as a basis for thermodynamic calculations leads to larger uncertainties in Δ r G '0 and greater network flexibility, as it does not account for thermodynamic interconnectivity between reactions with shared metabolites. As a result, the model-predicted feasible Δ r G '0 range through a linear combination of reactions considerably exceeds the group-contribution predicted Δ r G '0 range of the overall reaction. In the future, thermodynamic interconnectivity should be captured using the ΔG '0 of molecules or groups, and optimization techniques are needed to improve the runtime performance of the resulting thermodynamic models.
However, the use of Δ r G '0 as a basis for thermodynamic calculations is insufficient to explain why the iSfu648 model could still produce H 2 under monoculture growth conditions. As described in S1 Text, Δ f G '0 was used as a basis for thermodynamic calculations when attempting to find additional constraints which would prevent H 2 production. The failure to find such constraints indicates either that thermodynamics does not explain the absence of H 2 production in monoculture, or that current thermodynamic models cannot capture this phenomenon. If it is the latter, more accurate group contribution methods with smaller error estimates may eventually be able to explain the role of thermodynamics in syntrophic associations.

Formate and H 2 Transfer in Coculture
The thermodynamic coculture model of the syntrophic association between these species confirmed the role of formate and H 2 in electron transfer in the community, and led us to hypothesize that total H 2 consumption by the community indicates that M. hungatei cells are more abundant and/or faster growing than the S. fumaroxidans cells. The coculture model correctly predicted that both H 2 and formate could shuttle electrons between members of this community. Formate may be preferred over H 2 for electron transfer for thermodynamic reasons, as the Δ r G '0 of Eq 4 is more favorable (less positive) than that of Eq 3 in which formate is not exchanged. By exchanging formate in place of CO 2 and H 2 , S. fumaxoridans could sustain propionate oxidation at higher extracellular concentrations of formate than of H 2 . Formate exchange could also be preferred due to differences in kinetics, diffusion, and/or volatility. These scenarios could stabilize the syntrophic association by enabling faster shuttling of electrons to M. hungatei.
This work highlights some important obstacles to successful modeling of microbial consortia. In order for computational models of microbial communities to make meaningful predictions, individual species models must be integrated into a community model in a biologically relevant manner. Such integration will require an understanding of the objectives and constraints governing the behavior of each community member, and this work demonstrates that identifying the proper constraints and objectives requires extensive experimental characterization of the community. For example, neither maximization of growth rate nor maximization of ATP yield were sufficient for the iSfu648 model to predict the observed behavior of S. fumaroxidans. Additional constraints on reaction directions and flux ratios were required before iSfu648 could be combined with iMhu428 model to simulate the coculture. In addition, the iSfu648 model relied on data from gene expression and 13 C NMR experiments, suggesting that constraint-based approaches will complement traditional top down ('omics') approaches [60] by enabling a mechanistic understanding of microbial interactions [24].

Reconstruction of the iMhu428 Metabolic Model
The iMhu428 reconstruction of M. hungatei was built from the iMB745 reconstruction of M. acetivorans [41]. A preliminary draft reconstruction was built based on sequence homology (using the RAVEN Toolbox [42]), but the reconstruction contained less than 200 genes (S2 Dataset). Instead of performing extensive gapfilling, reactions from the iMB745 M. acetivorans reconstruction were copied into the M. hungatei reconstruction, with modifications to reflect key metabolic features of M. hungatei (see S1 Text). Results from the RAVEN Toolbox and the KEGG SSDB [43] were used to map genes in M. acetivorans to M. hungatei and identify those reactions which have genomic evidence (S2 Dataset). Finally, blocked reactions lacking genomic evidence were removed from the reconstruction.
The final iMhu428 reconstruction contains 720 reactions, 428 genes (associated with 493 reactions), and 639 metabolites. Of the 428 genes, 351 were added based on sequence homology, and 77 were added manually. The reconstruction is available in S2 Dataset and S3 Dataset in Excel and SBML formats.

Reconstruction of the iSfu648 Metabolic Model
The iSfu648 reconstruction of S. fumaroxidans was built from KEGG [43] (S6 Dataset and S7 Dataset) using the RAVEN Toolbox [42], which uses protein homology to identify the KEGG Orthology (KO) ID for each gene in a genome. The reactions and genes corresponding to that KO ID are then imported into the reconstruction. The resulting draft reconstruction was manually refined as described in S1 Text. The final iSfu648 reconstruction contains 874 reactions, 648 genes (associated with 770 reactions), and 893 metabolites. The reconstruction is available in S4 Dataset and S5 Dataset, in Excel and SBML formats.

Thermodynamics-Based Metabolic Flux Analysis (TMFA)
Flux-balance analysis (FBA) [27] is a constraint-based technique for predicting the state of a metabolic network consistent with physiochemical principles. FBA identifies a flux distribution which maximizes cellular growth (or some other objective function), subject to steady-state mass-balance and enzyme capacity constraints. Thermodynamics-Based Metabolic Flux Analysis (TMFA, [33,34]) extends FBA via the introduction of thermodynamic constraints, which require that the transformed Gibbs free energy of a reaction (Δ r G ' ) and its flux (v) have opposite signs. Estimates (D r G 0 0 est ) and uncertainties (SE D r G 0 0 est ) of Δ r G '0 for the reactions in the reconstructions were obtained using a group contribution method [61] via the von Bertalanffy 2.0 Toolbox [62]. TMFA was implemented as previously described [34], with additional details given in S1 Text. The mol files for metabolites in iMhu428 and iSfu648 are provided in S1 File and S2 File, respectively.

Parsimonious TMFA (pTMFA)
pFBA [54] is a constraint-based approach which maximizes cellular growth while also minimizing total flux through the network (a proxy for minimizing the total mass of enzymes required to sustain optimal growth through the network). pTMFA uses the same assumptions as pFBA while implementing the thermodynamic constraints of TMFA. pTMFA was implemented as a two-stage optimization process. In the first stage, growth rate is maximized via TMFA. In the second stage, the growth rate is fixed and the total flux through the network is minimized, subject to the same constraints as TMFA. Additional details on implementation can be found in S1 Text.

Coculture Model
For the coculture simulations, a community model of growth in a continuous stirred-tank reactor was developed that accounts for the biomass concentrations of each species. The model minimizes the species-weighted total flux through the metabolic networks subject to TMFA constraints for each species and mass balances around the entire reactor. Details on the specific implementation used in this work can be found in S1 Text. To avoid solving a mixed-integer non-linear program (MINLP), the dilution rate and biomass concentrations for each species were fixed, resulting in a MIP. To explore the community behavior under a variety of operating conditions, the reactor dilution rate was systematically changed, while allowing unlimited propionate uptake by the reactor.

Simulations
The uptake fluxes for carbon and other nutrients used in the simulations are given in S1 Table in  Supporting Information S1 Dataset. Supporting tables. This file contains Supporting Tables S1 to S7 as described in the text.