Deep Learning Enables Design of Multifunctional Synthetic Human Gut Microbiome Dynamics

Predicting the dynamics and functions of microbiomes constructed from the bottom-up is a key challenge in exploiting them to our benefit. Current ordinary differential equation-based models fail to capture complex behaviors that fall outside of a predetermined ecological theory and do not scale well with increasing community complexity and in considering multiple functions. We develop and apply a long short-term memory (LSTM) framework to advance our understanding of community assembly and health-relevant metabolite production using a synthetic human gut community. A mainstay of deep learning, the LSTM learns a high dimensional data-driven non-linear dynamical system model used to design communities with desired metabolite profiles. We show that the LSTM model can outperform the widely used generalized Lotka-Volterra model. We build methods decipher microbe-microbe and microbe-metabolite interactions from an otherwise black-box model. These methods highlight that Actinobacteria, Firmicutes and Proteobacteria are significant drivers of metabolite production whereas Bacteroides shape community dynamics. We use the LSTM model to navigate a large multidimensional functional landscape to identify communities with unique health-relevant metabolite profiles and temporal behaviors. In sum, the accuracy of the LSTM model can be exploited for experimental planning and to guide the design of synthetic microbiomes with target dynamic functions.


34
Microbial communities perform chemical and physical transformations to shape the properties of nearly every en-35 vironment on Earth from driving biogeochemical cycles to mediating human health and disease. These functions 36 performed by microbial communities are shaped by a multitude of abiotic and biotic interactions and vary as a func-37 tion of space and time. The complex dynamics of microbial communities are influenced by pairwise and higher-order 38 interactions, wherein interactions between pairs of species can be modified by other community members [1,2,3]. 39 In addition, the interactions between community members can change as a function of time due to variation in the 40 abiotic environment as well as environmental modification by the microbial community [4]. Therefore, flexible mod-41 eling frameworks that can capture the complex and temporally changing interactions that determine the dynamic 42 behaviors of microbiomes are needed. These predictive modeling frameworks could be used to guide the design of 43 precise interventions to manipulate community-level functions to our benefit.  Figure 1: LSTM model can predict the temporal changes in species abundance in a 12-member synthetic human gut community in response to dilution. (a) Proposed LSTM modeling methodology for the dynamic prediction of species abundance in a microbial community. The initial abundance information is an input to the first LSTM cell, the output of which is trained to predict abundance at the next time point. Consequently, the predicted abundance becomes an input to another LSTM cell with shared weights to predict the abundance at the subsequent time point. The process is repeated until measurements at all time points are available. (b) Scatter plot of measured (true) and predicted species abundance of a 12-member synthetic human gut community at 12 hr (N = 876, p-value = 2.44e − 257). (c) Scatter plot of measured (true) and predicted abundance at 24 hr (p-value = 6.51e − 257). (d) Scatter plot of measured (true) and predicted abundance at 36 hr (p-value = 7.42e − 257). (e) Scatter plot of measured (true) and predicted abundance at 48 hr (p-value = 1.66e − 227). (f ) Scatter plot of measured (true) and predicted abundance at 60 hr (p-value = 3.39e − 227). but also to identify what type of datasets are required for building predictive models of high richness community  LSTM framework based on our success in predicting community dynamics combined with the ease of incorporating 168 additional output variables. Therefore, we applied the LSTM framework to design health-relevant metabolite profiles 169 using synthetic human gut communities. 170 A core function of gut microbiota is to transform complex dietary substrates into fermentation end products such 171 as the beneficial metabolite butyrate, which is a major determinant of gut homeostasis [29]. In a previous study, 172 we designed butyrate-producing synthetic human gut microbiomes from a set of 25 prevalent and diverse human 173 gut bacteria using a hybrid gLV and statistical model. This hybrid model consists of a gLV model for predicting 174 community assembly and a linear regression model with interactions to predict butyrate production from species 175 absolute abundance at a given time point [25]. While the hybrid model approach was successful for predicting 176 butyrate concentration, designing community-level metabolite profiles rather than optimizing the concentration of 177 a single metabolite adds substantial complexity and limited flexibility using the hybrid modeling approach. Thus, 178 we leveraged the accuracy and flexibility of LSTM models to design the metabolite profiles of synthetic human gut 179 microbiomes. We focused on the fermentation products butyrate, acetate, succinate, and lactate which play important 180 roles in the gut microbiome's impact on host health and interactions with constituent community members [10]. 181 We used the species abundance and metabolite concentrations from our previous work [25] to train an LSTM 182 model. This model uses a feed-forward network (FFN) at the output of the final LSTM unit that maps the end-183 point species abundance to the concentrations of the four metabolites (Fig. 3a). The entire neural-network model 184 comprising LSTM units and a feed-forward network is learned in an end-to-end manner during the training process, 185 (i.e., all the network weights are trained simultaneously). Cross-validation of this model (Model M1, Table S1) on a set of hold-out community observations shows good agreement between the model predictions and experimental 187 measurements for metabolite concentrations and microbial species abundances (Fig. S1). Thus, we used this model 188 to design species-rich (i.e. >10 species) microbial communities with tailored metabolite profiles (Fig. 3a). 189 We first used the LSTM model M1 to simulate every possible combination of >10 species (26,434,916 commu-190 nities). The simulated communities separate into two regions: one centered around a dense ellipse of high butyrate 191 concentration characterized by communities containing the butyrate-producing species Anaerostipes caccae (AC) 192 and a second dense ellipse of low butyrate concentration characterized by communities lacking AC (Fig. 3b). This 193 bimodality due to the presence/absence of AC is consistent with our previous finding that AC is the strongest driver 194 of butyrate production in this system [25]. In addition, the strong negative correlation between lactate and butyrate 195 in the AC+ ellipse (R 2 = 0.72, p < 0.001, N=14,198,086) is consistent with the ability of AC to convert lactate into 196 butyrate [25]. These results demonstrate that the LSTM model can capture the major microbial drivers of metabolite 197 production as well as correlations between different metabolites. 198 We used our simulated metabolite production landscape to plan informative experiments for testing the capabili-199 ties of our model. First, we designed a set of "distributed" communities that spanned the range of typical metabolite 200 concentrations predicted by our model. To this end, we selected 100 communities that fell closest to the centroids of 201 100 clusters determined using k-means clustering of the 4-dimensional metabolite space. Second, we designed a set 202 of communities to test our model's ability to predict extreme shifts in metabolite outputs. To do so, we identified 203 four "corners" of the distribution in the lactate and butyrate space (Fig. 3b). We next examined the relationship 204 between acetate and succinate within each of these corners and found that the distributions varied depending on the 205 given corner (Fig. 3b, inset). The total carbon concentration in the fermentation end products across all predicted 206 communities displayed a narrow distribution (mean 316 mM, standard deviation 20 mM, Fig. S2). The production 207 of the four metabolites are coupled due to the structure of the metabolic networks and fundamental stoichiometric 208 constraints [30]. Therefore, the model learned the inherent "trade-off" relationships between these fermentation 209 products based on the patterns in our data. We chose a final set of "corner" communities for experimental validation 210 by choosing 5 communities from each combination of maximizing or minimizing each metabolite (80 communities 211 total, see Methods for details).

212
By experimentally characterizing the 180 designed communities, we found that the LSTM model M1 accurately 213 predicted the rank order of metabolite concentrations and microbial species abundances, substantially outperforming 214 a composite model (gLV and regression) trained on the same data for the majority (59%) of output variables 215 (Fig. S3a). Notably, the LSTM model prediction accuracy for the metabolites was similar for both the "distributed" 216 and "corner" communities ( Fig. S3b-e). These results indicate that our model is useful for designing communities 217 with a broad range of metabolite profiles that includes the extremes of the distributions. To understand how well 218 our model could separate groups of communities with extreme behaviors, we treated the "corners" as classes and 219 quantified the classification accuracy of our model. The model accurately classified the communities when considering 220 only butyrate and lactate concentrations. However, the model had poorer separation when acetate and succinate were 221 also considered in defining the classes (Fig. S3f ). The misclassification rate was higher for small Euclidean distances 222 between classes and decreased with the Euclidean distance (Fig. S3g). This implies that the insufficient variation 223 in concentrations due to fundamental stoichiometric constraints limited our ability to define 16 distinct classes that 224 maximized/minimized each metabolite. While model M1 accurately predicted metabolite concentrations and the 225 majority of species abundances, the predictions of several individual species were still quite poor (R 2 = 0 − 0.6, 226 Fig. S3a). Thus, we used the dataset to improve the model. To this end, we combined the new observations with 227 the original observations and randomly partitioned the data into 90% for training and 10% for cross-validation. The 228 resulting model (M2 , Table S1) was substantially more predictive of species abundances (R 2 > 0.5 for all but five 229 species FP, RI, CA, BA, CH (Fig. 3c).

230
One of the commonly noted limitations of machine learning models is their lack of interpretability for extracting 231 biological information about the system. Thus, we used our predictive LSTM model to decipher key relationships 232 among variables to deepen our biological understanding of the system. We used local interpretable model-agnostic instances to generate networks that provided key insights into microbe-metabolite (Fig. 3d) and microbe-microbe 237 (Fig. 3e) interactions. In general, these networks represent broad design principles for community metabolic output 238 by indicating which species have the most consistent and strong impacts on each metabolite and species abundance 239 across a wide range of sub-communities. For instance, the metabolite network highlights AC as having the largest 240 positive effect on butyrate production with additional positive contribution from EL and negative contribution from 241 DP, consistent with the previous hybrid gLV model of butyrate production by this community [25]. Additionally, the 242 number of microbial species impacting each metabolite in these networks trended with the number of microbial species 243 in the system that individually produced or consumed each metabolite (Fig. S4). For example, butyrate displayed 244 the fewest edges (3) and was produced by the lowest number of individual species (4). By contrast, acetate had the 245 most edges (6) and was produced by the largest number of individual species (19). The inferred microbe-metabolite 246 network consisted of diverse species including Proteobacteria (DP), Actinobacteria (BA, BP, EL), Firmicutes (AC, 247 ER, DL) and one member of Bacteroidetes (PC) but excluded members of Bacteroides. Therefore, while Bacteroides 248 exhibited high abundance in many of the communities, they did not substantially impact the measured metabolite 249 profiles but instead modulated species growth and thus community assembly (Fig. 3e).

250
The LIME explanations of inter-species interactions exhibited a statistically significant correlation with their 251 corresponding inter-species interaction parameters from a previously parameterized gLV model of this system [25] 252 (Fig. S5). The sign of the interaction was consistent in 80% of the interactions with substantial magnitude (> 0.05 253 in both the LIME explanations and gLV parameters). This consistency with previous observations suggests that the 254 LSTM model was able to capture the same broad trends in interspecies relationships as gLV (interpreted through 255 the average LIME explanation across all observed communities). The LSTM model captured more nuanced context-256 specific behaviors (interpreted as the LIME explanation for one specific community context) than the mathematically 257 restricted gLV model, which substantially improved the LSTM model's predictive capabilities. These results demon-258 strate that the LSTM framework is useful for developing high accuracy predictive models for the design of precise 259 community-level metabolite profiles. Our approach also preserves the ability to decipher different types of inter-  we computed the hold-out prediction performance (R 2 ) as a function of the size of the training set by sub-sampling 271 the data (Fig. 4a). We used 20-fold cross-validation to predict metabolite concentrations and species abundance.

272
Our results show that the ability to improve prediction accuracy as a function of the size of the training data set 273 was limited by the variance in species abundance in the training set (Fig. S6). For instance, certain species with 274 low variance (e.g. FP, EL, DP, RI) in abundance in the training set also displayed low sensitivity to the amount of 275 training data. The high sensitivity of specific metabolites (e.g. lactate) and species (e.g. AC, BH) to the amount of 276 training data indicates that further data collection would likely improve the model's prediction performance.

277
To determine how pairwise combinations of species impacted model prediction performance, we used 20-fold cross-278 validation to evaluate the prediction performance (R 2 ) on subsets of the total dataset, where subsets were selected 279 based on the presence of individual species or pairs of species (Fig. 4b). Using this approach, we identified individual 280 species and species pairs that had the greatest impact on the prediction performance of metabolite concentrations.

281
Sample subsets with poor prediction performance highlight individual species and species pairs whose presence 282 reduces the model's ability to make accurate predictions of final metabolite concentrations. Although the subsets 283 were much smaller than the total data set (n = 761), calculation of prediction performance was not limited by small 284 sample sizes, where the number of communities in each subset ranged from n = 77 to n = 478.

285
The interaction network shown in figure Fig. 3d shows the impact of individual species on each metabolite, but 286 does not provide information about whether the effect is due to individual species or pairwise interactions. To deter-287 mine whether pairwise interactions influence metabolite concentrations, we quantified how prediction performance 288 changed in response to the presence individual species and pairs of species. Specifically, if prediction performance 289 taken over a subset of communities containing a given species pair was markedly different than prediction performance 290 for the subsets corresponding to the individual species, this suggests the pairwise interaction impacts on metabolite 291 production. Using equation 4 (Methods), we found that the prediction performance of lactate and butyrate was 292 the least sensitive to species pairs (average decrease in prediction performance for subsets with species pairs of 0.72% 293 and 1.10% compared to corresponding single species subsets). However, the prediction performance of acetate and 294 succinate was the most sensitive to the presence of species pairs (increase in prediction performance of 6.68% for 295 acetate and a decrease of 2.951% for succinate). This difference in prediction performance suggests that pairwise 296 interactions influences the production of acetate and succinate, while the production of lactate and butyrate are 297 primarily driven by the action of individual species. The sensitivity of acetate and succinate to pairwise interactions is consistent with the inferred interaction network shown in Fig. 3d, which highlights multiple species-metabolite 299 interactions for acetate and succinate and a smaller number of strong species-metabolite interactions for butyrate 300 and lactate.

301
Pairs of certain Bacteroides and butyrate producers including BY-RI, BU-RI, and BY-AC resulted in reduced 302 prediction performance of acetate. This suggests that interactions between specific Bacteroides and butyrate pro-303 ducers were important for acetate transformations, which is consistent with the conversion of acetate into butyrate.

304
Based on the LIME analysis in Fig. 3d, AC, DP, and BP have the largest impact on lactate. Thus, the hold-out 305 prediction performance for lactate was primarily impacted by specific pairs that include these species. In sum, 306 these results demonstrate how the model can be used to identify informative experiments for investigating poorly 307 understood species and interactions between species, where collection of more data would likely improve prediction 308 performance.

309
Dynamic measurements of communities reveal design rules for qualitatively distinct 310 metabolite trajectories 311 We next leveraged the LSTM model's dynamic capabilities to understand the temporal changes in metabolite con-  (Fig. 5a). We analyzed the dynamic behavior of these communities using a clustering technique to extract 316 high level design rules of species presence/absence that determined qualitatively distinct temporal metabolite tra-317 jectories (i.e. broad trends consistent across a set of communities) and exploited the LSTM framework to identify 318 context-specific impacts of species on metabolite production (i.e. a more fine-tuned case-by-case analysis).

319
The temporal trajectories of species abundance and metabolite concentrations showed a wide range qualitatively 320 distinct trends across the 95 communities ( Fig. 5b-g). For example, some metabolites concentrations monotonically 321 increased (e.g. butyrate in Fig. 5b,c,e,g), monotonically decreased (e.g. lactate in Fig. 5b,c) or exhibited biphasic 322 dynamics (e.g. acetate in Fig. 5c). To determine if there were communities with similar temporal changes in 323 metabolite concentrations, we clustered communities using a minimal spanning tree [32] on the Euclidean distance 324 between the metabolite trajectories of each pair of communities (Fig. 5a). The resulting six clusters exhibited 325 high quantitative within-cluster similarity and qualitatively distinct metabolite trajectories ( Fig. 5b-g). Clusters 326 4 and 5 which contained the largest number of communities had a high fraction of "distributed" communities 327 (Fig. 3b). Clusters with a smaller number of communities contained a higher percentage of "corner" communities 328 (Fig. S7b,c). Therefore, the use of LSTM results from an initial experiment to identify "corner" communities 329 elucidated communities with qualitatively distinct temporal behaviors. These communities were unlikely to be 330 discovered via random sampling of sub-communities due to the high density of points towards the center of the 331 distribution and low density in the tails of the distribution (Fig. 3b). Additionally, some "corner" communities that 332 were similar in metabolite profiles when considering the end-point measurement separated into different clusters when 333 considering the dynamic data (e.g. Clusters 2 and 3, which have similar metabolite profiles at 48 hr but qualitatively 334 distinct dynamics (Fig. 5b). This demonstrates that using a community design approach to explore the extremes 335 of system behaviors with a limited time resolution enabled the identification of additional distinct behaviors when 336 the extreme communities were characterized with higher time resolution.

337
To identify general patterns in species presence/absence of these communities that could explain the temporal 338 behaviors of each cluster, we used a decision tree analysis to identify an interpretable classification scheme (Fig. S7d).

339
Using this approach, we observed that the large clusters were separated by relatively simple classification rules (i.e. 340 AC+ for cluster 4 and AC-for cluster 5), whereas the smaller clusters had more complex classification rules involving 341 larger combinations of species (3-7 species), all involving AC, DP, and DL (Fig. 5a). The influential role of DP 342 was corroborated by a previous study showing that DP substantially inhibits butyrate production [25]. In addition, 343 the inferred microbe-metabolite networks based on the LSTM model M2 demonstrated that the presence of DL 344 was linked to higher acetate and lower succinate production (Fig. 3d), consistent with its key role in shaping 345 metabolite dynamics in this system. The variation in the number of communities across clusters is consistent with 346 previous observations that species-rich microbial communities tend towards similar behavior(s) (e.g. Clusters 4 and 347 5 contained many communities). By contrast, more complex design criteria are required to identify communities 348 that deviate from this typical behavior (e.g. Clusters 1-3 and 6 contained few communities) [25].

349
While our clustering analysis identified general design rules for metabolite trajectories, there remained unex-350 plained within-cluster variation. Thus, we used the LSTM framework to identify those effects beyond these general 351 species presence/absence rules that determine the precise metabolite trajectory of a given community. Simultaneous 352 predictions of species abundance and the concentration of all four metabolites at all time points necessitates specific 353 modifications to the LSTM architecture shown in Fig. 1a. In particular, we consider a 29-dimensional input vector 354 whose first 25 components correspond to the species abundance, while the remaining 4 components correspond to the 355 concentration of metabolites (Fig. 5h). The 29-dimensional feature vector is suitably normalized so that the different  (Fig. 5i). The prediction accuracy of species abundance was lower than metabolite 366 concentrations, presumably due to the limited number of training set observations of each species (Fig. S8). 367 We used a a gradient-based sensitivity analysis of the LSTM model M3 to provide biological insights into the 368 contributions of each species on the temporal changes in metabolite concentrations (Fig. 5h,j, Methods). This 369 method involves computing partial derivatives of output variables of interest with respect to input variables, which are 370 readily available through a single backpropagation pass [33, 34]. As an example case, we applied this analysis approach 371 to the full 25-species community, which was grouped into Cluster 4, with the design rule "AC+" (Fig. 5a). Consistent 372 with this design rule, we observed strong sensitivity gradients between the abundance of AC and the concentrations of  Table S3). Dashed line indicates R 2 = 0.5, which is used as a cutoff for including a variable in the subsequent network diagrams. (d) and (e) Network representation of median LIME explanations of the LSTM model M2 from (c) for prediction of each metabolite concentration (d) or species abundance (e) by the presence of each species. Edge widths are proportional to the median LIME explanation across all communities from (b) used to train the model in units of concentration (for (d)) or normalized to the species' self-impact (for (e)). Only explanations for those variables where the cross-validated predictions had R 2 > 0.5 are shown. Networks were simplified by using lower thresholds for edge width (5 mM for (d), 0.2 for (e)). Red and blue edges indicate positive and negative contributions, respectively. (a) Sensitivity of metabolite prediction performance (R 2 ) to the size of the training dataset. Training datasets were randomly subsampled 30 times using 50% to 100% of the total dataset in increments of 10%. Each subsampled training set was subject to 20-fold cross-validation to assess prediction performance. Lineplot of the mean prediction performance over the 30 trials for each percentage of the data. Error bars denote 1 s.d. from the mean. (b) Schematic scatter plot representing how communities containing species A and B define a poorly predicted subsample of the full sample set (c) Heatmap of prediction performance (R 2 ) of acetate for each subset of communities containing a given species (diagonal elements) or pair of species (off-diagonal elements). (d) Heatmap of prediction performance for acetate, butyrate, lactate, and succinate. A sample subset containing a given species or pair of species included all communities in which the species were initially present. Predictions for each community were determined using 20-fold cross validation so that for each model the predicted samples were excluded from the training samples. N and p-values are reported in Table  S3.   X i (t n ) C j (t n+1 ) X i (t n ) C j (t n+1 ) h Figure 5: Community metabolite trajectories cluster into qualitatively distinct groups which can be classified based on presence and absence of key microbial species. (a) Schematic of experiment and network representing a minimal spanning tree across the 95 communities where weights (indicated by edge length) are equal to the Euclidean distance between the metabolite trajectories for each community. Node colors indicate clusters determined as described in the Methods. Red node with black outline annotated with "25" represents the community of all 25 species. Annotations indicate the most specific microbial species presence/absence rules that describe most data points in the cluster of the corresponding color as determined by a decision tree classifier (Methods). Communities that deviate from the rules for their cluster are indicated with a border matching the color of the closest cluster whose rules they do follow. Network visualization generated using the draw kamada kawai function in networkx (v2.1) for Python 3. . Therefore, the microbe-metabolite interaction network could be used to identify key species that 409 could be targeted for manipulating the dynamics of specific metabolites. 410 We performed time-resolved measurements of metabolite production and species abundance using a set of designed 411 communities and demonstrated that communities tend towards a typical dynamic behavior (i.e. Clusters 4 and 5).

412
Therefore, random sampling of sub-communities from the 25-member system would likely exhibit behaviors similar 413 to Clusters 4 and 5. We used the LSTM model to identify "corner cases" communities that produce metabolite 414 concentrations near the tails of the metabolite distributions at a single time point. Thus, the model allowed us to 415 identify unique sub-clusters with disparate dynamic behaviors. We demonstrated that the endpoint model predictions 416 were confirmatory (Fig. 3c) and also led to new discoveries when additional measurements were made in the time 417 dimension. Specifically, we found that some "corner cases" communities identified based on prediction of a single 418 time-point displayed distinct dynamic trajectories. For instance, Clusters 2 and 3 based on the decision tree classifier 419 displayed similar end-point metabolite concentrations (Fig. 5c,d). However, lactate decreased immediately over time 420 in Cluster 2 communities but remained high until approximately 30 hr and then decreased in Cluster 3 communities.

421
The design rule for Cluster 3 included the presence of lactate producers BU and DL (Fig. S4), suggesting that these 422 individual species' lactate producing capabilities enabled the community to maintain a high lactate concentration for    Table S1 and permanent stocks of each were stored in 25% glycerol at −80 • C as previously described [25]. Batches 474 of single-use glycerol stocks were produced for each strain by first growing a culture from the permanent stock in  Table S1 to a total volume of 5 mL (multiple tubes inoculated if more preculture volume 482 needed) for stationary incubation at 37 • C for the preculture incubation time listed in Table S1. All experiments were  Hasselbach equation to fresh media with a pH ranging between 3 to 11 measured using a standard electro-chemical 504 pH probe (Mettler-Toledo). We used (1) to map the pH values to the absorbance measurements.
The parameters b and pK a were determined using a linear regression between pH and the log term for the standards 506 in the linear range of absorbance (pH between 5.2 and 11) with A max representing the absorbance of the pH 11 507 standard, A min denoting the absorbance of the pH 3 standard and A representing the absorbance of each condition. Sequencing data were used to quantify species relative abundance as described previously [25,43]. Sequencing data 562 were demultiplexed using Basespace Sequencing Hub's FastQ Generation program. Custom python scripts were used 563 for further data processing as described previously [25,43]. Paired end reads were merged using PEAR (v0.9.10) [44] 564 after which reads without forward and reverse annealing regions were filtered out. A reference database of the 565 V3-V5 16S rRNA gene sequences was created using consensus sequences from next-generation sequencing data or   Similar to any recurrent neural network, an LSTM network, too, comprises of a network of multiple LSTM units, each representing the input-output map at a time instant. Fig. 1 shows the schematic of the proposed LSTM network architecture for abundance prediction. For a microbial community comprising of N species, each LSTM unit models the dynamics at time t using the following set of equations: where h t , c t , x t are the hidden state, cell state and input abundance at time t, respectively, and i t , f t , g t , o t are 580 input, forget, cell and output gates, respectively. σ is the sigmoid function, and denotes the Hadamard product.

581
The parameters {W mn , b mn } for m, n ∈ {f, g, h, i, o} are trainable and shared across all LSTM units. The output 582 gate o t is further used to generate the abundance for next time instant as: As shown in Fig. 1, y t is fed to the LSTM unit at the next timestep (t + 1), which in turn predicts the species 584 abundance at time t + 2. The process is repeated across multiple LSTM units in order to obtain x t final . The entire 585 architecture is trained to minimize the mean-squared loss between the predicted abundance x t final and true abundance 586x t final .

587
Using Teacher forcing for intermittent time-series forecasting 588 The end-goal for the proposed LSTM-network based abundance predictor is to accurately capture the steady-state 589 (final) abundance from initial abundance. In typical LSTM networks, the output of the recurrent unit at the previous 590 timestep y t−1 is used as an input to the recurrent unit at the current timestep x t . This kind of recurrent model, while 591 has the ability to predict final abundance, is incapable to handle he one-step-ahead prediction. The problem is even 592 more critical when one tries to anticipate more than a single timestep into the future. Teacher forcing [55] entails 593 a training procedure for recurrent networks, such as LSTMs, where 'true' abundances at intermittent timesteps are 594 used to guide (like a teacher) the model to accurately anticipate one-step-ahead abundance.

595
Teacher forcing is an efficient method of training RNN models that use the ground truth from a prior time step 596 as input. This is achieved by occasionally replacing the predicted abundance y t−1 from the previous timestep with  Table S2.

628
Using LSTM Model to Design Multifunctional Communities 629 We used the LSTM model trained on previous data (Fig. 3a) to design two sets of communities: a "distributed" 630 community set and a "corner" community set. For the "distributed" community set, we first took the predicted 631 metabolite concentrations for all communities with .10 species and used k-means clustering with k = 100 (Python 3, 632 scikit-learn v0.23.1, sklearn.cluster.Kmeans function) to identify 100 cluster centroids that were distributed across 633 all of the predictions. We then found the closest community to each centroid in terms of Euclidean distance in the 634 4-dimensional metabolite concentration space. These 100 communities constituted the "distributed" community set.

635
For the "corner" community set, we first defined 4 "corners" in the lactate and butyrate concentration space by Using LIME where R 2 i is the prediction performance taken over the subset of samples containing species i, and R 2 ij is the prediction 706 performance taken over the subset of samples containing species i and j. To generate the clusters from the dynamic community observations (Fig. 5) The decision tree shown in Fig. S5d and used to produce the annotations in Fig. 5a was generated using the Deci-    Table S3. (b-e) Prediction accuracy of model M1 for the indicated metabolites. Dashed line indicates the linear regression for all data points. Legends indicates the R 2 , Spearman ρ 2 , and RMSE for communities from the "corner" set (red) or "distributed" set (blue) for each variable. Solid black lines indicate x=y. (f) Confusion matrix for classification of the "corner" communities into their specified classes (shown in Figure  3b). Values indicate the fraction of communities from each predicted class whose metabolite concentrations were closest (Euclidean distance) to the centroid of each class (Measured Class). Colored boxes indicate "sub-classes" that fall within the 4 major classes determined in the lactate and butyrate concentration space as shown in Figure 3b. (g) Scatter plot of misclassification rate between each pair of classes (values from f, fraction of communities misclassified from one class to the other) versus the Euclidean distance between the centroids of that pair of classes. Heatmap representation of qualitative agreement/disagreement between specific inter-species interactions for the same comparison as in (a). Legend describes what each color represents. Figure S5. Sensitivity of species abundance prediction performance (R 2 ) to the size of the training dataset for each species. Training datasets were randomly subsampled 30 times using 50% to 100% of the total dataset in increments of 10%. Each subsampled training set was subject to 20-fold cross-validation to assess prediction performance. Sub-plots show the mean prediction performance (error bars indicate one standard deviation) over the 30 trials for each percentage of the dataset. Subplots were sorted according to the variance in species abundance taken over the total dataset. In general, prediction performance of low variance species was less likely to improve in response to more training data. N and pvalues are reported in Table S3. Figure S6. Characteristics of the dynamic community behaviors. (a) Minimal spanning tree of a graph representation of the 180 communities characterized in Figure 3 where each node represents a community the lengths of the edges connecting the nodes represent the Euclidean distance between a pair of communities in the 4-dimensional metabolite space. The even spread of the training and validation sets across this tree demonstrate that the subset of communities characterized in the dynamic experiment was representative of all 180 communities characterized in