Predicting microbial interactions with approaches based on flux balance analysis: an evaluation

Open Access © The Author(s) 2024. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. RESEARCH Josephetal. BMC Bioinformatics (2024) 25:36 https://doi.org/10.1186/s12859-024-05651-7 BMC Bioinformatics Predicting microbial interactions withapproaches based onflux balance analysis: anevaluation Clémence Joseph1†, Haris Zafeiropoulos1† , Kristel Bernaerts2 and Karoline Faust1* Abstract Background: Given a genome-scale metabolic model (GEM) of a microorganism and criteria for optimization, flux balance analysis (FBA) predicts the optimal growth rate and its corresponding flux distribution for a specific medium. FBA has been extended to microbial consortia and thus can be used to predict interactions by comparing in-silico growth rates for coand monocultures. Although FBA-based methods for microbial interaction prediction are becoming popular, a systematic evaluation of their accuracy has not yet been performed. Results: Here, we evaluate the accuracy of FBA-based predictions of human and mouse gut bacterial interactions using growth data from the literature. For this, we collected 26 GEMs from the semi-curated AGORA database as well as four previously published curated GEMs. We tested the accuracy of three tools (COMETS, Microbiome Modeling Toolbox and MICOM) by comparing growth rates predicted in monoand co-culture to growth rates extracted from the literature and also investigated the impact of different tool settings and media. We found that except for curated GEMs, predicted growth rates and their ratios (i.e. interaction strengths) do not correlate with growth rates and interaction strengths obtained from in vitro data. Conclusions: Prediction of growth rates with FBA using semi-curated GEMs is currently not sufficiently accurate to predict interaction strengths reliably. Keywords: Flux balance analysis, Metabolic modelling, Microbial interactions Background Microorganisms interact with one another, thereby forming complex interaction networks. Knowledge of these networks is necessary to understand community dynamics and steer it toward a desired behavior. However, it is a challenge to infer microbial interaction networks from abundance data [1], and only a few such networks have been fully resolved to date experimentally (e.g. [2]). In recent years, metabolic modelling has emerged as a new technique to address this problem [3–5]. In metabolic modeling, the knowledge about a cell’s biochemistry is represented in a concise format known as a genome-scale metabolic model (GEM). A GEM incorporates †Clémence Joseph and Haris Zafeiropoulos have contributed equally to this work. *Correspondence: karoline[email protected] 1 Department of Microbiology, Immunology and Transplantation, Rega Institute for Medical Research, Laboratory of Molecular Bacteriology, KU Leuven, 3000 Leuven, Belgium 2 Department of Chemical Engineering, Chemical and Biochemical Reactor Engineering and Safety (CREaS), KU Leuven, 3001 Leuven, Belgium Page 2 of 24 Josephetal. BMC Bioinformatics (2024) 25:36 a matrix that stores the stoichiometric coefficients of the substrates and products of all the known biochemical reactions of an organism. GEMs can be constructed automatically from annotated genomes [6–8] or built in a work-intensive process of manual curation [9]. The former approach is much faster but results in GEMs of lower quality. Low-quality GEMs can contain dead-end metabolites (i.e., metabolites that are neither substrates of internal reactions nor excreted), gaps (missing reactions), missing enzymereaction links, mass or charge imbalances or futile cycles (irreversible reactions coupled in a cycle). MEMOTE is a recent tool that checks GEM quality systematically [10]. Once a GEM has been obtained, it can be analyzed to learn more about the cell’s metabolic capabilities. Flux balance analysis (FBA) is a constraint-based optimization method and a popular technique to computes fluxes through the biochemical reactions assuming that intracellular metabolites are at a steady state, i.e., the production rate of each metabolite equals its consumption rate [11]. FBA maximizes a certain objective function and returns the flux of each reaction (flux vector) when the flux of the objective function equals its optimal. However, while the optimal value of the objective function is unique, the corresponding flux vector may have an infinite number of possibilities [12]. Usually, the objective function maximizes the flux through the (artificial) biomass reaction, which represents cell maintenance and growth [13]. While maximizing for the objective function, the number of non-zero fluxes can be minimized (parsimonious FBA, abbreviated as pFBA[14]). In this case, the objective value is the same as in the standard FBA but the number of required enzymes and, thus, the cost of enzyme production are minimized.Further constraints on flux values can be specified to implement thermodynamic constraints, such as irreversibility of reactions, or to represent limited metabolite concentrations in the medium. In FBA, the medium is defined constraints on fluxes through import reactions (uptake rates). Given the stoichiometric matrix represented by the GEM, FBA predicts a flux distribution and the growth rate of an organism, where the latter is the flux through the biomass reaction. This is useful in several applications, such as medium optimization, the design of knock-out mutants producing target metabolites, and the exploration of metabolic responses to altered conditions [15]. However, because of its steady-state assumption, FBA is static by default and thus cannot model batch processes. This restriction is circumvented in dynamic FBA, which combines FBA with the differential equations of a kinetic model [16]. The kinetic model describes the changes in biomass and metabolite concentrations, which are updated with growth, consumption, and production rates computed with FBA, in turn providing new uptake rates for the next FBA iteration. Recently, several tools have been developed that rely on static or dynamic FBA variants to model microbial communities. The definition of a community objective function is still an open problem [17]. As discussed previously [18], most of the tools developed so far can be divided into three groups depending on their solution to this problem, namely (1) introduction of a group-level objective function to optimize the community growth rate [19–21], (2) optimization of the growth rate of each species independently of the others [22–25], and (3) reliance on measured abundances to adjust species growth rates [3, 19]. Reverse ecology aims to learn as much as possible about the ecology of an organism or group of organisms from their genomes [26]. Community FBA is a powerful reverse Page 3 of 24 Josephetal. BMC Bioinformatics (2024) 25:36 ecology approach capable of predicting ecological interactions between species pairs from their genomes. This is achieved by building a GEM for each species and then computing growth rates alone and in the presence of another species in silico. Following Gause’s famous strategy [27], the comparison of growth rates in monoand in coculture then elucidates the interaction sign and strength. For instance, two species that compete will both grow better in monothan in co-culture, whereas a species crossfeeding on metabolites produced by another will grow better in cothan in monoculture. FBA-based interaction prediction was previously applied to study interactions as a driver of microbial co-occurrence [28], investigate the prevalence of interaction types [4] and their resilience to nutrient change and invasion [5], as well as the synthesis of novel metabolites in the presence of interaction partners [29]. Despite this range of applications, a systematic evaluation of its accuracy has not yet been performed to date. Thus, our goal here is to assess the accuracy of interaction prediction by community FBA tools. To evaluate interaction prediction accuracy, we collected interactions measured invitro in 6 studies investigating human and mouse gut microbial communities [2, 30– 34]. GEMs were obtained either from AGORA, which is a repository of semi-refined metabolic reconstructions for gut bacteria [35] or from the literature for high-quality reconstructions [36–39]. We focus on three tools that cover a range of approaches regarding the treatment of the (community) biomass function(s) (see Table1 and Fig.1). The Microbiome Modeling Toolbox (MMT) [3] implements a pairwise screen for microbe–microbe and/ or host–microbe metabolic interactions that are inferred by determining metabolic exchanges between them. For this, the biomass functions of the species under study are both included in a merged model. Using the merged model, one species is silenced while the growth rate for the other is optimized, and vice-versa (monocultures). Then, a third optimization is performed, maximizing both growth rates simultaneously (co-culture). Given the predicted growth rates in monoand co-culture, an interaction is reported if the ratio of growth rates in monoand co-culture is above or below a user-defined threshold. MMT also offers another approach (not evaluated here), where it incorporates relative abundances from sequencing data into the community model. In this case, the community model consists of three compartments (diet—lumen—fecal), and a community biomass reaction is built based on the biomass functions of each of several species present in a sample and on their corresponding relative abundances. In both approaches, the protein-demand reactions, i.e., reactions that enable the accumulation of enzymes, are coupled with their corresponding protein recycling/utilization reactions. This way, the dependency between protein synthesis and utilization is ensured [40]. Like the community modeling module of MMT, MICOM uses relative abundances derived from amplicon or metagenomic sequencing [19] as a proxy for dry-weight taxon abundances. MICOM can be considered as an extension of the multiobjective OptCom [21], and SteadyCom [20] approaches that simultaneously maximize both the individual and the community growth rates. OptCom acknowledged the existence of a tradeoff between growth at the level of the community and its members and addressed it by solving a bilevel optimization problem (one optimization problem embedded within another) where the inner problem pertains to the individual species’ biomass production and the outer to the community and its corresponding objective function. Using Page 4 of 24 Josephetal. BMC Bioinformatics (2024) 25:36 Table 1 Overview of the interaction prediction tools and their parameters Tool Input Community approach Optimisation aproach Community biomass function Output Special capabilities Method ID Tool parameters COMETS List of GEMs along with their corresponding initial biomass for each community Dynamic Maximizes each species biomass sequentially considering their initial biomass and total uptake then updates species’ biomass and extracellular compound concentrations No Growth curves = > Growth rates and maximal biomass Spatial dimension Time dimension Chemostat or batch H param = test_tube timeStep = 1 maxCycles = 10 initial_pop = 1 g/sp objective_style = MAXIMISE_OBJECTIVE_FLUX and MAX_OBJECTIVE_MIN_ TOTAL H/10 param = test_tube timestep = 0.1 maxCycles = 20 initial_pop = 0.002 g/sp objective_style = MAXIMISE_OBJECTIVE_FLUX and MAX_OBJECTIVE_MIN_ TOTAL Page 5 of 24 Josephetal. BMC Bioinformatics (2024) 25:36 Table 1 (continued) Tool Input Community approach Optimisation aproach Community biomass function Output Special capabilities Method ID Tool parameters MICOM List of GEMs from each community (samples) Cooperative tradeoff Maximizes community growth rate (sum of the individual growth rates) then, limits community growth to only a fraction of its maximum rate to create a trade-off between optimal community growth and individual growth rate maximization Yes Growth rates Like OptCom but in a more efficient, simpler implementation. Growth rate estimates dependent on species abundance in the community under study Tradeoff Fraction = from tradeoff() function pfba = True min_growth = 0 solver = gurobi OptCom Maximizes community growth rate (as described above) by maximizing each species growth at a steady state using a series of strategies (lagrangian, minimization of metabolic adjustment (MOMA), linear MOMA (LMOMA)) Yes Growth rates Does not assume that all taxa have the same growth rate; allows for taxon-specific dilution rates lMoma, Moma, Original Strategy = (moma, lmoma, original) min_growth = 0 solver = gurobi MMT List of GEMs handled in pairwise combinations Pairwise Maximizes biomass functions of both species simultaneously using a merged model of the 2 GEMs Yes Growth rates, interaction type User-provided threshold defining what is considered a significant difference between growth rate in monoand co-culture sigD = 0.1 c = 400 u = 0 solver = gurobi Page 6 of 24 Josephetal. BMC Bioinformatics (2024) 25:36 a computationally simpler approach, MICOM assumes a (constant) growth rate μi for each species and constrains the overall community growth rate μC, which is obtained by a weighted sum of the individual species growth rates using a trade-off parameter. As there are infinitely many combinations of weights even for a single value of μC, MICOM implements a regularization by implementing an additional function over the individual growth rates μi, requiring consistency with the observed abundances. As shown by the authors, quadratic regularization (L2) fulfills these requirements. Thus, the ‘cooperative trade-off’ approach in MICOM incorporates atrade-offbetween optimal community growth (maximizing μC) and individual growth rate maximization (L2 minimization). MICOM also supports OptCom-based approaches. More specifically, one can maximize the total community biomass subject to the maximization of every species’ biomass (“original” strategy). Alternatively, one can minimize the cooperative cost, meaning the relative loss of a species to benefit the community, subject to the maximization of Fig. 1 Overview of the tools used in the evaluation. Three tools were selected for the evaluation: MICOM with the cooperative tradeoff function and different OptCom versions, MMT simulatePairwiseinteractions function, and COMETS. The input parameters are further discussed in Table 3 in Methods. The “biomass maximization” panel represents the objective function maximization for all methods. For the OptCom implementation in MICOM, the original strategy of optimisation is shown Page 7 of 24 Josephetal. BMC Bioinformatics (2024) 25:36 the community’s total biomass (minimization of metabolic adjustment; “moma” strategy [41]); the cooperative cost, in this case, is based on the sum of the subtraction of each species’ growth rate from its optimal growth. Likewise, the “lmoma” strategy minimizes the cooperative cost considering its absolute value. The COMETS package [23, 42] introduces two dimensions not considered by MMT and MICOM, namely the physical space (in two or three dimensions) and time. In COMETS, each species has an initial biomass that can optionally beplaced in certain location(s). As the environment of a community changes over time due to the consumption or secretion of compounds by the species present, COMETS simulates these changes over time through dynamic FBA [43]. At each iteration, COMETS estimates each species’ uptake bounds based on the concentration of nutrients in the medium to check that there is a sufficient amount considering the total uptake rate and optimizes each model using standard FBA. The resulting fluxes are then provided as inputs to estimate the changes in the biomass of each species as well as the concentration of the extracellular metabolites for the next iteration. Therefore, contrary to the MMT module for community modeling and the MICOM approaches, COMETS does not assume a community biomass function. Here, we evaluated the performance of MMT, MICOM and COMETS on 100 experimentally quantified ecological interactions using 29 metabolic reconstructions (including four high-quality ones) of 25 mammalian gut bacteria (see Fig.1). Table1 provides details on the three tools. Results This study evaluated the performance of three different tools (MICOM, MMT and COMETS) with different parameters, two environments including invitro media (simulations of the ones used in the actual invitro experiments) and Western diet, and two sets of models (GEMs) with different levels of curation (AGORA models and refined models). First, growth rates were calculated for each monoculture in each condition and compared to the experimental growth rate of the corresponding species to determine a correlation (Fig.2A, Additional file1: TableS1 and Additional file2: Text). Then, coculture growth rates were predicted using 13 different tool settings across the four conditions (two environments and two curation levels). Interaction strength was quantified as the ratio between either the growth rates or the maximal abundance in co-culture versus monoculture. Predicted ratios were then compared to the experimental ratios of the species growth rate or abundance in co-culture versus monoculture (Fig.2B, Additional file1: Tables S2–S3 and Additional file2: Figs. S1–S2). Here, we present the results related to the AGORA reconstructions using the in-silico representations of the corresponding media used invitro, as this is common practice [35, 44]. Further analyses using (1) the AGORA models with the Western diet medium, (2) the refined models with the invitro media simulations and (3) the refined models with the Western diet are available in the Additional file2: Figs. S1–S2, Additional file1: Tables S3–S4). Monoculture growth rates When the AGORA reconstructions were used along with in silico media that matched the experimental media, there was no significant correlation between the predicted Page 8 of 24 Josephetal. BMC Bioinformatics (2024) 25:36 and the expected growth rate (Fig.3, Additional file1: Tables S1 and S2), whereas refined GEMs grown in silico on a Western diet resulted in a significant Spearman correlation (Additional file1: TableS1, see also Additional file2: Text), which was however driven by an outlier. We also evaluated the effect of different methods (COMETS, MICOM, MMT) on prediction accuracy in the two environmental conditions. For AGORA models, the agreement between observed and predicted growth rates is low for all tools. To evaluate prediction accuracy in a species-specific manner, we calculated the mean absolute difference between predicted and observed growth rates for each species across all methods (Additional file1: TableS2). As expected, our results show that the mean absolute difference per species is lower in simulations matching invitro media than in the Western diet (in 72% of the species). Interestingly, for a few AGORA models, the difference between predicted and observed growth rates was an order of magnitude smaller than for others, such as Desulfovibrio piger (0.455) and Prevotella copri (0.576), compared to Bacteroides thetaiotaomicron (3.282) or Clostridium inoculum (3.558). Furthermore, only the refined models of Bacteroides thetaiotaomicron and Faecalibacterium prausnitzii performed better than their AGORA equivalents for both conditions (Western diet and invitro media). The perfect agreement between MICOM and MMT was expected for monocultures as both carry out standard static FBA for single species. In case of co-cultures, where the Fig. 2 Overview of the process in this analysis: The three software tools (MICOM, MMT and COMETS) were validated using 18 semi-refined models from AGORA (Am) and four refined (i.e. manually curated) models (Rm). Predictions were carried out using in vitro media matching the composition of the media in the experiments with the uptake flux bound values (mmol/gDW/h). The predicted monoculture growth rates were compared to growth rates derived from experimental growth curves as the slope of the exponential phase in log scale. Finally, ratios of growth rates or maximal abundance (the latter only in COMETS) in coand monoculture were compared to corresponding experimental data in coand monoculture to evaluate the accuracy of interaction strength prediction Page 9 of 24 Josephetal. BMC Bioinformatics (2024) 25:36 approach for the community model varies (Table1), differing growth rates for each species, and thus interactions, were predicted from each tool. Interaction strengths Next, we assessed how well metabolic modelling tools predict interaction strength. For this, we quantified interaction strength either as the ratio of the growth rate in coversus mono-culture (MMT, MICOM, COMETS) or as the ratio of the maximal abundance in coversus mono-culture (COMETS). A positive effect of an interaction partner on a species results in a positive interaction strength (ratio above one), while a negative effect gives a negative one (ratio below one). We note that this measure can result in several interaction strengths per species (as many as there are interactions in which it participates). Additional file 1: Table S3 summarizes the expected versus observed interaction strengths and Fig.3B shows their corresponding correlations. For the refined GEMs (see also Additional file2: Text), predicted interaction strengths were significantly negatively correlated to measured ones for most methods. Overall, the effect sizes (correlation values) were too small to draw a conclusion on the impact of curation level or the medium. The effect of the methods used was also assessed. Overall, when all the models are taken into account, none of the four methods performed well, with correlation values being below 0.3 (Additional file1: TableS3). Additional file1: TableS4 shows the correlation between the expected and experimental ratios for each species across all methods. Fig. 3 Experimental versus predicted growth rates in simulations matching the media used in the experiments for: A mono-cultures and B co-cultures with pairs of AGORA reconstructions using several parameters sets for the 3 software tools (see Table 3). The black line indicates where perfect matches between predictions and observations would be positioned. Spearman correlation coefficient values between expected and predicted growth rates are shown in the respective color per method; black coloring denotes the global Spearman correlation coefficient Page 16 of 24 Josephetal. BMC Bioinformatics (2024) 25:36 its corresponding AGORA model was used for the analysis. In all but one case, the ANI score of the genome used was higher than 95%, indicating that the two genomes belong to the same species (Additional file1: Tables S9 and S10). As AGORA2 was recently published [60], we repeated the ANI score calculation using the reconstructions of it (Additional file1: TableS11). For all the query genomes a different genome was found with a higher ANI score than from the ones included in AGORA. Therefore, a different reconstruction was mapped for each strain under study. We performed the MMT analysis with those models and the interaction predictions differed from those using the AGORA models but were not in statistically significant agreement with the expected ones (Additional file1: TableS12). Refined models For the refined models, we identified four GEMs that matched species in the experimental data, but only a few were available for gut bacterial species. Two models were found for the human data set, Bacteroides thetaiotaomicron (BT_iAH991) [37] and Faecalibacterium prausnitzii (iFprauz) [36], and two for the mouse data set, namely Akkermansia muciniphila (Yakk_v2)[39] and Enterococcus faecalis (enterococcus_faecalisV583) [38]. These refined models were curated by different scientists, using different names for metabolites and reactions, which made them incompatible with each other. To resolve this issue, we modified the exchange reaction and metabolite names in the refined models to match those in the AGORA models. However, someentities in the Enterococcus faecalis model had no equivalent in the AGORA model and were left unchanged. Some fields (i.e. annotations) were added to the structure of some refined GEMs to meet the expectations of MMT, which was developed to be used with AGORA models and is designed to be consistent with their structure. Medium definition Accurately defining the environment is critical for reliable metabolic predictions. To assess the impact of environmental complexity, we selected two sets of media for our analyses. The first set is used in metabolic models of human gut microbiota, but it is not representative of the selected experimental conditions. The second set is designed to mimic the in-vitro selected experimental environment. MMT medium The MMT tutorial offers four media formulations that mimic the gut metabolic environment [35]. We chose the Western diet-like medium without oxygen as the baseline for our analysis since it is compatible with AGORA models (identical reactions and metabolite names). Furthermore, all models were able to grow on the Western diet-like medium except for the refined model of Akkermansia muciniphila. In vitro media We used a total of seven complex media, namely ABB, AF, mMCB, YCAG, YCGMS, YCFA, and YCGD, which were supplied in anaerobic conditions with nitrogen (5–10%), carbon dioxide (10%), and dihydrogen (5–7%) (Table3). The exact chemical composition of these media is unknown since they include components such as peptones, amino Page 17 of 24 Josephetal. BMC Bioinformatics (2024) 25:36 acids, yeast extract, brain heart infusion, and calf serum. To ensure consistency in the composition of peptones and other amino acid mixtures, we assimilated them to have the same composition per gram of amino acid and salt, as described by Microxpress® for soya peptone. Calf serum and heart infusion were supplemented with vitamins to ensure optimal bacterial growth. For tryptone and yeast extract, we followed the protocol described by Marinos etal. [63] for GEMs and converted the flux constraints according to the initial quantity of each in the respective medium. The concentration of other components was provided either in molar or gram per liter units, which we converted to molar mass. The molar concentration (mol/L) was used as flux value (mmol/gdW*h). We modified the invitro media to ensure optimal bacterial growth, using the Western diet as a reference for most bacterial models. The yeast extract and tryptone defined in silico contained oxygen, which was not compatible with the anaerobic conditions of the experiments, and so was removed. Additionally, although selected vitamins and ions were added to the composition of some complex components, such as calf serum and heart infusion (e.g., vitamin B12, thiamin, riboflavin, biotin), it was not enough to ensure the growth of each species in silico. Therefore, vitamins were added at a quantity judged not to be limiting for bacterial growth (1g/L) in amounts equal to the ones described in the MMT Western diet medium. Furthermore, some GEMs required other specific metabolites to initiate minimal growth, such as arabinose for Eggerthella lenta, N-Acetyl-Neuraminic Acid and 5,6-Dimethylbenzimidazole for Akkermansia muciniphila and for Enterococcus faecalis. These metabolites were also added to the media at a concentration of 1g/L, following the amount described in the Western diet medium. The complete description of each medium can be found in Additional file1: TableS13. Tools The scripts for the three methods are available through a GitHub repository. The scripts were executed using MATLAB R2022a or Python 3.10, and all methods were used with Gurobi 9.5.1. Table3 shows the complete set of parameters values used. Microbiome Modeling Toolbox (MMT) The “Computation and analysis of microbe-microbe metabolic interactions” tutorial [64] was followed to predict pairwise interactions with MMT. This method takes as input a list of GEMs and their paths, as well as a list of media or dietary conditions containing the exchange reactions, the lower and the upper bounds. Then, MMT performs all pairwise co-cultures and calculates both monoculture and co-culture growth rates. The minimum growth rate difference between mono and co-culture that is considered significantly different (sigD value) was initially set to 0.1 (10%). This allows identifying the type of interaction between the two species, but this parameter was not used here since only ratios were considered for this analysis. The coupling factor (c) was kept at its default value of 400, and the threshold (u), defining the flux value allowed through reactions if biomass is null, was set to zero. Page 18 of 24 Josephetal. BMC Bioinformatics (2024) 25:36 Table 3 Overview of the parameters used for each of the methods Method ID Name Description Value in experiment MMT c Coupling factor Define how the reaction fluxes are coupled to the biomass reaction flux and defined as: flux span = − (c * flux(biomass)) to + (c * flux(biomass)) 400 SigD Significant difference Value that counts as a significant difference between monoculture and coculture growth rate 0.1 u Threshold u Flux allowed in reactions if biomass flux = 0 0 mergeGenes Gene merging variable Boolean, wether the gene are added and merge in the community mode (Time consuming) FALSE MICOM cooperative_ tradeoff cooperative_tradeoff(= fraction) Minimum proportion of the maximal community growth to allocate to species growth rate Min value where both species are growing between 0.1 and 1 with 0.1 step pfba Parsimonious FBA Define if a parsimonious FBA is performed TRUE min_growth Minimum growth rate Minimum growth rate required for each species 0 Abundance Relative abundance of each species in the community. By default each species has the same abundance 0.5 and 0.5 OptCom Strategy Strategy used to solve the optimization problem MOMA lMOMA Original min_growth Minimal growth Minimal growth required for each species 0 pfba Parsimonious FBA Define if a parsimonious FBA is performed TRUE Page 19 of 24 Josephetal. BMC Bioinformatics (2024) 25:36 MICOM We used MICOM [65]to predict community growth rates. We provided the identifiers and filenames of the species in pairs and then created a community and furnished the medium to parameterize the model. The MICOM medium is defined as a list of exchange reaction names of the community associated with the available flux of the metabolite (positive values). The monoculture growth rates were obtained by creating a Table 3 (continued) Method ID Name Description Value in experiment COMETS H/10 initial_pop Gram of biomass in the environment 0.002 time_step Time step of a FBA problem in hour 0.1 maxCycles Number of steps max 20 H initial_pop Gram of biomass in the environment 1 time_step Time step of a FBA problem in hour 1 maxCycles Number of steps max 10 BOTH defaultVmax V max value per default in mmol/g. CDW/h 18.5 defaultKm KM value per default in M (molar conc.) 0.000015 SpaceWidth size of the cell in cm3 1 maxSpaceBiomass Capacity maximum in gr. cell dry weight 100 minspaceBiomass Capacity minimum in gr. cell dry weight 1.00E−11 obj_stype Objective type Define if the strategy used to solve the optimization problem MAX_OBJ_MIN_ TOTAL = > Pfba AND MAXIMIZE_OBJECTIVE_ FLUX = > pFBA Grid Grid size Number of boxes in the x and y axis to define the grid size [1,1] Static Related to metabolites in media definition, define if the metabolites are in limited amount TRUE Parameters not shown here were kept at their default value in the calculation Page 20 of 24 Josephetal. BMC Bioinformatics (2024) 25:36 community with only one species using the function optimize_single. The optcom function was used with the three strategies lMoma, Moma, and original, with minimal growth of 0 and pFBA. The community prediction was run with the Gurobi solver. The cooperative_tradeoff function was used for the community, with a fraction value that was defined using the tradeoff function and taking the smallest tradeoff value that allowed both species to grow, without consideration for the value of the growth rate, to minimize the cooperation between species. The minimal growth rate was also zero for this function and calculated with pFBA. Finally, the growth rates were extracted from the output. COMETS COMETS implements dFBA and takes as input the metabolic model of the community and the medium, defined as a dictionary of the metabolites in the medium and their corresponding molar concentration. We followed the comestpy "Growth in a test tube" tutorial [66]. COMETS was run with two settings for mono and co-cultures: For the H (hour) condition, an initial biomass of one gram was used with a time step of one hour. This approach was expected to yield comparable outcomes to other techniques that use a biomass function computed per gram of dry weight and hour. As for the H/10 condition, the initial biomass was set to 2mg with a time step of 0.1h. COMETS was run for both settings with both FBA (obj_style: MAXIMIZE_OBJECTIVE_FLUX) and pFBA (obj_style: MAX_OBJECTIVE_MIN_TOTAL). In general, the maximum number of cycles was set to 20 but was increased for some monocultures to ensure that the stationary phase was reached. The other community parameters were the same as described in the tutorial. Predictions All three tools predicted monoculture and co-culture growth rates for 30 AGORA monocultures and 89 pairs. Additionally, monoculture predictions were made for four refined models (four for the Western diet and eight for in-vitro media) and co-culture predictions were made for 76 co-cultures using one of the two available refined models. All predictions were made using the Western diet and one of the seven media that matched the experimental conditions. Finally, the ratios were calculated based on the growth rates obtained from these predictions. Statistics For each species, the effect of growth in co-culture was quantified by computing the ratio of its growth rate in coversus monoculture. When both a curated and a semicurated GEM were available for a species, both were included in the analysis. Statistical values were computed using R 4.1.3 and Python 3.10. Spearman correlation To assess correlations based on both values and ranks, Spearman’s correlations were computed for monoculture growth rates and co-culture ratios (predicted versus experimental). For each of the four groups of combinations (Western diet and AGORA models, Western diet and refined models, in-vitro media and AGORA models and in-vitro Page 21 of 24 Josephetal. BMC Bioinformatics (2024) 25:36 media and refined models), a global correlation value was calculated without separating the species or the methods, highlighting the impact of different media and levels of refinement of the prediction. The ‘cor.test’ function in R was used to calculate Spearman correlations and their p values. Wilcoxon test The Wilcoxon test was performed with the R function wilco16xon.test. For all conditions, the mean absolute difference was calculated between the predicted and experimental growth rate for monocultures and between the predicted and experimental ratio of growth rates for cocultures. A paired Wilcoxon test was then carried out to compare AGORA and refined models in both media conditions. ROC curves We used a multiclass ROC curve to handle three interaction sign classes: positive, neutral, and negative interactions. We arbitrarily defined these classes for non-neutral effects, at 20% of the difference between monoand co-culture. Consequently, ratios under 0.8 were designated negative, ratios over 1.2 were designated positive, and ratios in between were classified as neutral. Using the sklearn Python package’s roc_curve function, we computed the true positive rate (TPR) and false positive rate (FPR) for each threshold value. The threshold represents the value above which an interaction strength is classified as belonging to the positive class and below which it is classified as belonging to the negative class. The maximum number of thresholds is the number of ratios per method plus one. We employed the same range of threshold values to compute the TPR and FPR for each class within a method. Abbreviations Am AGORA model (semi-curated GEM obtained from AGORA) AUC Area under the curve FBA Flux balance analysis GEM Genome-scale metabolic model IVm In-vitro media MMT Microbiome Modeling Toolbox pFBA Parsimonious FBA GR Growth rate Rm Refined model (curated GEM) WD Western diet Supplementary Information The online version contains supplementary material available at https:// doi. org/ 10. 1186/ s1285902405651-7. Additional file1: TableS1. Overview of statistical values for growth rates predicted in monocultures. H: time step of 1 per hour, H/10: time step of 0.1 per hour. TableS2. Mean absolute difference per species across the four methods and for the different media. TableS3. Global Spearman correlation refers to the correlation across all methods and settings. Significant p values are highlighted in grey. GR: Interaction strength computed as the ratio of growth rates, H: time step of 1 per hour, H/10: time step of 0.1 per hour, MX: interaction strength computed as the ratio of maximal abundances, Pars: Parsimonious FBA. TableS4. Correlation between predicted and measured interaction strengths for each species, with those that have refined metabolic reconstructions highlighted in grey. Interaction strength is quantified by comparing the growth rate or maximal biomass in coand monoculture. TableS5. Wilcoxon test value for monocultures and cocultures with refined versus AGORA Models. TableS6. Products returned from FBA using the AGORA and a refined model of Faecalibacterium prausnitzii using in vitro media and alterations of those. TableS7. Experimental growth rates for monoculture species and bacterial load in media. Growth rates for the ABB, YCAG, YCFA, YCGD and YCGMS media, were obtained manually. Growth rates are expressed in h−1. TableS8. Description of the GEM used for each species with its structure description and MEMOTE score with species names from Page 22 of 24 Josephetal. BMC Bioinformatics (2024) 25:36 the literature used for experimental data and modelling and their updated names (nomenclature update 12/22). TableS9. Mouse strains and their matched AGORA metabolic models. The selection was based on the highest ANI score between each mouse strain genome and the AGORA models’ corresponding genomes. TableS10. ANI score for alignment of mice species against AGORA model species (only top 5 ANI scores shown). ANI scores below 80% were not considered in this analysis. TableS11. ANI scores of mice species against AGORA2 model species. With green highlight the genomes included in AGORA2 that have the highest ANI score with the genome of interest, with orange the corresponding ones included in AGORA v1. TableS12. Comparison of the growth rates returned from MMT using pairs of mice AGORA and AGORA2 models. R1 stands for the growth of species A and R2 for the one of species B. TableS13. In silico media description (IVm and WD) where values represent the lower bounds of uptake fluxes. The lower bounds are set as the mmol/L concentration of every metabolite in the experimental medium assuming unit values for biomass and time. Additional file2: Text. Findings regarding the experiments using the AGORA models along with the Western diet and the refined models with both the in vitro media and the western diet. Supplementary Figures are included. Acknowledgements We thank Daniel Rios Garza for the helpful discussions. Author contributions KF designed the study, CJ and HZ prepared GEMs, ran tools, and analyzed results, which were critically discussed with KF and KB. All four authors contributed to the manuscript. Funding This work was supported by funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program under Grant Agreement Number 801747 (EcoBox) and Grant Agreement Number 101000309 (3D’omics). Availability of data and materials Other informations can be found in the GitHub repository: https:// github. com/ ClemJos/ pairw ise_ tool_ comp which contains the experimental and predicted ratio and growth rate values, the scripts for MICOM, COMTES and MMT, the ROC curves produces for all species methods for all conditions. Declarations Ethics approval and consent to participate Not applicable. Consent for publication Not applicable. Competing interests The authors declare that there are no competing interests. Received: 23 March 2023 Accepted: 11 January 2024 References 1. Faust K. Open challenges for microbial network construction and analysis. ISME J. 2021;15:3111–8. 2. Venturelli OS, Carr AV, Fisher G, Hsu RH, Lau R, Bowen BP, et al. Deciphering microbial interactions in synthetic human gut microbiome communities. Mol Syst Biol. 2018;14:e8157. 3. Baldini F, Heinken A, Heirendt L, Magnusdottir S, Fleming RMT, Thiele I. The Microbiome Modeling Toolbox: from microbial interactions to personalized microbial communities. Bioinformatics. 2019;35:2332–4. 4. Freilich S, Zarecki R, Eilam O, Segal ES, Henry CS, Kupiec M, et al. Competitive and cooperative metabolic interactions in bacterial communities. Nat Commun. 2011;2:589. 5. Machado D, Maistrenko OM, Andrejev S, Kim Y, Bork P, Patil KR, et al. Polarization of microbial communities between competitive and cooperative metabolism. Nat Ecol Evol. 2020. https:// doi. org/ 10. 1101/ 2020. 01. 28. 922583. 6. Machado D, Andrejev S, Tramontano M, Patil KR. Fast automated reconstruction of genome-scale metabolic models for microbial species and communities. Nucleic Acids Res. 2018;46:7542–53. 7. Seaver SMD, Liu F, Zhang Q, Jeffryes J, Faria JP, Edirisinghe JN, et al. The ModelSEED biochemistry database for the integration of metabolic annotations and the reconstruction, comparison and analysis of metabolic models for plants, fungi and microbes. Nucleic Acids Res. 2021. https:// doi. org/ 10. 1093/ nar/ gkaa1 143. 8. Zimmermann J, Kaleta C, Waschina S. gapseq: informed prediction of bacterial metabolic pathways and reconstruction of accurate metabolic models. Genome Biol. 2021;22:1–35. 9. Thiele I, Palsson BØ. A protocol for generating a high-quality genome-scale metabolic reconstruction. Nat Protoc. 2010;5:93–121. 10. Lieven C, Beber ME, Olivier BG, Bergmann FT, Ataman M, Babaei P, et al. MEMOTE for standardized genome-scale metabolic model testing. Nat Biotechnol. 2020;38:272–6. 11. Orth JD, Thiele I, Palsson BØ. What is flux balance analysis? Nat Biotechnol. 2010;28:245–8. Page 23 of 24 Josephetal. BMC Bioinformatics (2024) 25:36 12. Smallbone K, Simeonidis E. Flux balance analysis: a geometric perspective. J Theor Biol. 2009;258:311–5. 13. Feist AM, Palsson BO. The biomass objective function. Curr Opin Microbiol. 2010;13:344–9. 14. Lewis NE, Hixson KK, Conrad TM, Lerman JA, Charusanti P, Polpitiya AD, et al. Omic data from evolved E. coli are consistent with computed optimal growth from genome-scale models. Mol Syst Biol. 2010;6:390. 15. Raman K, Chandra N. Flux balance analysis of biological systems: applications and challenges. Brief Bioinform. 2009;10:435–49. 16. Henson MA, Hanly TJ. Dynamic flux balance analysis for synthetic microbial communities. IET Systems Biol. 2014;8(5):214–29. 17. Perez-Garcia O, Lear G, Singhal N. Metabolic network modeling of microbial interactions in natural and engineered environmental systems. Front Microbiol. 2016;7:673. 18. Garza DR, Gonze D, Zafeiropoulos H, Liu B, Faust K. Metabolic models of human gut microbiota: advances and challenges. Cell Syst. 2023;14:109–21. 19. Diener C, Gibbons SM, Resendis-Antonio O. MICOM: metagenome-scale modeling to infer metabolic interactions in the gut microbiota. mSystems. 2020;5:e00606–19. 20. Chan SHJ, Simons MN, Maranas CD. SteadyCom: predicting microbial abundances while ensuring community stability. PLOS Comput Biol. 2017;13:e1005539. 21. Zomorrodi AR, Maranas CD. OptCom: a multi-level optimization framework for the metabolic modeling and analysis of microbial communities. PLoS Comput Biol. 2012;8:e1002363. 22. Bauer E, Zimmermann J, Baldini F, Thiele I, Kaleta C. BacArena: individual-based metabolic modeling of heterogeneous microbes in complex communities. PLOS Comput Biol. 2017;13:e1005544. 23. Harcombe WR, Riehl WJ, Dukovski I, Granger BR, Betts A, Lang AH, et al. Metabolic resource allocation in individual microbes determines ecosystem interactions and spatial dynamics. Cell Rep. 2014;7:1104–15. 24. Popp D, Centler F. μBialSim: constraint-based dynamic simulation of complex microbiomes. Front Bioeng Biotechnol. 2020;8:574. 25. Zhuang K, Izallalen M, Mouser P, Richter H, Risso C, Mahadevan R, et al. Genome-scale dynamic modeling of the competition between Rhodoferax and Geobacter in anoxic subsurface environments. ISME J. 2011;5:305–16. 26. Levy R, Carr R, Kreimer A, Freilich S, Borenstein E. NetCooperate: a network-based tool for inferring host-microbe and microbe-microbe cooperation. BMC Bioinform. 2015;16:164. 27. The struggle for existence. In: The science of the struggle for existence. Cambridge University Press; 2003. p. 1–26. 28. Zelezniak A, Andrejev S, Ponomarova O, Mende DR, Bork P, Patil KR. Metabolic dependencies drive species cooccurrence in diverse microbial communities. Proc Natl Acad Sci. 2015;112:6449–54. 29. Chiu H-C, Levy R, Borenstein E. Emergent biosynthetic capacity in simple microbial communities. PLoS Comput Biol. 2014;10:e1003695. 30. Cheng L, Kiewiet MBG, Logtenberg MJ, Groeneveld A, Nauta A, Schols HA, et al. Effects of different human milk oligosaccharides on growth of bifidobacteria in monoculture and co-culture with Faecalibacterium prausnitzii. Front Microbiol. 2020;11:569700. 31. D’hoe K, Vet S, Faust K, Moens F, Falony G, Gonze D, et al. Integrated culturing, modeling and transcriptomics uncovers complex interactions and emergent behavior in a three-species synthetic gut community. eLife. 2018;7:569700. 32. Das P, Ji B, Kovatcheva-Datchary P, Bäckhed F, Nielsen J. In vitro co-cultures of human gut bacterial species as predicted from co-occurrence network analysis. PLoS ONE. 2018;13:e0195161. 33. Weiss AS, Burrichter AG, Raj ACD, von Strempel A, Meng C, Kleigrewe K, et al. In vitro interaction network of a synthetic gut bacterial community. ISME J. 2021;16:1095–109. 34. Wang Y, LaPointe G. Arabinogalactan utilization by Bifidobacterium longum subsp. longum NCC 2705 and Bacteroides caccae ATCC 43185 in monoculture and coculture. Microorganisms. 2020;8:1703. 35. Magnúsdóttir S, Heinken A, Kutt L, Ravcheev DA, Bauer E, Noronha A, et al. Generation of genome-scale metabolic reconstructions for 773 members of the human gut microbiota. Nat Biotechnol. 2017;35:81–9. 36. Heinken A, Khan MT, Paglia G, Rodionov DA, Harmsen HJM, Thiele I. Functional metabolic map of Faecalibacterium prausnitzii, a beneficial human gut microbe. J Bacteriol. 2014;196:3289–302. 37. Heinken A, Thiele I. Systematic prediction of health-relevant human-microbial co-metabolism through a computational framework. Gut Microbes. 2015;6:120–30. 38. Veith N, Solheim M, van Grinsven KWA, Olivier BG, Levering J, Grosseholz R, et al. Using a genome-scale metabolic model of Enterococcus faecalis V583 to assess amino acid uptake and its impact on central metabolism. Appl Environ Microbiol. 2015;81:1622–33. 39. Ottman N, Davids M, Suarez-Diez M, Boeren S, Schaap PJ, dos Santos VAPM, et al. Genome-scale model and omics analysis of metabolic capacities of Akkermansia muciniphila reveal a preferential mucin-degrading lifestyle. Appl Env Microbiol. 2017;83:e01014–7. 40. Marashi S-A, Bockmayr A. Flux coupling analysis of metabolic networks is sensitive to missing reactions. Biosystems. 2011;103:57–66. 41. Segrè D, Vitkup D, Church GM. Analysis of optimality in natural and perturbed metabolic networks. Proc Natl Acad Sci USA. 2002;99:15112–7. 42. Dukovski I, Bajić D, Chacón JM, Quintin M, Vila JCC, Sulheim S, et al. A metabolic modeling platform for the computation of microbial ecosystems in time and space (COMETS). Nat Protoc. 2021;16:5030–82. 43. Mahadevan R, Edwards JS, Doyle FJ. Dynamic flux balance analysis of diauxic growth in Escherichia coli. Biophys J. 2002;83:1331–40. 44. Clasen F, Nunes PM, Bidkhori G, Bah N, Boeing S, Shoaie S, et al. Systematic diet composition swap in a mouse genome-scale metabolic model reveals determinants of obesogenic diet metabolism in liver cancer. iScience. 2023;26:106040. 45. Heirendt L, Arreckx S, Pfau T, Mendoza SN, Richelle A, Heinken A, et al. Creation and analysis of biochemical constraint-based models using the COBRA Toolbox v.3.0. Nat Protoc. 2019;14:639–702. 46. McLaren MR, Willis AD, Callahan BJ. Consistent and correctable bias in metagenomic sequencing experiments. eLife. 2019;8:e46923. Page 24 of 24 Josephetal. BMC Bioinformatics (2024) 25:36 47. Marrec L, Ghenu A-H, Bank C. Challenges and pitfalls of inferring microbial growth rates from lab cultures. Front Ecol Evol. 2023;11:1313500. 48. Schäfer M, Pacheco AR, Künzler R, Bortfeld-Miller M, Field CM, Vayena E, et al. Metabolic interaction models recapitulate leaf microbiota ecology. Science. 2023;381:eadf5121. 49. Lachance J-C, Lloyd CJ, Monk JM, Yang L, Sastry AV, Seif Y, et al. BOFdat: generating biomass objective functions for genome-scale metabolic models from experimental data. PLoS Comput Biol. 2019;15:e1006971. 50. Jansma J, El Aidy S. Understanding the host-microbe interactions using metabolic modeling. Microbiome. 2021;9:16. 51. Magnúsdóttir S, Heinken A, Fleming RMT, Thiele I. Reply to “Challenges in modeling the human gut microbiome.” Nat Biotechnol. 2018;36:686–91. 52. Scott WT Jr, Benito-Vaquerizo S, Zimmermann J, Bajić D, Heinken A, Suarez-Diez M, et al. A structured evaluation of genome-scale constraint-based modeling tools for microbial consortia. PLoS Comput Biol. 2023;19:e1011363. 53. Budinich M, Bourdon J, Larhlimi A, Eveillard D. A multi-objective constraint-based approach for modeling genomescale microbial ecosystems. PLoS ONE. 2017;12:e0171744. 54. Heinken A, Thiele I. Anoxic conditions promote species-specific mutualism between gut microbes in silico. Appl Env Microbiol. 2015;81:4049–61. 55. Kreimer A, Doron-Faigenboim A, Borenstein E, Freilich S. NetCmpt: a network-based tool for calculating the metabolic competition between bacterial species. Bioinformatics. 2012;28:2195–7. 56. Cao Y, Wang Y, Zheng X, Li F, Bo X. RevEcoR: an R package for the reverse ecology analysis of microbiomes. BMC Bioinform. 2016;17:1–6. 57. Goelzer A, Muntel J, Chubukov V, Jules M, Prestel E, Nölker R, et al. Quantitative prediction of genome-wide resource allocation in bacteria. Metab Eng. 2015;32:232–43. 58. Sánchez BJ, Zhang C, Nilsson A, Lahtvee P-J, Kerkhoven EJ, Nielsen J. Improving the phenotype predictions of a yeast genome-scale metabolic model by incorporating enzymatic constraints. Mol Syst Biol. 2017;13:935. 59. Niebel B, Leupold S, Heinemann M. An upper limit on Gibbs energy dissipation governs cellular metabolism. Nat Metab. 2019;1:125–32. 60. Heinken A, Hertel J, Acharya G, Ravcheev DA, Nyga M, Okpala OE, et al. Genome-scale metabolic reconstruction of 7,302 human microorganisms for personalized medicine. Nat Biotechnol. 2023;41:1320–31. 61. Goris J, Konstantinidis KT, Klappenbach JA, Coenye T, Vandamme P, Tiedje JM. DNA–DNA hybridization values and their relationship to whole-genome sequence similarities. Int J Syst Evol Microbiol. 2007;57:81–91. 62. Jain C, Rodriguez-R LM, Phillippy AM, Konstantinidis KT, Aluru S. High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries. Nat Commun. 2018;9:5114. 63. Marinos G, Kaleta C, Waschina S. Defining the nutritional input for genome-scale metabolic models: a roadmap. PLoS ONE. 2020;15:e0236890. 64. Computation and analysis of microbe-microbe metabolic interactions. http:// gibbs. unal. edu. co/ cobra doc/ cobra toolb ox/ tutor ials/ analy sis/ micro beMic robeI ntera ctions/ iframe_ tutor ial_ micro beMic robeI ntera ctions. html. Accessed 5 Sep 2023. 65. Micom documentation. https:// micomdev. github. io/ micom/. Accessed 5 Sep 2023. 66. Growth in a test tube—COMETS documentation. https:// segre lab. github. io/ cometsmanual/ test_ tube/. Accessed 5 Sep 2023. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.