Supplemental Materials
Table of Contents
Dynamic Multi-species Metabolic Modeling framework 2
The formulation of the DMMM framework 2
Sample Implementation of the DMMM framework 5
Advantage of the DMMM framework 8
DMMM vs. Monod-kinetic Model 8
DMMM vs. Stolyar Model 9
DMMM vs. Experimental Studies 10
Modeling Subsurface Iron-Reducing Community 12
Using DMMM to model Rhodoferax and Geobacter Species 12
Nutrient Uptake Parameters 13
Simulation of Cell-Death 16
List of Simulation Parameters 17
Acetate Uptake Parameters 17
Ammonium Uptake Parameters 17
Fe(III) Uptake Parameters 18
Simulations under Natural Conditions 18
Simulations During Acetate Addition 19
Parameter Sensitivity Analysis 20
Ammonium Affinity 20
Rhodoferax’s Fe(III) Affinity 20
Rhodoferax’s Acetate Affinity 21
List of Assumptions 22
DMMM and the Prediction of Metabolic States 25
Supplemental Figures 26
References 28
Dynamic Multi-species Metabolic Modeling framework
The formulation of the DMMM framework
A community metabolic model must account for the metabolic exchanges between species and with the environment, as well as the changes in biomass of the modeled species. The Dynamic Multi-species Metabolic Modeling (DMMM) framework extends the dynamic flux balance analysis (dFBA) formulation to the community realm. A mathematical formulation of the DMMM framework describes a community of N microbial species coexisting in an environment containing MEX metabolites.
The growth rate (dX/dt) of every microbial species in the community is given by:
(1)
The consumption/production rate (dS/dt) of every metabolite in the environment is given by:
(2)
Here, Xj is the biomass of the jth species in the community; Si is the concentration of the ith metabolite in the environment. The specific growth rate (μj) of the jth microbial species, and , the specific consumption/production rate of the ith metabolite due to the actions of the jth microbial species, are calculated using the FBA:
(3)
In Equation Set 3, cT is the objective function; in our simulations, growth maximization is used as the objective function. Aj is the stoichiometric matrix of the jth microbial species. vj is the reaction flux vector of the jth microbial species. vjmax and vjmin are the flux capacity constraints of the jth microbial species based on the corresponding genome-scale model. For the components of vjmax and vjmin corresponding to external metabolites, the uptake/production constraints (Vjmax,Vjmin) to these fluxes can be calculated based on the environmental concentration of these metabolites, using either the On/Off method (if the uptake kinetics are not available) or the Michaelis–Menten kinetics method (if kinetics are available) . It is important to note that a separate FBA model is used for each microbial species in the community – N FBA models are required for a community of N member species. The flux constraints and the reaction network itself forms a viable solution space, and the solution in this space where the objective function (specific growth rate) is maximal is selected by linear programming. If no viable optimal solution is found for species j, then it is concluded that species j does not have sufficient nutrients to survive at this time, and Equation 8 is used to calculate the death rate (as a negative specific growth rate).
By integrating the growth rates of all microbial species within the community, as well as the production/consumption rates of all metabolic species in the environment, the DMMM framework can dynamically predict the temporal changes in metabolite and biomass concentrations in a complex microbial community.
A small sample implementation of the DMMM framework is included in the next section to further illustrate the simulation process. The DMMM framework is implemented in MATLAB as an add-on to the COBRA toolbox (Becker et al., 2007). This allows the user to easily incorporate new species into the DMMM framework. The MATLAB code for the DMMM framework can be downloaded from Supplemental Materials.
Sample Implementation of the DMMM framework
A sample implementation of the DMMM framework is provided here to clarify the algorithm. This sample implementation assumes the existences of a community consisting of two microbial species, A and B, and three substrates, Met1, Met2, Met3 (Figure S1). The numbers of microbial and metabolic species are arbitrary in this example. The maximal numbers of microbial and metabolic species depend on the limitation of the computational hardware. After a simulation length (hrs) and a simulation step size (hrs) is given, this particular implementation of the DMMM framework execute the following routines during each simulation step:
R1: The concentration of Met1, Met2, and Met3 are used to calculate the upper and lower constraints to V1A, V2A, V3A, V1B, V2B, and V3B, where Vij represent the ith reaction flux of the species j. Two calculation methods are commonly used:
i. On/Off methods
if [Meti] > 0, then -∞< Vij <∞
if [Meti] <= 0, then Vij =0
j. Michaelis–Menten kinetics
R2: The FBA models of species A and B calculates the specific growth rates, μA and μB, as well as the reaction fluxes V1A, V2A, V3A, V1B, V2B, and V3B. The models can be described by the following equations.
The constraints to the reaction fluxes (including the uptake constraints) and the reaction network itself forms a viable solution space, and the point in this space where specific growth rate is maximized is selected using linear programming. If the linear programming method provides a solution, then continue to routine R3. If the linear program method fails for species j (indicating that the cells does not have the required nutrients to survive under the conditions at this time), then a special cell-death routine RD is executed to calculate the rate of death rd of the species j, and μj is set to rd.
R3: The rate of change in the biomass of A and B are calculated with the equation:
Calculate ∆Xj by integrating dXj/dt over the simulation step.
Calculate the new Xj at the end of the simulation step with Xj new = Xj old +∆Xj.
R4: The rate of change in the concentration of external metabolite i is with the equation:
Calculate ∆Si by integrating dSi/dt over the simulation step.
Calculate the new Si at the end of the simulation step with Sinew = Siold +∆Si.
Only the rates of change of external metabolites are integrated since all internal metabolites still follow the internal steady state assumption.
RD: The cell death RD routine is a special user-defined routine that is called upon when the environmental concentrations of substrates are insufficient to sustain the current biomass concentration of the organism. This routine is flexible and can be redefined for each individual organism to reflect its specific mechanism of cell-decay. The RD routine produces the rate of cell death rd, which should be considered a negative growth rate. The specific method used to calculate rd is described in the main text. After completion of RD, return to R3.
After routine R4 is completed, time advances one step-size. The routine R1 is initiated again. This continues until the simulation length is reached. At each time point, the flux constraints to each organism varies base on the substrate concentration at that particular time, and leading to dynamic variations in the growth rate. This procedure is illustrated in Figure S1.
Advantage of the DMMM framework
The usage of the DMMM framework to model microbial communities holds significant advantages over the traditional Monod-kinetic models, the Stolyar community metabolic model (Stolyar et al., 2007) mentioned in the main text, as well as pure experimental methods.
DMMM vs. Monod-kinetic Model
The Monod-kinetic model can model relatively simple metabolism quite well. For example, if an organism is treated as having only one limiting substrate, S, then the organism’s growth with respect to S can be modeled fairly accurately by:
If an organism is treated as having two limiting substrates, S1 and S2, then the organism’s growth can be modeled with:
For an organism with i limiting substrates, the organism’s growth is modeled with the General Monod Model:
For example, if we assume that each of the substrates are available at their saturation concentration (Ks), the Monod model would predict a growth rate of 1/2n of the maximum, even though it is clear that the substrate concentrations are sufficient to allow the organism to grow at ½ the maximum. Moreover, the yield in the presence of multiple substrates needs to be experimentally determined for each combination, as there are no known mechanisms to compute the yields accurately for different metabolic states. Since n is generally related to the metabolic complexity, the Monod-kinetic model does not scale up well with complex metabolisms. Realistically, prediction accuracy often becomes unacceptable when i > 2. Since the DMMM framework is based on constraint-based models, it is capable of handling significantly more complex metabolisms. For example, the Geobacter sulfurreducens model used in this paper contains 727 reactions, 55 of which are exchange reactions (i = 55). It is also important to note that the metabolic network of Geobacter has many modes of operation such acetate-limiting ammonium-utilization mode, iron-limiting ammonium-utilization mode, acetate-limiting nitrogen fixation mode, iron-limiting nitrogen fixation mode. The DMMM framework automatically selects the modes of operations based on the substrate concentrations at each time point. This is unachievable for traditional models.
DMMM vs. Stolyar Model
Stolyar et al. published the first constraint-based metabolic model of a mutualistic microbial community in 2007 (Stolyar et al., 2007). In the Stolyar model, the CBMs of Desulfovibrio vulgaris and Methanococcus maripaludis are directly connected. The authors assumed that the species are interdependent, thereby justifying their usage of a pair of arbitrary biomass flux weights (authors used 10:1, 1:1, 1:10) as the objective function. This objective function is inappropriate for most microbial communities. Another caveat of this model is that it cannot predict the dynamic shifts in population and their metabolite concentrations. While this approach may be appropriate when the microorganisms are inter-dependent, it is inappropriate in ecological settings where the community composition is dynamic.
In comparison, the DMMM framework does not rely on arbitrary objective function – instead, the FBA problems representing each organisms are solved separately, allowing the usage of more established objective functions such as maximization of biomass flux. The DMMM framework is not restricted to mutualistic communities; it is capable of predicting the dynamic shifts in community composition and metabolites’ concentration.
DMMM vs. Experimental Studies
While no model can ever achieve the accuracy of real experimental studies, usage of metabolic modeling has proved to be a valuable analytical tool. The DMMM framework allows scientists to consider thousands of possible reactions and metabolites simultaneously, this together with necessary experimental work, can provide significant insights that are otherwise overlooked due to the limitation on the processing capacity of the human brain. Furthermore, the DMMM framework allows scientists to observe the emergent properties of the complex microbial community, which are difficult (if not impossible) to predict without.
The DMMM framework allows complex experiments to be performed in silico prior to real experiments. These simulations can guide the formation of new hypothesis, which allows in situ and in vivo experiments to be performed selectively. Furthermore, this approach allows scientists to test otherwise difficult-to-prove hypothesis. For example, a scientist designing an environmental biotechnology strategy might want to study the effect of increasing groundwater flow rate, which would be very difficult to evaluate in situ.
Modeling Subsurface Iron-Reducing Community
Using DMMM to model Rhodoferax and Geobacter Species
The DMMM framework was used to evaluate competition between Geobacter and Rhodoferax species in anoxic environments in which Fe(III) is the primary electron acceptor available to these organisms (Figure 2). This represents conditions in uranium-contaminated subsurface environments, because Fe(III) is expected to be the predominant electron acceptor under anoxic conditions (Lovley, 1991) and the concentration of U(VI) contaminants in the groundwater is not significant enough to sustain growth as the sole acceptor (Finneran et al., 2002). While the concentration profiles of all metabolites that are either secreted into or consumed from the environment and the biomass of the organisms in the community can be calculated, our study focuses on acetate, and Fe (III), which are the electron donor and acceptors of these two organisms, as well as ammonium, which can significantly affect growth yield (Figure 2).
Geobacter sulfurreducens and Rhodoferax ferrireducens served as models for Geobacter and Rhodoferax species because genome-scale metabolic models of these organisms are available (Mahadevan et al., 2006; Risso et al., 2009; Sun et al., 2009). These models have been used to successfully predict the biomass yield of these two species under multiple physiological conditions (Mahadevan et al., 2006; Risso et al., 2009; Sun et al., 2009; Segura et al., 2008). The energetic of G. sulfurreducens and R. ferrireducens has been discussed in detail (Mahadevan et al., 2006; Risso et al., 2009). The ATP maintenance requirement value for G. sulfurreducens and R. ferrireducens was set to 0.45 mmol ATP/gDW (Mahadevan et al., 2006; Risso et al., 2009; Sun et al., 2009).The simulations were performed using the maximization of biomass yield as a cellular objective (Varma & Palsson, 1994; Varma & Palsson, 1995). As in the dynamic FBA (dFBA) formulation (Mahadevan et al., 2002), the nutrient uptake/production flux constraints (mmol/gDW/hr) used in the FBA models of the individual species are calculated using the Michaelis-Menten expressions.
Nutrient Uptake Parameters
The acetate transport kinetics of G. sulfurreducens have been determined by measuring [14C]-labeled acetate uptake (Richter et al., Submitted). Three different acetate uptake mechanisms, each with different saturation constants, Ks (mM), and specific uptake rates, Vmax, were identified. All three of these uptake systems were observed to be up-regulated under acetate-limited growth, suggesting that they are co-expressed (Risso et al., 2008). Hence, the overall acetate consumption rate of G. sulfurreducens was calculated using the experimentally identified uptake parameters for these uptake systems (Equation 4).