A full multiscale model of the metabolism of the human gut would need to include around 1014 individual bacteria belonging to hundreds of bacterial species, for which in many cases curated GEMs are unavailable. We thus necessarily resorted to a more coarse-grained approach, while maintaining some level of biological realism by constructing the model based on a validated, genome-scale metabolic network model of Lactobacillus plantarum [28]. Figure 1 gives an overview of the workflow of the paper. We first (1) constructed a metabolic model representing a subset of the gut microbiota, which we used for the dFBA model (2). We then asked to what extent cross-feeding can emerge in large communities of interacting and diversified bacteria, such as those found in the colon, using a dynamic multi-species metabolic modeling (DMMM) approach [18, 23, 29], which is an extension of the dynamic flux-balance analysis (dFBA) method [22, 30]. To this end, we constructed a well-mixed model of a bacterial consortium (3), by coupling 1000 of the dFBA models via a common, external exchange medium that allowed the bacteria to exchange a subset of the metabolites in the GEM. We initiated the exchange medium with a pulse of glucose, then observed the turn-over of glucose into a series of short-chain fatty acids (4), and quantified cross-feeding (5): the extent to which the bacteria exchanged metabolites via the common medium. Next we asked to what extent spatially diversified microbial communities can emerge in a tube-like environment (6), if the microbial communities are allowed to specialize to the local availability of metabolites. In the spatial model, the GEMs inside the bacteria were allowed to evolve. After running the model for a fixed time, we quantified how much the GEMs had diversified and performed local cross-feeding (7) and to what extent they had locally changed the external concentrations of metabolites (8), leading to stratification and niche formation.
Construction of a metabolic model representing a subset of the gut microbiota
We first constructed a hypothetical, but biologically-realistic “supra-organism” model [3, 31], called “metabacterium” here, that represents a sample of the gut microbial community in a single metabolic network model. For this preliminary, explorative study we used a GEM of Lactobacillus plantarum [28], a resident of the colon and a strain widely used for probiotics, and extended it with four key metabolic pathways of the intestinal microbial community: (1) propionate fermentation, (2) butyrate fermentation, (3) the acrylate pathway and (4) the Wood-Ljungdahl pathway. In future versions of our framework this network could be replaced by metabolic network models derived from metagenomic data [3] as they become available. The current, simplified network contains 674 reactions (Supplementary File 1), and compares well with consensus metabolic networks of carbohydrate fermentation in the colon [32, 33]. For a schematic overview of the key pathways including in the metabolic network, see Fig. 2
a.
The uptake and excretion rates of genome-scale metabolic networks can be calculated using constraint-based modeling. To represent diauxic growth, i.e., by-product secretion as a function of extracellular metabolite concentrations, we used an extension of FBA called Flux Balance Analysis with Molecular Crowding (FBAwMC) [34]. FBAwMC correctly predicts diauxic growth and the associated secretion of by-products in micro-organisms including E. coli, Saccharomyces cerevisiae [35], and L.plantarum [36]. As an additional, physiologically-plausible constraint FBAwMC assumes that only a finite number of metabolic enzymes fits into a cell, with each enzyme having a maximum metabolic turnover, V
max. For each reaction, FBAwMC requires a crowding coefficient, defined as the enzymatic volume needed to reach unit flux through that reaction. Each reaction is assigned a “crowding coefficient”, a measure of the protein cost of a reaction: Enzymes with low crowding coefficients have small molecular volume or catalyse fast reactions. Given a set of maximum input fluxes, FBAwMC predicts the optimal uptake and excretion fluxes as a function of the extracellular metabolite concentrations.
As FBAwMC optimizes growth rate, not growth yield as in standard FBA, it predicts a switch to glycolytic metabolism at high glucose concentrations at which faster metabolism is obtained with suboptimal yield. Its accurate prediction of diauxic growth together with by-product secretion as a function of extracellular metabolite concentrations make FBAwMC a suitable method for a microbial community model.
Metabolic diversity causes cross-feeding in a well-mixed system
To study the extent of cross-feeding emerging already from a non-evolving metabolic community of “metabacteria”, we first set up a simulation of 1000 interacting metapopulations, where each subpopulation was initiated with a set of crowding coefficients selected at random from an experimentally determined distribution of crowding coefficients of Escherichia coli [35, 36], for lack of similar data sets for L. plantarum. The simulation was initiated with pure glucose and was ran under anaerobic conditions. We then performed FBAwMC on all 1000 metapopulations, optimizing for ATP production rate as a proxy for growth rate. This yielded 1000 sets of metabolic input and output fluxes, F
i
, and growth rates, μ
i
for all 1000 metapopulations. These were used to update the extracellular concentrations, \(\vec {M}\) and metapopulation sizes, X
i
, by performing a single finite-difference step of [23, 29]
$$ \frac{d\vec{M}}{dt}=\sum_{i}X_{i}\vec{F_{i}} $$
(1)
and
$$ \frac{dX_{i}}{dt}=\mu_{i}X_{i}. $$
(2)
with a timestep Δ
t=0.1 h. After updating the environment in this way, we performed a next time simulation step.
Figure 2
b shows how, in the simulation, the metabacteria modified the environment over time. The secondary metabolites that were produced mostly are acetate, butyrate, carbon dioxide, ethanol, formate, lactate and propionate. This compares well with the metabolites that are actually found in the colon [37] or in an in vitro model of the colon [38]. In the first 30 min of the simulation, the initial pulse of glucose is consumed, and turned over into acetate (red), lactate (grey), formate (brown), and ethanol (yellow). These are then consumed again, and turned over into proprionate (purple) via pathways 3 and 4 (Fig. 2
a) and into butyrate (blue) via pathways 7 and 11. CO
2 is also increasing due to the turnover of pyruvate into acetyl co-A via pathway 5 (pyruvate dehydrogenase). After about two hours of simulated time, proprionate and CO
2 levels drop again due to the production of butyrate (blue): proprionate is consumed reversing reaction 3 and 4; CO
2 is consumed in pathway 9 that produces acetate from formate. The conversion of acetate back to acetyl-coA then drives the production of butyrate; a surplus of acetyl-coA is turned over into acetaldehyde and ethanol in pathway 10. Interestingly, formate and CO
2 are produced at the same time; this rarely occurs in any single organism but does occur in this microbial consortium.
To test to what extent these results depend on the ability of the individual FBAwMC models to represent metabolic switching and overflow metabolism [34, 36], we also simulated the model using standard flux-balance analysis [39]. In this case, all glucose was converted into ethanol, whereas lactate and propionate did not appear in the simulation (Additional file 1: Figure S1). To test to what extent the results rely on cross-feeding, we also checked if any of the single-species simulations could also produce so many metabolites. Out of 100 single-species simulations none produced as many or more excreted metabolite species than the interacting set of species.
Quantification of cross-feeding
Most of the metabolites were only transiently present in the medium, \(\vec {M}\), suggesting that the metabolites were re-absorbed and processed further by the bacteria. To quantify the amount of such cross-feeding in the simulations, we defined a cross-feeding factor, C(i), with i a species identifier. Let
$$\begin{array}{@{}rcl@{}} F_{\mathrm{up,tot}}(i,j)&\equiv& \int_{t=0}^{t_{\text{max}}} B(n,t)F_{\text{up}}(i,j,t)dt\\ F_{\mathrm{ex,tot}}(i,j)&\equiv& \int_{t=0}^{t_{\text{max}}} B(n,t)F_{\text{ex}}(i,j,t)dt \end{array} $$
(3)
be the total amount of metabolite j that species i consumes and excretes during the simulation. B(i,t) here equals the biomass of species i at time t. The amount of carbon species i gets via cross-feeding then equals,
$$ \begin{aligned} C(i)&=\sum\limits_{j}c_{\mathrm{C}}(j)\text{max}(F_{\mathrm{up,tot}}(i,j)-F_{\mathrm{ex,tot}}(i,j),0)\\ &\quad-6F_{\mathrm{up,tot}}(i,\text{glucose}). \end{aligned} $$
(4)
Here, c
C(j) is the molar amount of carbon atoms per mol metabolite j (e.g., c
C(glucose)=6). If species i during the fermentation consumes more of metabolite j than it has produced, species i has cross-fed on metabolite j. We subtract the amount of glucose from the sum, because glucose is the primary nutrient source that is present at the start of the simulation. Now we can calculate the total amount of carbon the population acquires via cross-feeding, relative to the total amount of carbon taken up by the population
$$ C_{\text{rel}}=\frac{\sum_{i}C(i)}{\sum_{i}\sum_{j} c_{\mathrm{C}}(j)F_{\mathrm{up,tot}}(i,j)}. $$
(5)
If C
rel=0, there is no cross-feeding. In that case, every species only consumes glucose as carbon source or only consumes as much carbon from other metabolites as it has secreted itself. Conversely, if C
rel=1 all carbon the species has consumed during the simulation is from non-glucose carbon sources the species has excreted itself. For the whole simulation C
rel=0.39±0.02, indicating that 39% of all carbon consumed by the bacteria comes from cross-feeding. Cross-feeding was largest on lactate, CO
2, acetate, ethanol, formate and propionate. Many of these metabolites are known to be involved in bacterial cross-feeding in the colon or cecum (for interconversion between acetate and lactate, see Ref. [40]; and for interconversion between acetate and butyrate in the murine cecum, see Ref. [41]). In the original L. plantarum model we also find cross-feeding, but only on lactate and acetaldehyde (Additional file 2: Figure S2). Taken together, in agreement with previous computational studies that showed cross-feeding in pairs of interacting E. coli [23], these simulations show that cross-feeding interactions occur in coupled dynamic FBAwMC models.
Spatially explicit, evolutionary model
The well-mixed simulations showed that cross-feeding appears in populations of interacting metabacterial metabolic networks. However, this does not necessarily imply microbial diversity, because it is possible that the same metabacterium secretes and reabsorbs the same metabolites into the substrate, in which case there would be no true cross-feeding. Furthermore, the previous section did not make clear whether cross-feeding will be ecologically stable under conditions where subpopulations of the supra-organisms are lost. In a spatially explicit model, cross-feeding possibly arises more easily and is more easy to detect, as different metabolic functions can be performed at different locations [42]. We therefore developed a spatially explicit, multiscale evolutionary model of gut microbial metabolism. We initiate the simulation with a population of metapopulations of bacteria that can perform all metabolic functions under anaerobic conditions, just as in the well-mixed simulation. We then let the systems evolve and study if meta-populations of bacteria with specific metabolic roles evolve.
Model description
Figure 3 sketches the overall structure of our model. The model approximates the colon as the cross-section of a 150 cm long tube with a diameter of 10 cm. The tube is subdivided into patches of 1 cm
2, each containing a uniform concentration of metabolites, and potentially a metapopulation of gut bacteria (hereafter called “metabacterium”) (Fig. 3
a). Each metabacterium represents a small subpopulation (or ’metapopulation’) of gut bacteria with diverse metabolic functions, and is modeled using a metabolic network model containing the main metabolic reactions found in the gut microbiota, as described above (Fig. 2
a). Based on the local metabolite concentrations, \(\vec {c}(\vec {x},t)\), the metabolic model delivers a set of exchange fluxes F
i,n and a growth rate, \(\mu (\vec {x})\), which is assumed to depend on the ATP production rate (Fig. 3
b; see “Methods” for detail). The metabolites disperse to adjacent patches due to local mixing, which we approximate by a diffusion process (Fig. 3
c), yielding
$$ \frac{d\vec{c}(\vec{x},t)}{d t}=\vec{F}(\vec{x},t)B(\vec{x},t) + \frac{D}{L^{2}} \sum\limits_{\vec{i}\in\text{NB}(\vec{x})} \left(\vec{c}(\vec{i},t)-\vec{c}(\vec{x},t)\right), $$
(6)
where \(\vec {F}(\vec {x},t)\) is the flux of metabolites between the medium and the metabacterium, and the sum runs over the four nearest neighbors \(\text {NB}(\vec {x})\); dispersion is approximated by Fick’s law, where D is a diffusion coefficient and L=1 cm the interface length between two adjacent patches. The local density of metabacteria, \(B(\vec {x})\) is given by
$$ \frac{d B(\vec{x},t)}{dt}=\mu(\vec{x},t)B(\vec{x},t). $$
(7)
To mimic meals, a pulse of glucose of variable magnitude enters the tube once every eight hours (Fig. 3
d). The metabolites move through the tube via a simplified model of advection: At regular intervals, all metabolites are shifted one patch (Fig. 3
e). Metabolites continuously leave the tube at the end through an open boundary condition (Fig. 3
f). To mimic peristaltic movements that locally mix the gut contents together, metabacteria randomly hop to adjacent lattice sites (Fig. 3
g) and leave the gut only via random hops over the open boundary condition. In a subset of simulations, accelerated bowel movements are simulated by advecting the metabacteria together with the metabolites. To a first approximation, the boundaries are impermeable to the metabolites, a situation reflecting many in vitro models of the gut microbiota (reviewed in Ref. [43]); later versions of the model will consider more complex boundary conditions including absorption of metabolites [44].
When the local biomass in a patch, \(B(\vec {x},t)\), has grown to twice its original value, the metapopulation expands into the second position on the grid point (Fig. 3
h). To mimic a local carrying capacity, the metapopulation does not spread out or grow any further if both positions in the patch are occupied. In the visualizations of the simulations, full patches are shown in green, singly occupied patches are shown in yellow, and empty patches are shown in black (Figs. 3
i and 4). During expansion, changes in the relative abundance of species may enhance or reduce the rate of particular reactions, or even delete them from the metapopulation completely. Similarly, metabolic reactions can be reintroduced due to resettling of metabolic species, e.g., from the gut wall mucus [45]. To mimic such changes in species composition of the metapopulation, during each expansion step, we delete enzymes from the metabolic network at random, reactivate enzymes at random, or randomly change crowding coefficients such that the metapopulation can specialize on one particular reaction or become a generalist.
The crowding coefficients, as they appear in the flux-balance analysis with molecular crowding (FBAwMC) method that we used for this model, give the minimum cellular volume filled with enzymes required to generate a unit metabolic flux; they are given by the V
max of the enzyme and enzyme volume [34]. Equivalently, in our metapopulation model, the crowding coefficient of a reaction is the minimum intracellular volume averaged over all bacteria in the patch that must be filled with enzymes in order to generate a unit flux through the reaction. It depends on the density of the enzyme in the bacteria and also on the corresponding values of V
max. Because the V
max of a reaction can differ orders of magnitudes between species (see for example the enzyme database BRENDA [46]), the evolutionary dynamics in our model could drive the metabacteria to reduce all crowding coefficients concurrently, producing a highly efficient generalist. To introduce a biologically more realistic trade-off between metabolic rate and cost in terms of volume, we therefore included an experimentally observed trade-off between growth rate and growth yield among micro-organisms [47, 48]: Micro-organisms that grow fast have low growth yield and vice versa. We take this trade-off into account explicitly by assuming a maximal growth rate given the carbon uptake rate of the cells. This trade-off prevents the metabacteria from growing infinitely fast by mutating their crowding coefficients.
As an initial condition, we distribute metabacteria over the grid, each containing all available metabolic reactions, i.e., each metabacterium initially contains all bacterial “species” that the complete metabacterium represents. To reflect variability in the relative abundances of the bacterial species in each metabacterium the crowding coefficients are drawn at random from an experimental distribution as described above (Fig. 3
a).
Evolution of diversity due to metabolic cross-feeding
To evaluate the behavior of our model, we performed ten independent simulations. These show largely similar phenomenology; therefore we will first describe the progression of one representative simulation in detail, and then discuss differences with the other simulations. Figure 5
a shows the average number of metabolic reactions present in the metabacteria over time in the simulation. At t=0 all metabacteria still have all 674 reactions, but over time the number of available reactions gradually drops to below 200. This reduction of the number of metabolic genes could indicate a homogeneous population that is specialized, e.g., on fermentation of glucose where most of the metabolic network is not used. An alternative explanation is that each of the metapopulation retains a subset of the full network, an indication of cross-feeding. The amount of cross-feeding will likely change over the tube: The metabacteria in the front have direct access to glucose, whereas the metabacteria further down in the tube may rely on the waste-products of those in front. We therefore determined a temporal average of the cross-feeding factors, C
rel
(Eq. 5), at each position in the tube over t=3500 to t=4000, a time range at which most genes have been lost. The first observation to note is that in the spatial evolutionary simulations, the average cross-feeding factor C
rel has a higher value than in the well-mixed simulations. In this particular simulation, the spatial average cross-feeding factor at t=4000 is C
rel=0.65±0.09, compared with C
rel=0.39±0.02 in the well-mixed case (n=10). The cross-feeding factor for individual cells (C(i), Eq. 4), showed large population heterogeneity. As Fig. 5
b shows, the cross-feeding factor in the tube front is close to 0, indicating the presence of primary glucose consumers, while cross-feeding slowly increases towards the distal end until it almost reaches 1, indicating complete cross-feeding. Thus in the proximal end the bacteria rely mostly on the primary nutrient source, while near the distal end cells of the tube rely on cross-feeding. This observation is consistent for all simulations (see Additional file 3: Figure S3).
Emergence of metabolic stratification
We next investigated the mechanism by which such cross-feeding emerges in the simulation. Additional file 4: Figure S4 plots the metabolite concentrations over evolutionary time for the simulation of Fig. 5. In this particular simulation, the concentrations of formate and lactate initially rise rapidly, after which they drop gradually. The butyrate concentrations increase over evolutionary time. In all simulations, the metabolite concentrations change gradually, but not necessarily following the same temporal pattern.
Figure 6 shows the spatial distribution of a set of key metabolites averaged over 2000 h to 4000 h of the representative simulation. Interestingly, the flow of metabolites through the colon in interaction with the bacterial population creates a spatially structured, metabolic environment. The proximal end is dominated by the primary carbon source glucose (Fig. 6
a), with the peak in the average glucose concentration due to the periodic glucose input. Further down in the tube we find fermentation products, including lactate and ethanol, whereas the distal end contains high levels of acetate and CO
2, showing that the metabacteria convert the glucose into secondary metabolites. Among these secondary metabolites, the levels of acetate (Fig. 6
b), ethanol (Fig. 6
e), formate (Fig. 6
f), lactate (Fig. 6
g) and propionate (Fig. 6
h) drop towards the distal end off the tube, so they are further metabolized by the metabacteria. In this particular simulation, butyrate and CO2 are not consumed and their concentrations increase monotonically towards the end. The small drop at the very distal end is caused by the metabolite outflow. The profiles of the other simulations were consistent with this representative simulation (Additional file 5: Figure S5). In all simulations, the proximal end was dominated by glucose. Further towards the end of the tube, zones of fermentation products developed as in the representative simulation, but the precise location of each product was different and not all products were present. Most notably, in two out of ten simulations, butyrate was absent and in two other simulations proprionate was absent. Also, in three out of ten simulations lactate was more confined to the front of the tube (up to around 50 sites) than in the representative simulation.
Metabacteria specialize on local metabolic niches
These results demonstrate that the metabacteria spatially structure their metabolic environment, generating a stratified structure of metabolic “niches” along the tube, each offering a separate set of metabolites. Therefore, we next asked if this environmental structuring gives rise to metapopulations uniquely adopted to the microenvironment. We took computational samples of all metabacteria found in the tube between 3500 h and 4000 h, to average out the variations at the short timescale. We tested the growth rate of these samples (consisting of on average n≈1100 metabacteria) in six, homogeneous metabolic environments, containing uniform concentrations of pure (1) glucose, (2) acetate, (3) formate, (4) lactate, and (5) propionate, and (6) a mixture of of CO
2 and H
2. Figure 7 shows the average and standard deviation of the growth rates of the metabacteria in each of these six environments, as a function of the position from which they were sampled from the tube. Strikingly, the metabacteria near the distal end of the tube have lost their ability to grow on glucose (Fig. 7
a), indicating that they have specialized on secondary metabolites, including acetate (Fig. 7
b) and lactate (Fig. 7
e). Interestingly, in support of the conclusion that the metabacteria specialize on the metabolic niches generated by the population as a whole, the metabacteria sampled from the distal end on average grow faster on acetate and lactate than the metabacteria sampled from the front of the tube. Acetate and lactate are produced in the proximal colon and flow to the distal part of the tube where the metabacteria can metabolize it; in the front of the tube acetate and lactate concentrations are lower, such that neutral drift effects can safely remove the corresponding metabolic pathways from the metabacteria. Remarkably, the metabacteria also grow on CO
2, because of the presence of hydrogen gas, that allows growth on CO2 via the Wood-Ljungdahl pathway [49]. To further characterize the alternative metabolic modes occurring in the model, we clustered the population present at the end of the simulation t = 4000 h with respect to their maximum growth rates in the six environments (Fig. 8). Clearly, different metabolic “species” can be distinguished. One “species” can metabolize glucose, a second “species” can metabolize most secondary metabolites and a third “species” has specialized on acetate. Thus in our simulation model a number of functional classes appear along the tube, each specializing on its own niche in the full metabolic network.
Increased flux through the tube makes diversity collapse
From the results in the previous section, we conclude that the inherent spatial structuring of the colon results in separate niches. This allows the population to diversify, such that different “species” have different metabolic tasks. A recent population-wide metagenomics study of stool samples from the Flemish and Dutch population [50] showed that, among a range of life-style related factors and medicine use, the diversity of the human gut microbiota correlates strongest with the Bristol stool scale (BSS), a self-assessed indicator of the “softness” of the stool. The analysis showed that for softer stools (higher stool index, indicative of faster transit times [51]), the diversity of the gut microbiota was reduced [52]. To investigate whether transit time could also be correlated with reduced diversity in our model, we studied the effect of increased fluxes through the tube (“diarrhea”), by assuming that the supra-bacteria flow through the tube at the same rate as the metabolites do. Strikingly, the maximal growth rate of the cells has become independent of the position (Fig. 9). Again, we clustered the population present at the end of the simulation with respect to their maximum growth rates in glucose, acetate, H2 and CO2, formate, lactate and propionate (Fig. 10). In contrast to the simulations without cell flow, the population does practically not diversify. All supra-bacteria can grow on glucose, acetate and H2 and CO2. Thus, our simulations suggest that increased transit times may contribute to a reduction of microbial diversity, by reducing the spatial heterogeneity in the gut and, consequently, the construction of ecological niches and cross-feeding interactions.