Integrated culturing, modeling and transcriptomics uncovers complex interactions and emergent behavior in a synthetic gut community

Whereas the composition of the human gut microbiome is relatively well resolved, predictive understanding of its response to perturbations such as diet shifts is still lacking. Here, we followed a bottom-up strategy to explore human gut community dynamics. We established a synthetic community composed of three representative human gut isolates in well-controlled conditions in vitro. We then explored species interactions by performing all mono- and pair-wise fermentation experiments and quantified with a mechanistic community model how well tri-culture dynamics was predicted from mono-culture data. With the model as a reference, we demonstrated that species grown in co-culture behaved differently than in mono-culture and confirmed their altered behavior at the transcriptional level. In addition, we showed with replicate tri-cultures and in simulations that dominance in tri-culture sensitively depends on initial conditions. Our work has important implications for gut microbial community modeling as well as ecological interaction detection from batch cultures.


Introduction
The human gut microbiome is a complex, spatially heterogeneous and dynamic system consisting of hundreds of species interacting with each other and with the human host. It is a daunting task to develop predictive models for such a system, yet the potential rewards are high and would, for instance, enable targeted interventions to shift dysbiotic communities towards more healthy states. Two conditions need to be fulfilled for predictive models to be successful: first, the system has to be sufficiently well characterized to build the model and second, the dynamics should be generally deterministic. First successes in modeling the behavior of gut microbial communities give reason for cautious hope (Buffie, Bucci et al., 2015, Cremer, Arnoldini et al., 2017, Muñoz-Tamayo, Giger-Reverdin et al., 2016, Stein, Bucci et al., 2013. Most of these studies took a top down approach, where the change in composition of an entire community in vivo is modeled. For instance, Cremer and colleagues predicted the ratio of Firmicutes and Bacteroidetes in fecal samples as a function of estimated water content and nutrient influx using a diffusion model (Cremer et al., 2017). Others have fitted population models to time series of taxon (mostly genus) abundances obtained from 16S rRNA gene sequencing. For instance, one study fitted a variant of the generalized Lotka-Volterra (gLV) model to a cecal gut time series of mice exposed to a pathogen (Clostridium difficile), an antibiotic or both, thereby inferring the interactions between different genera (Stein et al., 2013). The time series predicted with the parameterized gLV model agreed well with observations. The same approach was also used to predict species that inhibit C. difficile growth in murine and human microbiota, one of which significantly lowered mortality when transferred to mice before infection with C. difficile (Buffie et al., 2015).
Despite these successes, the gLV model and its variants have several drawbacks that limit their widespread application. Importantly, they assume that community dynamics can be predicted from pair-wise interactions and that the interaction mechanism can be ignored. These assumptions have recently been challenged both experimentally and computationally: Friedman and colleagues experimentally quantified the accuracy reached when predicting the behavior of more complex soil communities from species pairs (Friedman, Higgins et al., 2017), whereas Momeni and colleagues systematically compared gLV models of metabolite-mediated species interactions to their mechanistic counterparts in silico (Momeni, Xie et al., 2017). While the authors in the former case concluded that the behavior of larger communities could to a considerable extent be predicted from that of smaller ones, the latter showed that the (extended) gLV model cannot describe accurately several common types of interaction mechanisms.
An alternative to the gLV model and its variants are mechanistic models, which account for nutrient-mediated interactions by explicitly describing the dynamics of the produced and consumed compounds (see (Momeni et al., 2017) and references therein). They thus require more system knowledge than generic gLV and related models do. However, most members of the gut community have not been thoroughly characterized, and little is known about their responses to different nutrients, pH values and interaction partners -even for those that have been studied more closely. It is challenging to obtain this type of biological knowledge and to resolve interaction mechanisms in vivo. However, in vitro studies not only allow acquiring detailed knowledge of the microorganisms' pH and nutrient preferences, but also of their behavior in the presence of other microorganisms. In vitro studies of human gut microorganisms have a long tradition and have been carried out in several manners. Classical mono-and coculture studies in batch and chemostat fermentors have explored nutrient preferences and interaction mechanisms (Falony, Calmeyn et al., 2009a, Falony, Vlachou et al., 2006, Moens, Verce et al., 2017, Moens, Weckx et al., 2016, Rivière, Selak et al., 2016. Artificial gut systems, such as the TNO In Vitro Model of the Colon (TIM-2) (Venema, 2015) and the Simulation of the Human Intestinal Microbial Ecosystem (SHIME) (T. Van de Wiele, P. Van den Abbeele et al., 2015), seek to reproduce as closely as possible the conditions of the human gastrointestinal tract in a well-controlled manner. They are composed of several compartments connected with peristaltic pumps that vary in their pH to mimic the pH gradient in the human intestines. In addition, the SHIME system can be extended with a module dedicated to the investigation of human-microbiome interactions (HMI) (Marzorati, Vanhoecke et al., 2014). The gut community has also been studied at smaller scales, in minibioreactor arrays (Auchtung, Robinson et al., 2015) and with gut-on-chip microfluidic devices (Kim, Huh et al., 2012, Shah, Fritz et al., 2016. However, in most cases, gut simulators are inoculated with fecal material. In the range from top-down to bottom-up approaches of gut microbial community dynamics, these can be considered as intermediate cases, where the host is eliminated, but the community is not further simplified. The goal of these studies is usually to quantify the behavior of the entire community under different conditions. In the case of HuMiX and SHIME's HMI module, the interaction of particular gut microorganisms with epithelial cells is targeted (Marzorati et al., 2014, Shah et al., 2016. Since the exact composition of fecal material (which also includes bacteriophages and fungi) is difficult to resolve, it is hard to track each member in such a community. Although the in vitro dynamics of colon (Kettle, Louis et al., 2015) and rumen (Muñoz-Tamayo et al., 2016) communities has been described with mechanistic models previously, these models did not account for the behavior at species level, but grouped species with similar metabolic activities into guilds. Whereas it is of interest to model guild dynamics, the resolution of guild-level models may be insufficient to understand microbial community dynamics in the gut. Species in the same guild do not necessarily respond in the same manner to altered environments and perturbations. Guild definitions are arbitrary to an extent, and gut bacteria with flexible metabolic strategies may change their guild membership. In addition, the concepts of tipping elements (Lahti, Salojãrvi et al., 2014) and strongly interacting species (Gibson, Bashan et al., 2016) suggest that particular species can have a disproportionate impact on gut community dynamics.
In our opinion, experiments with defined communities of known composition, grown in well-controlled conditions, are crucial to learn more about the interactions of gut species and how these shape community dynamics. Wellcontrolled in vitro experiments are also necessary for the development and validation of predictive models of gut microbial communities. However, only few in vitro experiments with defined gut communities have been reported to date (Newton, Macfarlane et al., 2013, Trosvik, Rudi et al., 2008, Trosvik, Rudi et al., 2010 and none have to our knowledge employed mechanistic models to predict community dynamics. The objective of the present study was to establish a defined gut community composed of representative human gut strains, to study their interactions under well-controlled conditions in vitro and to validate a quantitative mechanistic model by predicting community behavior in a tri-culture with parameters from mono-culture data. Whereas mechanistic models have been tested in this manner before for a cystic fibrosis community (Schmidt, Riedele et al., 2011), such an approach has not yet been applied to a synthetic gut community.
To reach our objective, we established three gut bacterial strains, described below, in an optimized mMCB medium (Moens et al., 2016) in 2-L laboratory fermentors run in batch mode. To measure kinetic parameters, we carried out mono-cultures and then cultivated all three strains together in a tri-culture. In addition, we quantified the dynamics of each of the three possible pair-wise combinations and performed all experiments in at least two biological replicates. For each cultivation experiment, we took samples at 15 to 18 time points over two days. For all samples, we quantified optical density, counted cells with flow cytometry, measured the concentration of fructose and short-chain fatty acids (formate, acetate, propionate, butyrate and lactate) and quantified bacterial abundance with TaqMan qPCR using species-specific primers and probes. In addition, we quantified the concentration of CO2 and H2 in the headspace during fermentation in a semi-continuous fashion. Finally, we sequenced the total RNA in selected samples. Figure 1 summarizes our approach. To our knowledge, this is the first study that investigates a synthetic gut community with a combination of mono-and co-cultures, mechanistic modeling and gene expression analysis.
Our synthetic community was composed of Faecalibacterium prausnitzii A2-165 (Duncan, Hold et al., 2002b) (FP), Roseburia intestinalis L1-82 (Duncan, Hold et al., 2002a) (RI) and Blautia hydrogenotrophica S5a33 (Bernalier, Willems et al., 1996) (BH), all three of which were isolated from human feces. We carefully selected our community members according to several criteria. First, we targeted abundant and typical members of the human gut microbiome. According to the Flemish Gut Flora Project (FGFP), for which fecal samples of 1,106 healthy Flemish adults were sequenced, Faecalibacterium is the topabundant genus and a Faecalibacterium OTU the top-abundant OTU in the healthy gut (Falony*, Joossens* et al., 2016). Roseburia and Blautia are both among the top 10 most abundant gut genera. In addition, Faecalibacterium and Roseburia OTUs were present in all sequenced samples, whereas Blautia OTUs were present in 99.9% of the samples. The selected Faecalibacterium and Roseburia strains were also among the strains with more than 1% abundance in several of the whole-genome-shotgun sequenced stool samples from the Human Microbiome Project (Kraal, Abubucker et al., 2014). Thus, our community is composed of representative human gut strains. Second, we specifically targeted butyrate producers (RI and FP). Butyrate is an important energy source for gut epithelial cells and is beneficial for gut health (Geirnaert, Calatayud et al., 2016, Rivière et al., 2016. Butyrate producers are often found to be enriched in healthy as compared to dysbiotic gut microbiota (Antharam, Li et al., 2013, Rivera-Chávez, Zhang et al., 2016. Thus, high butyrate production will likely be a quality criterion for bacterial cocktails designed for therapeutic purposes. Our synthetic community allows studying more closely the factors that influence butyrate production. Third, we designed our community with a large number of potential interactions. As Figure 2 illustrates, our community contains multiple cross-feeding and competitive interactions. For instance, FP needs acetate, which is produced by BH. Vice versa, BH grows on the formate generated by FP. Both cross-feeding relationships together constitute a mutualistic interaction. In parallel, both species compete for the energy source, fructose. The same relationships also occur between RI and BH (which can exchange hydrogen and CO2 in addition to formate), while RI and FP compete for fructose and acetate. This system also constitutes a rare example of two species pairs that simultaneously compete and mutually cross-feed. Finally, we only considered human gut species for which at least a draft genome was available, to ease primer design and the interpretation of RNA-seq data.

Blautia hydrogenotrophica is metabolically versatile
We first confirmed the cross-feeding interactions postulated for BH with smallvolume screening experiments, where the pH was not kept constant and the atmosphere contained 10% CO2 and 10% hydrogen gas. We found that under these conditions, BH was able to grow heterotrophically on formate and autotrophically as described previously on CO2 and hydrogen gas (Bernalier et al., 1996), presumably in both cases via the Wood-Ljungdahl pathway, of which all genes were located in BH's genome (Rey, Faith et al., 2010). We also found in agreement with (Rey et al., 2010) that BH grew on fructose, oligofructose and glucose, albeit more slowly than on formate. In agreement with (Bernalier et al., 1996), we detected lactate in addition to acetate for these substrates. We thereby confirmed the potential for fructose competition between BH on the one hand and RI or FP on the other in our medium. BH also consumed small concentrations of galactose, but did not consume fucose, inulin or lactate. For BH grown on fructose, we confirmed slow growth and lactate production in addition to CO2, H2 and acetate in the fermenter. Notably, when growing BH on formate in the fermenter, carbon dioxide and hydrogen gas were produced besides acetate, but no lactate.
In conclusion, our work showed that Blautia hydrogenotrophica is a surprisingly versatile member of the human gut ecosystem.

Mono-culture dynamics does not follow standard kinetics
We employed pH-controlled monocultures to characterize the properties and growth kinetics of the individual strains in our model. Table 1 provides an overview for all fermentation experiments carried out, whereas Supplementary Table 1 gives additional information for each experiment. RI consumed fructose and produced butyrate, CO2 and hydrogen gas, as described previously (Falony, Verschaeren et al., 2009c) as well as small amounts of lactate and formate ( Figure 3A). Interestingly, there was no net consumption of acetate when more fructose than acetate was provided. While net acetate consumption has been found to correlate negatively with hydrogen production (Falony et al., 2009c), we saw here that it also depends on the ratio of initial fructose and acetate. When given in equal concentrations, RI partially consumed acetate. In all further experiments, when adding acetate, it was added at the same concentration as fructose. FP in mono-culture produced formate, less CO2 and butyrate than RI and no hydrogen gas, but did not entirely consume fructose ( Figure 3B). After having excluded exposure to oxygen (by adding oxygen gas via sterile water), redox potential (by continuously adding the oxidizing agent potassium ferrocyanide trihydrate), pH (lowered to 5.8), a threshold requirement for fructose (halving the fructose concentration did not stop its consumption) or end-product inhibition (by adding initial butyrate) as explanations, we found that doubling the concentration of yeast extract lowered fructose concentrations. Since adding fresh but autoclaved medium during the fermentation did not lower fructose concentrations, we assume that FP was growth-limited by (a) heat-labile cofactor(s) present in the yeast extract. A recent flux balance analysis with a manually curated metabolic reconstruction suggests that FP requires several amino acids (L-alanine, L-cysteine, L-methionine, L-serine and L-tryptophan) and the co-factors biotin (vitamin B7), cobalamin (vitamin B12), folic acid (vitamin B9), hemin, nicotinic acid, pantothenic acid and riboflavin (vitamin B2) for growth (Heinken, Khan et al., 2014). With the exceptions of cobalamin and externally supplied hemin, these nutrients should be present in yeast extract according to the metabolic reconstruction of Saccharomyces cerevisiae iMM904 (Mo, Palsson et al., 2009) and the amino acids should furthermore be present in other medium components (bacteriological peptone, soy peptone and tryptone). Heinken and colleagues also predicted that FP can grow in the presence of oxygen, which is in agreement with our observation that addition of low concentrations of oxygen does not alter FP's growth curve (Heinken et al., 2014). BH produced acetate, hydrogen, CO2 and small concentrations of lactate, while consuming formate almost entirely ( Figure 3C). It also consumed fructose, but did not deplete it. Whereas the carbon recovery for the other mono-cultures was close to 100%, it only reached 60% for BH in mono-culture on formate and fructose. These unexpected behaviors defy simple kinetic models and necessitate further adjustment of the equations.

Prediction accuracy of model parameterized on mono-cultures is speciesdependent
We designed a model that described the dynamics of each species and of key compounds (including fructose, formate, acetate, butyrate, hydrogen gas and CO2) with ordinary differential equations assuming Monod kinetics (see Material and Methods). The model differentiated between substrates required for growth and co-substrates such as acetate that enhanced growth but were not required. It also took species-specific differences in lag phases into account. Since we observed that FP did not deplete fructose, presumably because of a lack of cofactors, we introduced a dependency on an undefined metabolite referred to as "unknown compound". We parameterized this model on selected mono-culture experiments and then predicted mono-culture dynamics (Figure 4 A-C, Supplementary Figure 1). The model reached high prediction accuracy for FP and RI, but did not describe well the experimental data for BH (see Table 1). More precisely, the model showed that BH did not consume formate and fructose as quickly as expected if its growth would follow Monod kinetics. We confirmed culture homogeneity by analyzing 16S rRNA gene sequencing data of the last sample (Supplementary Figure 2). A yeast contaminant (S. cerevisiae S288c) detected in RNA-seq data of BH mono-culture samples (Supplementary Figure 3) does not explain the incongruence between growth and energy source consumption, since i) no contamination was observed on plates inoculated with bioreactor samples and incubated in anaerobic and aerobic conditions, ii) S. cerevisiae would consume fructose and iii) no ethanol production was measured. We also found only small concentrations of potential peptide degradation products (isobutyric acid and isovaleric acid). We therefore assume that BH in mono-culture initially grew on undefined medium components and only later switched to formate and fructose, but the time resolution was too low to model this potentially biphasic growth. We also compared model performance for RI with and without product inhibition by hydrogen gas. Since we found no differences in model performance, we removed an initial hydrogen gas inhibition term.

Blautia consumes formate produced by Roseburia and Faecalibacterium
When growing FP and BH together, we observed that fructose was entirely depleted and acetate, butyrate, hydrogen gas, CO2 and small concentrations of lactate were produced ( Figure 3F). Interestingly, there was an initial production of formate, which was then consumed, confirming that formate was cross-fed from FP to BH. Formate consumption was also observed without initial acetate ( Figure 3H). In the bi-culture of RI and BH, CO2, hydrogen, butyrate and small concentrations of lactate were produced, whereas fructose and a small amount of acetate were consumed ( Figure 3D). The same fermentation products were also obtained in the absence of initial acetate ( Figure 3G). In contrast to RI in mono-culture, no formate was detected, suggesting that it was entirely cross-fed to BH. It was unclear whether CO2 and hydrogen gas produced by RI reached sufficiently high concentrations to be cross-fed to BH. Finally, when RI and FP were co-cultivated, fructose and acetate were consumed and butyrate, formate, hydrogen gas and CO2 were produced ( Figure 3E). The finding that formate reached lower concentrations than in FP mono-culture already hints at a negative effect of RI on FP.

Comparison of mono-and co-culture data suggests ecological interactions
Since Gause's early work on competition between yeast and Paramecium species (Gause, 1932, Gause, 1934, growth rates in mono-and bi-culture experiments have been compared to determine ecological interactions (e.g. (de Vos, Zagorski et al., 2017, Freilich, Zarecki et al., 2011, Wang, Wei et al., 2017). The rationale is that growth rates in bi-culture should increase for mutualistic organisms as compared to mono-culture growth rates, whereas bi-culture growth rates should decrease for competitors.
When comparing maximal abundances, cross-feeding and competitive interactions were already apparent. Both FP and BH reached significantly higher maximal bacterial counts in FP-BH bi-cultures and in tri-cultures with FP dominance ( Figure 3F, 3H and 3J) than they did in mono-culture ( Figure 3B and 3C), suggesting a mutualistic relationship (unpaired two-sided Wilcoxon FP shift: 0.4, FP 95% confidence interval: 0.12 and 0.55, FP p-value: 0.03; BH shift: 0.5, BH 95% confidence interval: 0.33 and 0.69, BH p-value: 0.017). The maximal cell number of FP tended to be lower when competing with RI ( Figure 3E) than when grown alone (unpaired two-sided Wilcoxon shift: 0.47, 95% confidence interval: -0.03 and 1.42, p-value: 0.11). Interestingly, there was no difference in maximal bacterial counts for RI alone versus grown with FP in bi-cultures and tri-cultures with RI dominance (unpaired two-sided Wilcoxon shift: 0.07, 95% confidence interval: -0.39 and 0.31, p-value: 0.69), so that formally, their relationship could be described as amensalism (one organism is affected negatively whereas the other is not affected). Finally, according to maximal bacterial counts, BH benefited more from the presence of FP than from RI (unpaired two-sided Wilcoxon shift: 0.29, 95% confidence interval: 0.06 and 0.93, p-value: 0.008).

Model needs bi-culture data to accurately predict tri-culture dynamics
When growing all three gut bacterial strains together, fructose was consumed and butyrate, acetate, CO2, hydrogen gas as well as lactate were produced. Formate was produced initially, but quickly consumed ( Figure 3I and 3J). We performed the tri-culture six times with varying species proportions in the inoculum and found that in all tri-cultures, BH was always co-dominant, together with either RI or FP. In two out of the six cases, RI won over FP as the codominant partner of BH, whereas in the four other cases, FP won. The result mattered for the final butyrate concentrations, which averaged 37.5 mM when RI won and 23.5 mM when FP won. We attempted to describe tri-culture dynamics with the model parameterized on mono-cultures, but failed to obtain a good fit (see Table 1 and Supplementary  Figures 4 and 5). After a series of tests, we concluded that incorporating coculture data was necessary to describe tri-culture dynamics. We finally selected two FP mono-cultures and the RI-BH and FP-BH bi-cultures with initial acetate to parameterize our model. As a validation, we predicted the behavior of RI-BH and FP-BH bi-cultures in the absence of initial acetate, which resulted in a good fit ( Figure 5G and 5H, Supplementary Figure 7). The model parameterized on mono-and bi-cultures fitted the tri-culture data better than the model parameterized on mono-cultures only (Table 1, Figure 4I and 4J, Figure 5I and 5J, Supplementary Figures 5 and 8).
When inspecting the differences between the two parameterizations, we found that the model parameterized on mono-cultures predicted lower abundances for all three species in bi-and tri-cultures than they actually reached ( Figure 4D-J, Supplementary Figures 4 and 5). Vice versa, the model parameterized on monoand bi-cultures predicted too high abundances for RI and BH in mono-culture ( Figure 5A and 5C, Supplementary Figure 6; the FP mono-culture was included in the parameterization). According to the higher maximal cell counts in tri-culture predicted with the second than the first parameterization (unpaired Wilcoxon: 0.002), BH did significantly better in tri-culture than expected based on its mono-culture growth.
The fact that the model describes well both mono-and tri-culture dynamics is a sign of emergent behavior in the presence of interaction partners. When looking at the parameters inferred from mono-and bi-cultures (given in Supplementary  Table 2), BH's consumption rates for formate and fructose and RI's fructose consumption rate were lower compared to their values obtained from monoculture parameterization, whereas their maximal growth rates were not much affected (BH) or increased (RI). Thus, according to this analysis, less of the energy source is needed in the presence of an interaction partner than in monoculture.

Initial abundance and lag phase predict species dominance in tri-culture
Next, we tested whether dominance in tri-culture could be predicted from lag phase and initial abundance. Towards this aim, we computed the FP/RI ratio in simulations with varying lag phase and initial abundance. Experimental observations of dominance agreed well with model predictions ( Figure 6A and 6B). Our systematic investigation also showed that there was a non-linear relationship between initial RI abundance and FP dominance ( Figure 6C-F). Thus, even when knowing initial abundances, lag phases and species interactions, it is hard to predict the winner (and thereby resulting butyrate concentration) intuitively without a model in hand. The final abundances of the three strains in simulations are also non-linearly dependent on other parameters, including BH's growth rate, its fructose consumption rate and its fructose half-saturation constant (Supplementary Figure 9). These results underline that besides kinetic parameters, initial conditions and lag phase can determine strain abundances in co-culture in a non-linear way.

Altered gene expression in response to interaction partners provides first insights into emergent behavior
To further investigate the emergent behavior, we sequenced RNA for the three mono-cultures and the tri-culture where FP co-dominated for three time points and two biological replicates and assessed significantly differential gene expression across all samples in mono-versus tri-cultures for all three species (Supplementary Table 3). In total, 6.7%, 7.9% and 1.6% of RI's, FP's and BH's proteins respectively were significantly differentially expressed (protein numbers taken from UniProt (The UniProt Consortium, 2017)). Interestingly, in tri-culture, FP down-regulated a series of enzymes needed for vitamin B12 coenzyme biosynthesis. Since cobalamin (vitamin B12) was one of the co-factors suspected to limit FP growth in mono-culture, this finding may mean that FP benefited from greater co-factor availability in tri-culture. In tri-culture, FP also up-regulated enzymes involved in amino acid and oligopeptide transport and amino acid and protein biosynthesis. BH likewise up-regulated amino acid biosynthesis in tri-culture. For RI, which reached lower abundances in the selected tri-cultures than in mono-culture, the transcription response was mixed: some amino acid biosynthesis enzymes were down-regulated, others upregulated (including enzymes involved in ornithine biosynthesis). However, expression of ribosomal proteins was lower than in RI mono-culture, in agreement with its long lag phase in the selected tri-cultures. In summary, the analysis of differential gene expression uncovered a number of metabolic changes in the presence of interaction partners, thus supporting further the altered behavior detected through modeling.

Discussion
Here, we investigated the dynamics of a well-defined small but representative gut microbial community in vitro. We found that BH was metabolically versatile and grew as fast as primary fermenters such as RI. We demonstrated experimentally that formate was cross-fed between BH on the one hand and FP and RI on the other and confirmed mutualistic as well as competitive interactions between these three bacterial strains. When growing on formate, we identified BH to be a net producer of both hydrogen and carbon dioxide, in contrast with its traditionally assumed role in the gut ecosystem. While formate is rarely highlighted as a key intermediate in gut cross-feeding interactions, it has been reported to be an end-product of primary polysaccharide degradation by both Bifidobacterium and Lactobacillus spp. (Falony, Lazidou et al., 2009b, Moens et al., 2017. Hence, our results invite a re-evaluation of the ecological niche of BH in relation with microbial formate production potential. The model, which encodes our knowledge of the system, is not only important for predictions, but also as a reference. We gained insights from its agreements as well as from its disagreements with our observations. For instance, we assumed initially that RI would be inhibited by the hydrogen gas it generated. However, a hydrogen gas inhibition term was not required to accurately describe RI behavior in mono-culture, which implied that hydrogen gas inhibition did not affect RI growth at the concentrations reached in our experiments. We also needed the model to ascertain that changes from mono-to bi-or tri-cultures were not just due to variations in the inoculum composition or the lag phase, but that there was a true change in the dynamics that the model parameterized on mono-cultures alone could not capture. We confirmed this emergent behavior with RNA-seq, which revealed significantly different gene expression in triculture as compared to mono-culture, especially for FP and BH. The downregulation of FP's vitamin B12 coenzyme biosynthesis pathway in tri-culture is of particular interest, since it suggests that dependency on co-factors changes with interaction partners. It has been posited that the majority of gut microorganisms in need of B12 precursors is unable to synthesize them (Degnan, Taga et al., 2014). If this need is altered by the presence of interaction partners, it cannot be exploited as easily for selective manipulation as suggested in (Degnan et al., 2014).
Although species-level kinetic models parameterized on mono-cultures may in some cases describe bi-culture dynamics correctly (Van Wey, Cookson et al., 2014), our example shows that this is not a general property. This means that models of multi-species communities will probably have to take the internal metabolism and its response to interaction partners into account. Gut bacteria such as BH have flexible metabolic strategies that they employ according to circumstances. Emergent metabolism in the presence of interaction partners has been described in theoretical work before (Chiu, Levy et al., 2014), but has been rarely investigated experimentally ((Aharonovich & Sher, 2016)). Constraintbased modeling approaches, which take the entire metabolism into account (Orth, Thiele et al., 2010), require high-quality metabolic reconstructions for each community member, which take months of curation effort to obtain . Thus, scaling species-level quantitative models to larger communities will be a formidable challenge.
Mono-and bi-cultures are increasingly carried out in batch in a high-throughput fashion to determine ecological interactions and to quantify their strengths (de Vos et al., 2017, Sher, Thompson et al., 2011. Such systematic quantification is an important step forward, but there are challenges to tackle. Our work showed that dominance in batch may sensitively (i.e. non-linearly) depend on initial conditions such as the lag phase and the initial abundance, both of which are hard to control experimentally. Thus, a growth experiment performed in biological replicates but with the same inoculum may identify one species as the winner and another as the loser. However, a replicate with a slightly different inoculum composition may reach the opposite conclusion. Such a dependency on the initial conditions (albeit with larger abundance differences) has also been reported in several competition experiments for Streptomyces species (Wright & Vetsigian, 2016) and may thus be a common case. To ascertain that bacteria change their behavior in response to an interaction partner, RNA-seq can be carried out on mono-and bi-cultures (Aharonovich & Sher, 2016, Plichta, Juncker et al., 2016. Here, we showed that a model can also reveal emergent behavior by its failure to describe co-culture dynamics when parameterized on monocultures only. It is an open question how to scale up such approaches to achieve the high-throughput needed for systematic measurements of interaction strengths. This is to the best of our knowledge the first in vitro investigation of a defined bacterial community that combines mutualism with competition in two cases. In the case of FP and BH, mutualism appears to supersede competition, leading to increased maximal bacterial numbers coupled with up-regulation of biosynthesis for both interaction partners. For RI and BH, we had no such clear experimental evidence, since RNA-seq was performed on tri-cultures dominated by FP, but the comparison of maximum bacterial numbers across mono-, bi-and tri-cultures suggests that RI and BH do not benefit as much from each other as FP and BH do. Since the model described RI-BH bi-culture dynamics well without taking CO2 and hydrogen gas cross-feeding into account, we assume that due to their low partial pressure gasses are less efficiently cross-fed to BH than formate, although both are likely metabolized via the same pathway (Wood-Ljungdahl). Thus, interactions that look similar on paper can play out differently, depending on the environment. In the two replicates of the RI-FP bi-culture, RI and FP both survived (albeit FP in far lower numbers) despite competing for the same substrate, presumably because in our experimental set-up, the time until nutrient depletion was too short for the competitive exclusion principle to apply.
Our tri-culture experiments also demonstrated the importance of initial conditions on fermentation end-products. According to our model, the initial abundance and lag phase determined whether butyrate reached high or low concentrations in the tri-culture fermentations. Since these are likely to be relevant parameters in the gut environment and difficult to control, cocktail communities will have to be designed such that they will carry out their function across a wide range of initial conditions. While our work highlighted a number of challenges to microbial community modeling, the model's ability to predict tri-culture dynamics from bi-cultures gives hope that with sufficient knowledge, we will ultimately be able to model more complex microbial communities.

Cultivation experiments in stationary bottles
Monoculture cultivation experiments for BH were performed in stationary glass bottles without controlling the pH (screening). The bottles contained 50 mL of heat-sterilized pH 6.8 mMCB medium, supplemented with either 50 mM of Dfructose (Merck), D-glucose (Merck), D-galactose (Merck), L-fucose (Merck), sodium formate (VWR International), sodium acetate trihydrate (Merck), DL lactic acid (VWR International), oligofructose (Raftilose P95; Beneo-Orafti NV, Tienen, Belgium; (Falony et al., 2009b)) or inulin (OraftiHP; Beneo-Orafti; (Falony et al., 2009b)) as the sole energy sources. Additional cultivation experiments were performed in medium devoid of any main energy source to test autotrophic growth. For the cultivation experiments in bottles, stock solutions of fructose, glucose, galactose, fucose, sodium formate, sodium acetate trihydrate, and lactic acid were initially made anaerobically through autoclaving at 210 kPa and 121 °C for 20 min. The solutions were subsequently filtersterilized and transferred into glass bottles, which were sealed with butyl rubber septa that were pierced with a Sterican needle (VWR International) connected with a Millex-GP filter (Merck) to assure sterile conditions. For the cultivation experiments with lactate, the pH was adjusted to 6.8 under anaerobic conditions, using sterile solutions of sodium hydroxide (Merck). Stock solutions of oligofructose and inulin were made sterile by membrane filtration. The inocula were prepared as follows. Cells of the strains under study were transferred from −80°C to test tubes containing 10 mL of RCM that were incubated anaerobically at 37 °C for 24 h. Subsequently, the strains were propagated for 12 h in glass bottles containing 50 mL of heat-sterilized pH 6.8 mMCB medium, supplemented with the energy source under study, always at a final concentration of 50 mM fructose equivalents. These pre-cultures were finally added to the glass bottles aseptically. During the inoculum build-up, the transferred volume was always 5% (vol/vol). All bottles were incubated anaerobically at 37°C in a modular atmosphere-controlled system (MG anaerobic work station; DonWithley Scientific Ltd., West Yorkshire, United Kingdom) that was continuously sparged with a mixture of 80% N2, 10% CO2, and 10% H2 (Air Liquide, Paris, France). Samples for further analyses were withdrawn after 0, 6, 12, 24, 48, and 100 h. All experiments were performed at least in duplicate.

Fermentation experiments
To prepare inocula, cells were transferred from −80 °C to test tubes containing 10 mL of RCM, and incubated at 37 °C for 24 h. Subsequently, the strains were propagated twice for 12 h in glass bottles containing 100 mL of mMCB medium (with acetate in the case of RI and FP, and with formate in the case of BH), supplemented with fructose. All incubations were performed anaerobically in a modular atmosphere-controlled system (MG anaerobic workstation) that was continuously sparged with a mixture of 80% N2, 10% CO2, and 10% H2 (Air Liquide). The inocula were finally added aseptically to the fermentors. During the inoculum build-up, the transferred volume was always 5% (vol/vol). Fermentations were carried out in 2-L Biostat B-DCU fermentors (Sartorius) containing 1.5 L of mMCB medium supplemented with the co-substrates (acetate and/or formate) if necessary and 50 mM of D-fructose as the energy source. Anaerobic conditions during fermentations were assured by continuously sparging the medium with N2 (PraxAir, Schoten, Belgium) at a flow rate of 70 mL min −1 . The fermentation temperature was kept constant at 37 °C. A constant pH of 6.8 was imposed and controlled automatically, using 1.5 M solutions of NaOH and H3PO4. To keep the medium homogeneous, a gentle stirring of 200 rpm was applied. Temperature, pH, and agitation speed were controlled online (MFCS/win 2.1 software, Sartorius). Fermentations were followed for 48 h, with samples taken at 10 min and 2 h, 3 h, 5 h, 6 h, 7 h, 9 h, 10 h, 11 h, 13 h, 14 h, 15 h, 17 h, 18 h, 24 h, 30 h and 48 h after inoculation. At selected time points (3 h, 9 h and 15 h after inoculation), subsamples were treated for RNA extraction by adding 5 vol of RNAlater® (xxx company). All mono and tri-culture fermentations were performed in triplicate. All bi-culture fermentations were performed in duplicate, except for the bi-culture fermentations using medium lacking acetate with FP and BH that were performed only once.

Quantification of bacterial abundance
During all experiments, the optical density at 600 nm (OD600) was measured against ultrapure water as blank with a VIS spectrophotometer (Genesys 20; Thermo Scientific, Waltham, MA, USA). Each measurement was performed in triplicate. Total bacterial abundance was also measured by flow cytometry, using an Accuri C6 flow cytometer (BD Biosciences, Erembodegem, Belgium), as described previously (Moens et al., 2016). All samples were diluted in filtersterilized water (Vittel, France) to obtain a concentration between 1.0 × 10 3 and 5.0 × 10 6 cells mL −1 . Flow cytometric analysis was performed by mixing 500 μL of sample with 5 μL of a 100× SYBR Green I solution (Sigma-Aldrich) and 5 μL of a 500 mM ethylenediaminetetra acetic acid (EDTA) solution (Sigma-Aldrich). Afterwards, samples were left in the dark at room temperature for 15 min. Flow cytometric counts were obtained using an Accuri C6 flow cytometer (BD Biosciences), equipped with a 50 mW solid state laser (488 nm). Green fluorescence was measured in the FL1 channel (530 ± 15 nm) and all data were processed with the Cflow Plus software (Accuri). Gating was performed to distinguish signals from noise. All data were collected as a FL1/SSC density plot with a primary threshold of 10,000 on the FL1 channel. Measurements were performed in triplicate. qPCR assays with species-specific TaqMan primers and probes were performed to quantify the abundance of each species separately. For this, 2 mL of fermentation sample was centrifuged at 20,570 × g for 20 min. Cells were washed in 2 mL of physiological solution (NaCl, 8.5 g L −1 ) and centrifuged again at 20,570 × g for 20 min to obtain washed cell pellets. Subsequently, these cell pellets were resuspended in 2 mL of physiological solution and diluted 20 times for DNA extraction. Direct DNA extractions by alkaline thermal lysis were performed based on (Girish, Haunshi et al., 2013, Rudbeck & Dissing, 1998, modified as follows: 100 μL of the sample was mixed with 100 μL of 0.2 M NaOH in a sterile micro centrifuge tube. The mixture was vortexed and heated at 90°C for 10 min, after which eight volumes (1600 μL) of 0.04 M Tris HCl pH 7.5 (Thermo Fisher Scientific) was added for pH neutralization. 4 μL of final mixture was used for qPCR. The extracted genomic DNA was stored at -20°C until qPCR amplification. Calibration curves were obtained by initially growing all strains in RCM for 24 h, and two-fold propagation in medium for 12 h, as described above. From each of these grown cultures, separate fourfold decimal and nine-fold binary dilution series were prepared. The generation of cell pellets, direct extraction of DNA, and subsequent quantification of cell concentrations by flow cytometry were performed as described above, with the exception that, prior to DNA extraction, samples for calibration were diluted less than fermentation samples, to accommodate a wider qPCR quantification range. Primers and oligoprobes (Supplementary Table S4) were manually designed using the online Primer3Plus software (Untergasser, Nijveen et al., 2007) and the genome sequences of the strains. Melting temperatures and presence of hairpins, self-dimers, and pair-dimers were double-checked using the online OligoCalc software (Kibbe, 2007). Secondary structures of the generated amplicons was investigated using the online Mfold program (Zuker, 2003). Primers and probes were synthesized by Thermo Fisher Scientific. Strain specificity of primers and probes was confirmed in silico by Primer-BLAST (Ye, Coulouris et al., 2012) and in vitro by PCR and qPCR analysis on genomic DNA of the strains (Supplementary Table S1). qPCR assays were carried out using a 7500 FAST Real-Time PCR system (Applied Biosystems, Carlsbad, CA, USA), equipped with 96-well plates. Each qPCR assay mixture of 20 μL contained 10.0 μL of TaqMan® Fast Universal PCR Master Mix (2X), no AmpErase® UNG (Thermo Fisher Scientific), 2.0 μL of each primer (3.0 μM), 2.0 μL of the TaqMan probe (1.5 μM), and 4.0 μL of extracted genomic DNA solution or sterile nuclease-free water (Thermo Fisher Scientific). The qPCR amplification program consisted of an initial denaturation step at 95°C for 20 s, followed by 45 two-step cycles at 95°C for 3 s and at 60°C for 30 s. In each run, negative (sterile nuclease-free water without genomic DNA) and positive controls (with extracted genomic DNA from the relevant strains) were used. The cycle threshold (Ct) values were determined using the automatically determined thresholds from the 7500 software v2.0.6 (Applied Biosystems). Finally, during a re-analysis of all qPCR runs, Ct values were normalized using an inter-plate calibrator to account for differences among qPCR runs. The above-described generation of cell pellets, direct extraction of DNA, and qPCR assays were performed in triplicate. Contamination was checked by aerobic and anaerobic plating on RCM agar and 16S rRNA gene amplicon sequencing of end point fermentation samples (48 h). Sequencing was performed as described previously (D'hoe, Conterno et al., 2018).

Metabolite profiling
Concentrations of fructose, as well as concentrations of formate, acetate, butyrate, lactate, and ethanol were determined through high-performance liquid chromatography (HPLC) with refractive index (RI) detection, using a Waters chromatograph (Waters, Milford, MA, USA) equipped with an ICSep ICE ORH-801 column (Transgenomic North America, Omaha, NE, USA), and applying external standards, as described previously (Falony et al., 2009b). Briefly, the mobile phase consisted of 5 mM H2SO4 at a flow rate of 0.4 mL min -1 . The column temperature was kept constant at 35°C. Sample preparation involved a first centrifugation (4,618 x g for 20 min at 10°C) for removal of cells and debris, followed by the addition of an equal volume of 20% (mass/vol) trichloroacetic acid for protein removal. For determining oligofructose and inulin consumption, samples were incubated at room temperature for 24 h to assure complete hydrolysis of polysaccharides. Subsequently the samples were centrifuged (21,912 x g, 20 min, 4°C) and filtered (pore size of 0.2 μm; Uniflo 13 Filter Unit; GE Healthcare, Little Chalfont, United Kingdom), prior to injection (30 μL) into the column. Samples were analyzed in triplicate. For an additional screening experiment with BH grown in the presence of 350 mM formate, the concentrations of ethanol, acetoin, acetic acid, propionic acid, butyric acid, isobutyric acid and isovaleric acid produced were determined via gas chromatography with flame ionization detection (GC-FID), using a FocusGC chromatograph (Interscience, Breda, The Netherlands) equipped with a Stabilwax-DA column (Restek, Bellefonte, PA, USA), and applying external standards, as described previously (Moens, Lefeber et al., 2014). The samples were analyzed in triplicate.

Model definition
We modeled change of species abundances over time with the following three ordinary differential equations: where X denotes species abundance, S metabolite concentration and Q a lag phase parameter. The growth rates are then defined as nonlinear growth functions as in (Grivet, 2001, Smith & Waltman, 1995, which assume Monod kinetics (Monod, 1950): where K is the Monod (half-saturation) constant, μ is the maximal specific growth rate and ω a weight parameter. Nutrient dependency can be either obligatory (growth without nutrient is not possible) or facultative (growth without nutrient is possible). For instance, the fructose uptake is multiplied with RI's maximal growth rate, whereas its acetate uptake is modeled with an additive term. Therefore, in the absence of fructose, RI's growth rate is zero, but not in the absence of acetate. The weight parameter adjusts how strongly a facultative substrate contributes to the overall growth rate. The unknown compound models the dependency of FP on an undetermined co-factor.
The lag phase function is defined as in (Baranyi & Roberts, 1994): , where i stands for RI, FP or BH.
The Qi variables follow exponential growth: Thus, the larger the initial value of Qi, the shorter the lag phase.
The changes of metabolite concentrations are then modeled as follows: The α and ν parameters are production and consumption rates, respectively. Species abundance is measured in 10 8 bacterial counts/mL, metabolite concentration in mM, the unit of μ is 1/h, the unit of K is mM, the unit of α and of ν is mM/(10 8 bacterial counts/mL) and ω is dimensionless. The model assumes that death rates are negligible. CO2 and hydrogen consumption by BH is not included in the final version of the model. We tried to account for CO2 consumption with a multiplicative term in BH's growth rate. However, this did not improve the model fit. Since the model without CO2 describes RI-BH bi-culture dynamics well, we assume that BH grows mostly heterotrophically on fructose and on the formate produced by RI and that the hydrogen gas and CO2 produced by RI did not reach sufficient concentrations in the head space to allow autotrophic growth as observed during the screening, where the atmosphere contained 10% of hydrogen gas and 10% of CO2. The model definition is available as Source Code File in python (Model definition).

Model parameterization
We parameterized our model on mono-cultures alone (parameterization 1) and on mono-and bi-cultures (parameterization 2). The goodness of fit of both parameterizations are summarized in Table 1. The model was fitted using the function fmin() from the scipy Python package (Jones, Oliphant et al., 2001), to minimize the normalized root mean square error. An initial estimate of the parameters was obtained by manually fitting the data iteratively. The initial concentration of the unknown compound was set to 30 mM. Samples after the end of the log phase, when the bacterial counts started to decline, were omitted from the fitting. Parameterization 2 consisted of several steps, as fitting all parameters at once did not lead to convergence, due to the nonlinear growth rates. For this, FP was first fitted on two FP mono-cultures and cross-validated on the third one. The consumption parameters of BH were obtained from FP-BH bi-cultures with initial acetate, afterwards the maximal specific growth rates and half-saturation constants were obtained from the same bi-cultures. RI was fitted on a RI-BH bi-culture with acetate. Lag phases were calculated as the time to reach Γ i (Q i ) = 0.5 : (0) was estimated by visual inspection of the log plots. Model parameters obtained and maximal abundances predicted with both parameterizations as well as estimated lag phases are provided in Supplementary Table 2. Data and model fits were plotted with Python's matplotlib (Hunter, 2007).

RNA extraction and sequencing
Total RNA was extracted from RNAlater®-treated samples using the phenol-free total RNA purification kit coupled with DNase I treatment (VWR International) according to the manufacturer's protocol for Gram-positive bacteria. A secondary DNAse digestion was performed using the Ambion® TURBO DNA-free™ DNase Treatment and Removal Reagents Kit (Thermo Fisher Scientific), after which the samples were purified using the RNA Clean & Concentrator™-25 kit (Zymo Research, Irvine, CA, USA) according to the manufacturer's instructions.
The eluted RNA was stored at -80°C. The absence of DNA contamination was evaluated using PCR (35 or 40 cycles) and gel electrophoresis. The concentrations of the samples were determined with a Nanodrop, and with a Qubit 2.0 fluorometer using the Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific). RNA integrity, expressed as the RNA integrity number (RIN), and yield were determined using RNA Nano/Pico 6000 LabChips (Agilent Technologies, Santa Clara, CA, USA) that were run in an Agilent 2100 Bioanalyzer (Agilent Technologies). Whereas most of the RINs were above 7, RINs of three BH mono-culture samples at 3 h, 9 h and 15 h were around 2.6 and in four cases (BH mono-culture at 15 h, FP mono-culture at 9 h, and for both tri-culture replicates at 3 h) the RINs could not be determined. However, by pooling over three extraction rounds, sufficient RNA for sequencing (minimum of 536 ng and median of 2800 ng) was obtained for all samples. Library preparation encompassed the use of Ribozero rRNA removal for Grampositive Bacteria and the Illumina TruSeq stranded mRNA Library preparation kit (IIlumina, San Diego, CA, USA). Library preparation was performed without the mRNA purification step, according to the manufacturer's instructions. The enriched libraries were sequenced on an Illumina NextSeq 500 instrument (paired-end, 2×76 bp reads, Mid output kit, Illumina). From the Illumina platform, paired-end reads in FASTQ format (CASAVA 1.8, Phred + 33) were obtained and separated into distinct files for each single-end read and for each sample.

RNA-seq analysis
The analysis of the raw sequencing reads was performed as follows: reads were trimmed using Trimmomatic (Bolger, Lohse et al., 2014) with the following parameters "CROP:74 HEADCROP:10 SLIDINGWINDOW:4:15 MINLEN:51", to remove initial and last bases which had biases in their nucleotide content as reported by FastQC (Andrews, 2010), to remove stretches of low-quality bases and to keep reads with at least 51 bases after trimming. FastQC was re-run on the trimmed data to ensure the previous biases were corrected. SortMeRNA (Kopylova, Noé et al., 2012) was used with default parameters and included databases to remove rRNA reads.
With the remaining non-rRNA reads, we ran MetaPhlAn2 (Truong, Franzosa et al., 2015) with default parameters and database, and mash screen (Ondov & Philippy, 2017, Ondov, Treangen et al., 2016 with default parameters against the complete RefSeq genomes and plasmids database, to search for potential contaminants. Both MetaPhlAn2 and the top hits from mash screen found the correct bacterial genomes from the three strains used in this study, together with reads from yeast (S. cerevisiae S288c). Additionally, low amounts of the phage PhiX174 were reported by mash screen. To accurately quantify the presence of these potential contaminants in our samples, and to quantify gene expression from the cultured bacteria, we mapped the non-rRNA reads to these five species using Bowtie2 (Langmead & Salzberg, 2012) with default parameters except for "-X 800" to allow for longer insert sizes. The reference genomes used are the following: GCF_000156535.1_ASM15653v1_genomic.fna (RI), GCF_000157975.1_ASM15797v1_genomic.fna (BH), GCF_000162015.1_ASM16201v1_genomic.fna (FP), CF_000146045.2_R64_genomic.fna (yeast) and NC_001422.1 (PhiX174). Gene expression was quantified using the htseq-count Python script (Anders, Pyl et al., 2015) (with parameter -a 2 to exclude multimapping reads) for all species using their available *.gff reference annotation files. Given the small size of the PhiX174, we quantified the reads mapping to its entire genome rather than its gene expression. Differential gene expression analysis of the three cultured strains was performed with DESeq2 (Love, Huber et al., 2014). To remove the effect of the different bacterial compositions in the tri-culture samples, we extracted the reads from each strain prior to the differential expression analysis, and analyzed each strain separately. In the DESeq2 design formula we included two factors: type of culture (mono-or tri-culture) and time (3 h, 6 h and 15 h). The results of the differential expression analyses were computed using a Wald test of the triculture versus the mono-culture samples.
For each strain, we extracted the genes significantly changing their expression (with Benjamini-Hochberg adjusted p-value < 0.05) in tri-culture and mapped them to different functional annotations downloaded from the IMG database (Markowitz, Chen et al., 2012): COG categories, COG numbers and KO numbers. The RNA-seq data processing code is available on GitHub (https://github.com/vllorens/syntheticGutCommunity).    (D-F) The dependency of the end point abundances of the three species on the lag phase and initial abundance of FP and RI is shown with simulations. This dependency is nonlinear, especially for initial abundances of FP and RI, illustrating that dominance in batch is highly sensitive to initial conditions. The simulations were carried out with the model parameterized on mono-and bi-culture data. For the simulations in (A), the initial abundances of RI, FP and BH were set to 0.58, 0.04 and 0.21, respectively, whereas for the simulations in (B), the lag phases for RI, FP and BH were set to 0.33, 0.08 and 0.1, respectively. These initial abundance and lag phase values represent the averages of observed initial abundances and estimated lag phases across all tri-culture experiments.  stand for Roseburia intestinalis, Faecalibacterium prausnitzii and Blautia hydrogenotrophica, respectively. The taxon Lachnospiraceae incertae sedis contains BH. High relative abundances of potential contaminants (> 10%) were found in one RI monoculture (RI_16), one RI-FP co-culture (RI_FP_9) and one triculture (RI_FP_BH_10).

Supplementary Figure Legends
Supplementary Figure 3 Test for viral, prokaryotic and eukaryotic contamination in RNA-seq data. The RNA of two non-prokaryotic organisms reached noticeable abundances: the bacteriophage phiX174, which is used as a control in Illumina sequencing, and the yeast S. cerevisiae S288c, which probably came from the yeast extract employed in the medium. In most of the samples, these potential contaminants have transcriptome-size-corrected abundances below 5%, but in one sample (Blautia hydrogenotrophica at 3 h) the yeast RNA abundance reached 18%.
Taxonomic assignment was carried out with MetaPhlAn2 and mash screen against the complete RefSeq genomes and plasmids database. Total read counts were corrected for transcriptome size (genome size in the case of the bacteriophage). The bacterial contamination observed in the Roseburia intestinalis mono-culture with 16S rRNA gene sequencing was not confirmed with RNA-seq. Row 1 refers to experiments RI_16, FP_14, BH_16 and RI_FP_BH_14, while row 2 refers to experiments RI_15, FP_15, BH_15 and RI_FP_BH_15. Supplementary

Supplementary Table Legends
Supplementary Table 1 A summary of all fermentation experiments is provided. Fermentation experiments selected for RNA-seq or model parameterization are marked accordingly.
Supplementary Table 2 The table provides parameter values after model parameterization on monocultures only and on mono-and bi-cultures (sheet 1), lag phases estimated from experiments (sheet 2) and predicted maximal abundances for the two parameterizations (sheet 3).
Supplementary Table 3  Table 3 lists genes that are significantly up-or down-regulated in tri-culture as compared to mono-culture for all three strains across all time points.
Supplementary Table 4 The sequences for species-specific primers and probes as well as their specificities are provided.