# Community metabolic modeling of host-microbiota interactions through multi-objective optimization

Anna Lambert, Marko Budinich, Maxime Mahé, Damien Eveillard, Samuel Chaffron

► To cite this version:

Anna Lambert, Marko Budinich, Maxime Mahé, Damien Eveillard, Samuel Chaffron. Community metabolic modeling of host-microbiota interactions through multi-objective optimization. iScience, 2024, 27 (6), 10.1016/j.isci.2024.110092. hal-04674843

HAL Id: hal-04674843

https://hal.science/hal-04674843v1

Submitted on 21 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Community metabolic modeling of host-microbiota interactions through multi-objective optimization

<figure><img src="image_2.png" /></figure>

Anna Lambert,
Marko Budinich,
Maxime Mahé,
Samuel Chaffron,
Damien Eveillard

anna.lambert@univ-nantes.fr

(A.L.)

damien.eveillard@univ-nantes.

fr (D.E.)

## Highlights

Computational score to predict metabolic interactions among organisms

Among gut bacteria, probiotics tend to get higher interaction scores with enterocyte

Cross-feeding of choline predicted between LGG and the enterocyte

Minimal gut ecosystem presence predicted to favor enterocyte maintenance
<figure><img src="image_3.png" /></figure>

## Article

# Community metabolic modeling of host-microbiota interactions through multi-objective optimization

Anna Lambert,<sup>1,*</sup> Marko Budinich,<sup>1</sup> Maxime Mahé,<sup>2,3,4</sup> Samuel Chaffron,<sup>1</sup> and Damien Eveillard<sup>1,5,*</sup>

## SUMMARY

The human gut microbiota comprises various microorganisms engaged in intricate interactions among themselves and with the host, affecting its health. While advancements in omics technologies have led to the inference of clear associations between microbiome composition and health conditions, we usually lack a causal and mechanistic understanding of these associations. For modeling mechanisms driving the interactions, we simulated the organism's metabolism using *in silico* genome-scale metabolic models (GEMs). We used multi-objective optimization to predict and explain metabolic interactions among gut microbes and an intestinal epithelial cell. We developed a score integrating model simulation results to predict the type (competition, neutralism, mutualism) and quantify the interaction between several organisms. This framework uncovered a potential cross-feeding for choline, explaining the predicted mutualism between *Lactobacillus rhamnosus* GG and the epithelial cell. Finally, we analyzed a five-organism ecosystem, revealing that a minimal microbiota can favor the epithelial cell's maintenance.

## INTRODUCTION

The human gut microbiota is a complex ecosystem that significantly impacts host health and disease across the life course.<fnref n="1" /> Therefore, clarifying its role and associated mechanisms in shaping health is an essential fundamental research goal. Environmental omics enable a better understanding of this ecosystem by studying its composition and variation associated with host phenotypes.<fnref n="2" /> This has led to identifying bacteria promoting health (e.g., Faecalibacterium Prausnitzii is depleted in patients with inflammatory bowel disease<fnref n="3" />) or disease (e.g., class Betaproteobacteria is enriched in diabetic individuals<fnref n="4" />). However, correlation does not imply causality, and each strain's mode of action is yet to be understood. Disentangling the impact of bacteria within a diverse ecosystem structured by heterogeneous environmental and host factors is challenging.

Gut microbiota metabolism is central to host physiology as its participation in digestion modulates nutrient availability to the host.<fnref n="5" />,⁶ Moreover, the crosstalk among bacteria shapes the ecosystem,⁷ consequently modulating its function.⁸ To investigate mechanistic metabolic interactions between organisms, computational models are useful tools.⁹,¹⁰ Researchers can reconstruct metabolic networks based on an organism's complete genome, biochemical databases, and literature knowledge. Leveraging stoichiometry and thermodynamics information, this network can be formalized into a genome-scale metabolic model (GEM), facilitating the simulation and examination of its metabolic phenotype.¹¹ Indeed, GEMs predict the metabolic phenotype in a given condition and explain this prediction by revealing important nutrients and activated metabolic pathways.¹² Many GEMs are now available,¹³,¹⁴ and automated reconstruction processes for microorganisms have emerged,¹⁵⁻¹⁷ paving the way to analyze various potentially unexplored bacterial metabolic mechanisms. Simulating the metabolism of individual organisms is a well-established and extensively employed approach for expanding fundamental knowledge. For instance, it aids in closing gaps in metabolic knowledge by identifying differences between predictions and experimentations.¹¹ Additionally, this approach contributes to enhancing industrial capacity. An example of this is genetic manipulation, which boosts the production of molecules of interest.¹⁸ To orient the model toward biologically relevant metabolic behavior, its biomass production, representing its growth, is maximized,¹⁹ emphasizing phenotypes where the organism undergoes replication. For a community, this means maximizing the overall biomass produced by all organisms in the ecosystem. While interesting, this method favors the organism with the better yield, possibly outcompeting the others. We expect a more nuanced compromise in biological systems, driven by metabolic trade-offs and resource competition.²⁰ To avoid this bias, the ecosystem biomass is generally weighted based on the relative abundance of each species in a specific condition (e.g., the Microbiome Modeling Toolbox [MMT]²¹), aligning their growth rate to their effective presence in the ecosystem. To go further, Diener et al. introduced

<footnote n="1">¹Nantes Université, École Centrale Nantes, CNRS, LS2N, UMR 6004, 44000 Nantes, France</footnote>

<footnote n="2">²Nantes Université, Inserm, TENS UMR1235, The Enteric Nervous System in Gut and Brain Diseases, IMAD, Nantes, France</footnote>

<footnote n="3">³Division of Pediatric General and Thoracic Surgery, Cincinnati Children's Hospital Medical Center, Cincinnati, OH, USA</footnote>

<footnote n="4">⁴Center for Stem Cell and Organoid Medicine, Cincinnati Children's Hospital Medical Center, Cincinnati, OH, USA</footnote>

<footnote n="5">⁵Lead contact</footnote>

<footnote>*Correspondence: anna.lambert@univ-nantes.fr (A.L.), damien.eveillard@univ-nantes.fr (D.E.)</footnote>

<footnote>https://doi.org/10.1016/j.isci.2024.110092</footnote>
MICOM,²² a framework inferring the growth rate from relative abundance data, refining the weights of the ecosystem biomass. This method calibrates the model to fit observed species proportions inferred from omics data. Parametrizing the model to better simulate a known condition allows for satisfactory predictions with simpler models, but the dependence on environmental data is troublesome. Here, we aimed to achieve non-parametric modeling, striving for accurate predictions through refined model analysis rather than from model fitting.

We used multi-objective linear programming to predict and explore potential interactions between several organisms without restricting ourselves to a unique configuration. This approach allows for the optimization of independent objectives, revealing trade-offs between each organism's biomass as a Pareto front. This is usually applied in bi-objective modeling of microbial ecosystems to infer the type of interaction between bacteria of interest, which is made accessible through the MMT.²¹ Here, we deployed a multi-objective modeling framework to gain further insights into the interaction between gut microbes and small intestinal epithelial cells. We inferred an interaction score predicting the type and level of interaction among organisms within an ecosystem based on the analysis of the Pareto front. This analysis highlighted known probiotics when applied to 331 bacterium-enterocyte ecosystems. Furthermore, we delved into the mechanisms underlying the mutual interaction between *Lactobacillus rhamnosus* GG (LGG) and the enterocyte, uncovering a potential cross-feeding relationship involving choline. Finally, we integrated four bacterial models with the enterocyte model to explore the intricate metabolic interplays within a more complex ecosystem. In this context, our findings reveal that the presence of gut bacteria significantly influences and supports the enterocyte's objective (i.e., maintenance of the cell without replication, see STAR Methods).

# RESULTS

## Scoring metabolic interaction of bacterium-enterocyte ecosystem models using bi-objective optimization

To explore metabolic interactions between the human gut and the microbiota, and quantify them as a score, we built ecosystem models defined as the integration of GEMs of small intestinal epithelial cells (i.e., enterocyte) and one or several bacteria through a pool compartment (Figure 1A). Classically, metabolic modeling allows the maximization of a single objective: the growth of the bacteria or the maintenance of the enterocyte (Figure 1B). The plurality of objectives must be addressed when joining several organisms in an ecosystem. To simulate ecosystem behaviors, we consider that each organism independently maximizes its objective, shaping the competition and cooperation for available nutrients.²³ It results in a problem one can solve using multi-objective linear programming. The set of all optimal solutions for the ecosystem is a Pareto front. It reflects the trade-offs between all involved biological objectives. Herein, we advocate that the shape of the Pareto front describes the nature and level of metabolic interaction between the organisms and can be summarized as a score.

As a first approximation, we worked with simple ecosystems consisting of a single bacterium and the enterocyte. In this context, the Pareto front comprises two dimensions: one for each objective. Maximal growths of each organism alone are used to normalize each axis as they standardize growths under interaction regimes (i.e., studying the interaction between organisms without an influence of the original growth). This transformation is essential to make our interaction score comparable and interpretable between ecosystems built using different bacteria. When both objectives negatively or positively impact the other, both organisms are in competition or mutualism, respectively. When neither objective value is affected by the other, it is considered neutralism (see Figure 1C for illustration of three representative Pareto front shapes). We propose an ecosystem interaction score (S) based on the normalized Pareto front's area under the curve (AUC). The AUC of the "non-interaction front," defined by the neutral interaction, is subtracted from the AUC of the normalized Pareto front. It results in a positive score for mutualism, a null score for neutrality, and a negative score for competition (see STAR Methods for details).

Analyzing the ecosystem interaction dynamics via the organisms objective trade-offs, or Pareto front, represents a progression from methods that compare growth and co-growth rates, such as the one outlined in Magnúsdóttir et al.¹⁷ (Figure S1). Although the MMT²¹ incorporates Pareto front analysis, we advance this technique by quantifying the ecosystem interaction potential.

## Interaction score for 331 gut bacterium-enterocyte ecosystems under three different dietary conditions

331 strain-level gut bacteria models from the European Molecular Biology Laboratory (EMBL) GEMs¹⁶ were selected as they were described as gut microbes based on the Virtual Metabolic Human (VMH) database.¹⁴ Their pairwise interaction score with the enterocyte was computed under three nutritional conditions (Figure 1A, Table S1 for complete results). The first condition, an unconstrained diet, represents a synthetic environment where all modeled nutrients are unlimited in the lumen. The other two conditions, the Western diet (WD) and the protein diet (PD), are biologically representative diets. The WD is high in fat and simple sugars but low in fibers, while the PD is rich in protein but poorer in fats and simple sugars than the WD²⁴ (see STAR Methods). An increased number of non-dominated points on the Pareto front reveals more factors affecting the interaction, characterizing its complexity.

Overall, ecosystems had an interaction score close to zero under the unconstrained diet but exhibited negative and positive scores under the WD and PD (Figure 2A). This differential interaction scoring in constrained media was expected, illustrating either the mutual reliance on a common metabolite leading to competition or the adaptation through cross-feeding, usually favored in nutrient-limited environments.²⁵ In other words, these findings predict how metabolic interactions are likely driven by restricted nutritional environments. Additionally, this result serves as a point of comparison for the impact of changing nutrients availability in the medium.

The sensitivity of the interaction score to diets is demonstrated by observed differences in scores between PD and WD for a given ecosystem. This highlights how nutrient availability influences ecosystem interactions, which is consistent with the known impact of diet on the gut microbiota composition and metabolism.²⁶ Furthermore, ecosystems with positive scores generally had an increased score under
<figure><img src="image_4.png" /></figure>

A

<figure><img src="image_5.png" /></figure>

B

<figure><img src="image_6.png" /></figure>

C

Ecosystem interaction score (S) = Pareto AUC - non-interaction AUC

<figure><img src="image_7.png" /><figcaption>**Figure 1. Schematic rationale of the genome-scale modeling and its application in the metabolic interaction between two systems within an ecosystem**</figcaption></figure>

(A) Representation of an ecosystem model built from a bacterium (blue) and the enterocyte's (red) GEMs joined together through a lumen compartment (gray). Arrows represent metabolic reactions. Black arrows are transport reactions, enabling the transit of metabolites between compartments. Three different diets (unconstrained, western diet [WD] and protein diet [PD]) were imposed on the ecosystem by defining corresponding dietary constraints (see STAR Methods).

(B) Solving a GEM implies identifying its solution space and optimal solutions when maximizing its objective function, often resumed by its growth rate. The solution space of a bacterium (left) described here in three dimensions illustrates how the optimal solutions for the bacterium were identified by maximizing its growth. Similarly, the optimal solutions for the enterocyte are identified by maximizing its cell maintenance (cellular membrane maintenance, proteins and energy production, no replication). Both objectives must be accounted for in an ecosystem model, which is possible using multi-objective optimization.

(C) Schematic representation of three different Pareto fronts and how it can be used to define the ecosystem interaction score (S). The shape and area of the bi-objective solution space can be used to define competition (left panel), neutralism (center panel, defining the non-interaction front), and mutualism (right panel). Dotted lines link the added "original growth points" and the strict Pareto front.

PD compared to WD, while those with negative scores showed a decrease in score under WD compared to PD (Figure 2B). In other words, the PD appeared to enhance positive interactions between bacteria and the enterocyte, while the WD is predicted to favor negative interactions. Interestingly, there is no case of an ecosystem yielding a score of opposite sign between the two diets, implying that the diet influence was insufficient to induce a shift in ecosystem predicted interactions from competition to mutualism or vice versa.

**Further characterization of the Pareto front into "interaction types"**

The shapes of the Pareto fronts are more diverse than depicted in Figure 1C. They can inform us about the nature of biotic interactions. To leverage this information and enrich the metabolic interaction score, we defined discrete categories as "interaction type" (Figure S2, see STAR Methods). Indeed, a score can be positive but only favoring one organism, as represented by the interaction types "Favors bacteria" and "Favors enterocyte." In some cases, the presence of the other can favor both organisms, but this advantage disappears as the other organism gets to maximize its objective, as observed in "Limited mutualism." "Competition," "Neutralism," and "Mutualism" interaction types correspond to the shapes described in Figure 1C.

To determine how S varied based on diet in each interaction type, we distributed the ratio calculated previously ($S_{PD}/S_{WD}$) by interaction types (Figure 2D). Most positive interactions, specifically those favoring the enterocyte (favors enterocyte, limited mutualism), were increased in PD. However, interaction scores from the mutualism interaction type were higher in WD. Interestingly, the increase of S in WD for the mutualistic ecosystem is not explained by reaching a higher objective value but because the interaction partially compensates for the lower objective value of each organism in WD compared to PD (Table S2). More than a beneficial effect of WD, this observation illustrates the potential of metabolic interactions to improve fitness in restricted conditions, resulting in likely increased metabolic cross-feedings in nutrients-limited environments.²⁵ While WD enhanced some interactions favoring bacteria, it mostly favored competitive interactions.
<figure><img src="image_8.png" /><figcaption>**Figure 2.** Prediction of metabolic interaction scores and types for 331 enterocyte-bacterium ecosystems</figcaption></figure>

(A) Interaction score (S) 2D distribution and density curve among 331 ecosystems subject to three diets (unconstrained, WD, PD), with a zoom on negative values. Here, the number of points that constitute the Pareto front approximates the topological complexity of organismal interactions. Thus, higher number of points, as illustrated in Figure S2, will be associated with higher number of interaction scenarios or type of interactions.

(B) Evaluation of the effect of the diet on the score through a $S_{WD}^{SD}$ distribution colored by the sign of the score.

(C) Assessment of the predictive potential of S on the ecological type (probiotics $n = 18$, pathogens $n = 46$, commensals $n = 123$, uncharacterized $n = 144$) colored by diet (WD or PD). A Mann-Whitney U-test was used to determine statistical significance. (*: p-value < 0.05)

(D) Evaluation of the impact of diet on the score in interaction types using $S_{WD}^{SD}$ ratio distribution.

(E) Proportions of each Pareto inferred interaction type and exploration of how ecological types place themselves on the interaction types. Specification of bacterial strain predicted as mutualist with the enterocyte.

See also Figures S1 and S2.

## Integration of metabolic interaction score and type to predict "host-beneficial" microbes

The VMH database¹⁴ provided the ecological type of the 331 bacterial GEMs. Microorganisms were categorized into three ecological types: probiotic, pathogen, commensal. Any microorganism that had no defined ecological type was categorized as uncharacterized. Probiotics had a significantly higher score compared to pathogens in PD (Mann-Whitney, *p* = 0.033) and commensals in WD (Mann-Whitney, *p* = 0.040) (Figure 2C). Overall, we observed a tendency for higher scores in known probiotics (Figure 2C). This result raises interest in some uncharacterized bacteria displaying a high interaction score with the enterocyte as *Cetobacterium somerae* ATCC BAA-474 ($S_{PD}$ = 23.64, $S_{WD}$ = 12.10), *Klebsiella aerogenes* KCTC 2190 ($S_{PD}$ = 16.92, $S_{WD}$ = 8.70), or *Morganella morganii* subsp *morganii* KT ($S_{PD}$ = 11.42, $S_{WD}$ = 5.80).

Next, we explored the ecological type distribution among interaction types (Figure 2E). A third of the simulated bacteria (33%), encompassing all ecological types, were predicted to have a neutral interaction with the enterocyte. Many bacteria favored the enterocyte (22%), largely dominated by commensals and uncharacterized bacteria. Ecosystems where the bacteria's growth was favored (20%) contained more pathogens and probiotics than the ones favoring the enterocyte. 17% of the bacteria were engaged in a competitive interaction with the host. Limited and high mutualisms were identified as the least common interaction types, accounting for only 5% and 2% of the analyzed bacteria, respectively. Probiotics corresponded to half of the mutualistic bacteria ($N$ = 6) with two lactic acid bacteria and a short-chain fatty acid (SCFA) producer (LGG, *Lactobacillus paracasei* ATCC 334, and *Bifidobacterium scardovii* JCM 12489). Additionally, unexpected mutualistic interactions were observed for *Corynobacterium kroppenstedtii* DSM 44385 (categorized as a pathogen), *Barnesiella*
<figure><img src="image_9.png" /></figure>

*intestinihominis* YIT 11860 (classified as a commensal), and *Peptoniphilus indolicus* ATCC 29427 (tagged as uncharacterized) with the enterocyte.

## Metabolic exchanges driving mutualism between the enterocyte and LGG reveal a potential cross-feeding of choline

LGG is a very well-studied probiotic bacteria. Among other properties, it shows good adherence to the intestinal epithelial layer and supports the survival of intestinal epithelial cells.²⁷ Our computational analysis predicts LGG to engage in a mutualistic relationship with the host, as reflected by its interaction scores of 8.68 in WD and 6.60 in PD. Their Pareto front formed a spike from which the peak allowed the highest biomass production for both organisms and, consequently, for the ecosystem (See Figure S3). To unravel the underlying mechanisms driving this mutualistic interaction by identifying essential reactions to reach this optimum, we sampled 100,000 solutions on the Pareto front in WD. We focused on exchange reactions that exhibited a strong correlation (corr ≥0.95, Spearman) with the overall biomass production of the ecosystem (i.e., summed objective values, here, the peak of the Pareto front), as depicted in Figure 3A.

The exchange reactions identified through this analysis enabled both organisms to reach higher objective values when interacting, helping to predict valuable nutrients and cross-feeding pathways. The intricate metabolic interplay within the ecosystem involves the shared utilization of various metabolites, contributing to the overall ecosystem functioning. For the ecosystem to reach optimal biomass, we predicted the uptake of many amino acids and dipeptides to be essential, mainly for the enterocyte. However, the digestion of leucine-leucine and glycine-phenylalanine dipeptides by the enterocyte's enzymes provided the bacteria with leucine and phenylalanine.⁷ Oligo-elements and vitamins were more important for LGG, but the enterocyte consumed biotin and riboflavin (7) and released them in the blood. Leucine, phenylalanine, and riboflavin exchanges between the enterocyte and the bacteria were further supported by their identification as exchanged metabolites (See STAR Methods, Figure S4).

The secretion of choline by LGG was found to be strongly correlated with its uptake by the enterocyte (corr = 0.96, Figure S4), and both these transports are associated with a higher biomass for the ecosystem. Specifically, as illustrated in Figure 3B, our model predicted that LGG uptakes choline sulfate (1), a compound found in human food such as plants, algae, and numerous fungi such as *Aspergillus oryzae*, a ferment used for sake, miso, and soy sauce.²⁸,²⁹ LGG hydrolyzes choline sulfate (2) in sulfate ions (SO₄²⁻), which are useful for LGG's growth, and choline, secreted in the intestinal lumen. LGG reduces the excess sulfate to H₂S (3) and excretes it in the lumen. The choline, now available in the lumen, is uptaken by the enterocyte (4). In the enterocyte, the absorbed choline transforms into CDP-choline (5) which, along with diacylglycerol, further converts to phosphatidylcholine. Phosphatidylcholine is the principal component of cell membranes and participates in lipid metabolism,³⁰ therefore regarded as an essential brick for cell maintenance.

Choline is an essential nutrient and source of methyl. However, less than half of the tested populations (adult men and women) reached the recommended intake.³¹ Therefore, it is consistent with the literature to expect its supplementation, through cross-feeding with LGG, to favor the enterocyte's maintenance. Phosphatidylethanolamine (PE), another abundant phospholipid,³⁰ was available in excess in the environment when the enterocyte was simulated on its own. However, since phosphatidylcholine was limited in availability, the enterocyte's maintenance and PE intake were constrained. In the community with LGG, as the choline availability increased, the enterocyte also exhibited an elevated absorption of PE (6), conjointly leading to improved maintenance for the enterocyte.

## Modeled minimal gut microbiome metabolism greatly favors the enterocyte

An ecosystem is formed of more than two organisms, and, to move toward modeling realistic gut-microbiota interactions, more organisms must be included in the modeled community. The bacterial strains from Shetty et al.'s minimal microbiome³² were modeled using CarveMe,¹⁶ and four were selected because of their distinct metabolic potential and predicted phenotypes: *Akkermansia muciniphila* ATCC BAA-835, *Bacteroides xylanisolvens* HMP 2_1_22, *Faecalibacterium prausnitzii* A2-165, and *Ruminococcus bromii* ATCC 27255. It should be noted that, while four bacteria were chosen for this analysis, it is computationally applicable to up to 10 bacteria and the enterocyte (Figure S5). They were integrated with the enterocyte in an ecosystem model, which was analyzed using multi-objective linear programming. Here, as the Pareto front is in five-dimensional space, its description was reduced to its extreme points for analytical purposes, constituting the extreme solutions of the ecosystem.

Solutions from the Pareto front represent extreme community phenotypes where all available nutrient usage is optimized. To envision a more realistic (i.e., suboptimal) set of community phenotypes, three thousands additional solutions were randomly sampled within the entire solution space (i.e., solutions within the volume embedded by the Pareto front, such as the gray surface in Figure 1 when applied in 2 dimensions), combining optimal and random solutions results in an inclusive set of potential phenotypes for the ecosystem. Here, solutions are described by a 5-dimensional vector, each value being an organism's objective. To visualize how each organism's objective impacts the rest of the ecosystem, we performed a principal-component analysis (PCA) on the objective values (Figure 4A). The PCA revealed that PC1, which explains 30% of the variance, strongly aligns with an increase in ecosystem biomass and is coherently bordered by extreme solutions. Notably, all organisms contributed to the increase in ecosystem biomass. PC2 (21% of explained variance) refined this interpretation and distinguished distinct groups within the ecosystem: *Bacteroides xylanisolvens* and the enterocyte exhibited a similar trend, while *Akkermansia muciniphila* showed an opposing direction. *Ruminococcus bromii* and *Faecalibacterium prausnitzii* fell between these groups.

Inferring the interaction score for the ecosystem (see STAR Methods) is informative of its overall interaction. Here, an ecosystem interaction score of 1.70 implies that the modeled organisms are dominated by mutualist interaction compared to competition. However, this lacks pre-cisions about the role of each actor in the interaction. As a strategy to highlight the impact of a given organism on the ecosystem's interaction potential, we calculated the score of smaller ecosystems after the removal of each organism. Therefore, we measured the score for five
A

Uptaken

<figure><img src="image_10.png" /></figure>

Secreted

<figure><img src="image_11.png" /></figure>

B

<figure><img src="image_12.png" /><figcaption>**Figure 3. Exploring the Pareto front of the LGG-enterocyte interaction reveals a possible cross-feeding involving choline**</figcaption></figure>

(A) We sampled the Pareto front of the interaction between LGG (blue) and the enterocyte (red). We inferred the correlations between transport reaction and ecosystem biomass to identify metabolic cross-feedings driving the interaction. Metabolites for which the exchanges are correlated with the ecosystem biomass, therefore driving the interaction, are described in boxplots by the distribution of their possible fluxes among the sampling. A metabolite is categorized as uptaken (or secreted) when its absorption (or excretion) by the relevant organism is associated with ecosystem biomass (see STAR Methods). The color of the boxplot refers to the relevant organism (blue for LGG, red for the enterocyte). White boxplots refer to metabolites transported into the blood from the enterocyte. Metabolites involved in the metabolism depicted in (B) are in bold. (The box of a boxplot extend from quartile one to quartile three of the data. The median is marked by a line. The whiskers extend from the box to the farthest data point lying within 1.5x the inter-quartile range from the box).

(B) Schematic representation of the predicted cross-feeding of choline. The arrows represent metabolic reactions. Dashed arrows represent the exchange reactions (i.e., the consumption of metabolites from the diet). Blue and red arrows illustrate the use of a metabolite for the objective of the LGG and the host, respectively. The metabolic pathway colored in green is the predicted cross-feeding of choline. All reactions pictured in this figure are obligatory (See STAR Methods) in the most mutualistic solution of the Pareto front.

See also Figures S3 and S4.
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Figure 4. Multi-objective community metabolic modeling of five organisms in an ecosystem (enterocyte, *Akkermansia muciniphila*, *Bacteroides xylanisolvens*, *Faecalibacterium prausnitzii*, and *Ruminococcus bromii*)</figcaption></figure>

(A) Principal-component analysis of solutions defined by the objective value of each organism. The point's size is proportional to the sum of this solution's objective values. Gray points result from a random sampling in the solution space of the ecosystem. Black and colored points are extreme solutions defining the Pareto front. Colored points are a selection of extreme solutions from the Pareto front where every organism (objective) has a non-null objective value, representing cohabiting metabolic phenotypes for the ecosystem where each organism can grow or maintain itself.

(B) Upset diagram plot of the impact of the absence of each organism on the ecosystem interaction score. Scaling the interaction score from bi-objective to four-dimensional ecosystems is required to calculate hypervolumes instead of the area under the curve (compared to Figure 1) but follows the same principle (see STAR Methods).

(C) Visualization of the four cohabiting metabolic phenotypes for the ecosystem (solutions where every organism had a non-null objective value). The axes are objective values for the corresponding organism, normalized by its highest observed objective value (in the ecosystem or alone).
See also Figure S5.
reduced ecosystems (removing a different organism each time) to better understand the role of each organism in the community's interplay (Figure 4B). Removing the enterocyte or *B. xylanisolvens* from the ecosystem yielded negative scores (−0.44 and −0.18). This suggests that they are involved in a beneficial relationship, raising the overall score of the ecosystem. The score remained positive when removing *F. prausnitzii* or *R. bromii* (1.17 or 1.26). The highest interaction score was observed when *A. muciniphila* was absent from the ecosystem (2.16). Since one cannot mathematically compare these four-dimensional scores to a five-dimensional score when all organisms are present in the ecosystem, we cannot conclude whether *Akkermansia muciniphila* exhibits a weaker but still beneficial interaction or if it decreases the overall score of the ecosystem. Still, it is the bacteria with the least favorable effect on the score. The impact of each organism on the ecosystem score is coherent with their disposition on the PCA (Figure 4A), concluding that PC2 is correlated with the interaction score.

Four extreme solutions were selected where no organism had a null objective value (Figure 4C). The maintenance of the enterocyte was highly favored by the presence of bacteria, with an improved objective value in all selected solutions. However, it is noteworthy that the growth of the bacteria within the ecosystem is reduced compared to their maximum growth potential when in isolation. The objective values of the enterocyte and the *B. xylanisolvens* are correlated across solutions. Besides, *R. bromii*, *F. prausnitzii*, and particularly *A. muciniphila* seemed to follow an opposite trend. In this situation, *B. xylanisolvens* seems to be involved in a mechanism essential to the enterocyte's improved maintenance. *B. xylanisolvens* is a complex polysaccharide degrader,³³,³⁴ making carbohydrates and folate available to the enterocyte, which may explain this positive interaction. However, when reaching optimum, objectives of the enterocyte and *B. xylanisolvens* were antagonists to the growth of *A. muciniphila*. *A. muciniphila* is a next-generation probiotic associated with intestinal and systemic health, thought to preserve gut barrier function in the intestinal tract.³⁵ The known metabolic action of *A. muciniphila* on the gut epithelium is through SCFAs production.³⁶ However, SCFA's impact is prevalent in the colon³⁷ and is not modeled in the enterocyte used in this study.

DISCUSSION

Deciphering the interactions between epithelial cells and gut microbes is a crucial step toward a better understanding of human health. This motivates their exploration through the use of GEMs to uncover potential metabolic interplays. This study uses multi-objective modeling to capture the trade-offs between the host's and gut bacteria's objectives and predict putative metabolic interaction mechanisms. We summarized this complex information by a generic score that indicates overall collaborative, competitive, or neutral interactions between organisms, as well as the quantification of this potential interaction. Known probiotics¹⁴ were associated with a higher interaction score with the enterocyte, raising interest in the score's potential as a screening tool to assess interactions between organisms. Interaction types were computed from the Pareto front shape (i.e., points that build the Pareto front), further characterizing the metabolic interaction. Indeed, the topological complexity of the Pareto front approximates the diversity of metabolic phenotypes within the frame of interactions.³⁸

The interaction score is dependent on the diet imposed on the ecosystem. Notably, the PD favored more mutualistic interactions, while the WD increased competition. However, for sustaining bacterial growth, the nutritional constraints applied to model the diets were under-constrained. Specifically, for all diets, when nutrient concentrations are known, they are adequately restricted. However, when nutrient concentrations are unspecified, we must make them freely available in the lumen. Indeed, considering them as absent yields to an absence of growth, which confirms the difficulty of building a set of constraints that reflect nutritional constraints.³⁹,⁴⁰

Given that under-constrained ecosystems still highlight the influence of the diet on metabolic interactions, we anticipate that further research, with more targeted subjects and strictly defined diets, will yield predictions of higher accuracy. Overall, the present predictions are dependent on the models' quality. We, therefore, expect that, as scientific effort lessens classic modeling limitations (See Limitations of the study), the exactitude of the predictions will improve.

LGG, a known probiotic,²⁷ has been identified in silico as such through a high interaction score and a mutualistic interaction type. This predicted interaction was explored to identify the metabolic mechanisms behind it. This study identified nutrients driving the beneficial interaction between LGG and the host, highlighting the importance of vitamins, amino acids, and oligo-elements, already recognized as central cross-feeding metabolites.⁷ Moreover, cross-feeding was predicted, implying the desulfatase of choline sulfate by LGG, making choline available for the enterocyte while using sulfate ions to grow. As shown before, choline metabolism into phosphatidylcholine participates in cell membrane synthesis and lipid metabolism.³⁰

Choline can be converted to trimethylamine (TMA) by the colonic microflora, which, once absorbed by the colonocytes and transported to the liver, can be oxidized in trimethylamine N-oxide (TMAO).⁴¹ TMAO is associated with various diseases, such as cardiovascular disease, when found in high concentrations.⁴² However, due to the use of overapproximations for nutritional constraints, the predicted mechanisms are also prone to be over- or under-estimated. In this context, as choline sulfate availability is modeled in excess, we anticipate that the exchange of choline would occur in smaller quantities than initially predicted, thereby preventing the accumulation of high concentrations of TMAO in the liver. Presently, in vitro experiments are essential to confirm this cross-feeding.

The score and analysis developed in this study allow an automatic exploration of ecosystem metabolic interplays. Here, we demonstrated their scalability to an ecosystem of five organisms, paving the way toward the description of more intricate ecosystem interactions. However, in practice, such modeling remains limited by available computational resources. In its current form, our largest ecosystem scales up to eleven organisms, including the enterocyte. Beyond the sole computing power, expanding the model size requires complexification of result representations. Extending ecosystem size calls for further developments and optimizations for giving access to more representative ecosystem models.

The application of the analysis described in this paper to communities of interest, such as a consortium of gut bacteria associated with a health condition, could participate in identifying pathogenic or therapeutic pathways. Similarly, the health impact of dietary products with known
<figure><img src="image_15.png" /></figure>

bacterial composition, like cheese or yogurt,⁴³ could be further explored. Moreover, assessing which type of epithelial cell interacts the most with a species could enrich our knowledge of how different gut locations build different communities.⁴⁴ As an illustration, we expect *Akkermansia muciniphila* to reach a higher interaction score with a colonocyte than the small intestine epithelial cell used in this work. While this framework was demonstrated here on metabolic interactions between the human gut epithelium and the microbiota, its application is generic enough to be relevant for various ecosystems and contexts. Tumors are accompanied by bacteria, which promote or suppress cancer based on the situation.⁴⁵ Applied in this tumor-microenvironment ecosystem, a positive score would highlight pathogenic interaction, opening the way for adapted treatments, like targeted antibiotics. Finally, this can be applied to non-health-related fields, such as the study of ocean and soil ecosystems. Specifically, it can allow the study of uncultured strains through the modeling of metagenome-assembled genomes (MAGs).

In this work, we converge with previous results of Heinken and Thiele,²⁴ using multi-objective optimization to provide a mathematical proof of symbiosis. Indeed, we demonstrate here that two organisms in nutritional co-limitation can attain better fitness than alone. This observation justifies the development of tight multi-species communities in various ecological contexts, as their cohabitation leads to better survival. Additionally, we show that the diet modulates the importance of this mutualism. Compared to an open environment such as the ocean, the gut is an enclosed habitat where diet is the principal nutritional input. On the other hand, the microbiota regulates food intake (in terms of quantity and quality) through the gut-brain axis.⁴⁶ Therefore, our nutrition forges our microbiota's interactions and fitness, shaping its composition,⁴⁷ while the microbiota influences our nutrition.⁴⁶ Controlling the microbiota's early assembly and ecology through an adapted diet is key to favoring a self-maintaining healthy composition.⁴⁸ This encourages pursuing research on the gut ecosystem metabolism and how microbiota and nutrition interact and impact human health.

## Limitations of the study

As previously stated in the discussion with nutritional constraints, models are approximative representations and present limitations. The models used in this study were automatically reconstructed from a genome using the CarveMe pipeline.¹⁶ Due to their dependence on genome completion and lack of curation based on the literature on the specific strain, these models are more susceptible to diverging from chemical knowledge, potentially leading to the erroneous representation of metabolic pathways.

Here, the cross-feeding predicted between LGG and the enterocyte relies on the ability of LGG to metabolize choline sulfate into choline and sulfate. Choline sulfatase is a protein synthesized from the gene BetC. This gene was found in many bacteria but not in LGG.⁴⁹ This is a reaction added through gap filling, meaning that it was necessary for the growth of LGG *in silico* although not found in its genome. From this information, we can hypothesize that LGG is able to perform an analogous reaction utilizing an alternate enzyme. Another possibility is that another pathway is necessary for LGG to grow, but a lack of information resulted in the addition of the choline desulfatase pathway to compensate for its absence. Overall, this reveals a gap in knowledge in the metabolism of LGG, requiring further investigation.

Additionally, the prediction of this method highly depends on the defined objectives. Indeed, the model would benefit from including more than maintenance in the enterocyte's objective function, as it is a differentiated cell with a specific role for the human body. Consequently, with the ongoing evolution of objective definitions and comprehensive curation, a corresponding improvement is expected in the score's accuracy.

## STAR★METHODS

Detailed methods are provided in the online version of this paper and include the following:

* **KEY RESOURCES TABLE**

* **RESOURCE AVAILABILITY**

* Lead contact

* Materials availability

* Data and code availability

* **METHOD DETAILS**

* Genome-scale metabolic models

* Ecosystem model: Pool compartment, diet

* Multi-objective linear programming: Pareto front

* Linear problem formalisms

* Interaction score

* Interaction types

* Sampling of the Pareto front

* Exchanged metabolites and obligatory reaction

* Building a five-dimension ecosystem

* Five-dimensional principal component analysis (PCA)

* **QUANTIFICATION AND STATISTICAL ANALYSIS**

## SUPPLEMENTAL INFORMATION

Supplemental information can be found online at https://doi.org/10.1016/j.isci.2024.110092.
## ACKNOWLEDGMENTS

A.L. benefits a doctoral fellowship managed by Bba Milk Valley, dairy industrial association. We thank Regional Councils of Bretagne (grant no. 19008213) and Pays de la Loire (grant no. 2019-013227) for their financial support through the interregional project PROLIFIC. The bioinformatics core facility of Nantes (BiRD, Biogenouest), Nantes Université, France, provided computational support.

## AUTHOR CONTRIBUTIONS

Conceptualization, A.L., D.E., S.C., and M.M.; methodology, A.L., D.E., and S.C.; software, A.L. and M.B.; formal analysis, A.L., investigation, A.L.; writing – original draft, A.L.; writing – review and editing, A.L., D.E., S.C., M.M., and M.B.; visualization, A.L.; supervision, D.E., S.C., and M.M.

## DECLARATION OF INTERESTS

The authors declare no competing interests.

Received: October 30, 2023

Revised: March 12, 2024

Accepted: May 21, 2024

Published: May 23, 2024

## REFERENCES

1. Thursby, E., and Juge, N. (2017). Introduction to the human gut microbiota. Biochem. J. 474, 1823-1836. https://doi.org/10.1042/BCJ20160510.

2. Rothschild, D., Weissbrod, O., Barkan, E., Kurilshikov, A., Korem, T., Zeevi, D., Costea, P.I., Godneva, A., Kalka, I.N., Bar, N., et al. (2018). Environment dominates over host genetics in shaping human gut microbiota. Nature 555, 210-215. https://doi.org/10.1038/nature25973.

3. Cao, Y., Shen, J., and Ran, Z.H. (2014). Association between Faecalibacterium prausnitzii Reduction and Inflammatory Bowel Disease: A Meta-Analysis and Systematic Review of the Literature. Gastroenterol. Res. Pract. 2014, e872725. https://doi.org/10.1155/2014/872725.

4. Larsen, N., Vogensen, F.K., van den Berg, F.W.J., Nielsen, D.S., Andreasen, A.S., Pedersen, B.K., Al-Soud, W.A., Sørensen, S.J., Hansen, L.H., and Jakobsen, M. (2010). Gut Microbiota in Human Adults with Type 2 Diabetes Differs from Non-Diabetic Adults. PLoS One 5, e9085. https://doi.org/10.1371/journal.pone.0009085.

5. Flint, H.J., Scott, K.P., Duncan, S.H., Louis, P., and Forano, E. (2012). Microbial degradation of complex carbohydrates in the gut. Gut Microb. 3, 289-306. https://doi.org/10.4161/gmic.19897.

6. Jandhyala, S.M., Talukdar, R., Subramanyam, C., Vuyuru, H., Sasikala, M., and Nageshwar Reddy, D. (2015). Role of the normal gut microbiota. World J. Gastroenterol. 21, 8787-8803. https://doi.org/10.3748/wjg.v21.i29.8787.

7. Culp, E.J., and Goodman, A.L. (2023). Cross-feeding in the gut microbiome: Ecology and mechanisms. Cell Host Microbe 31, 485-499. https://doi.org/10.1016/j.chom.2023.03.016.

8. Khoruts, A., Dicksved, J., Jansson, J.K., and Sadowsky, M.J. (2010). Changes in the Composition of the Human Fecal Microbiome After Bacteriotherapy for Recurrent Clostridium difficile-associated Diarrhea. J. Clin. Gastroenterol. 44, 354-360. https://doi.org/10.1097/MCG.0b013e3181c87e02.

9. Oberhardt, M.A., Palsson, B.Ø., and Papin, J.A. (2009). Applications of genome-scale metabolic reconstructions. Mol. Syst. Biol. 5, 320. https://doi.org/10.1038/msb.2009.77.

10. Thiele, I., Heinken, A., and Fleming, R.M.T. (2013). A systems biology approach to studying the role of microbes in human health. Curr. Opin. Biotechnol. 24, 4-12. https://doi.org/10.1016/j.copbio.2012.10.001.

11. Bordbar, A., Monk, J.M., King, Z.A., and Palsson, B.O. (2014). Constraint-based models predict metabolic and associated cellular functions. Nat. Rev. Genet. 15, 107-120. https://doi.org/10.1038/nrg3643.

12. Passi, A., Tibocha-Bonilla, J.D., Kumar, M., Tec-Campos, D., Zengler, K., and Zuniga, C. (2021). Genome-Scale Metabolic Modeling Enables In-Depth Understanding of Big Data. Metabolites 12, 14. https://doi.org/10.3390/metabo12010014.

13. King, Z.A., Lu, J., Dräger, A., Miller, P., Federowicz, S., Lerman, J.A., Ebrahim, A., Palsson, B.O., and Lewis, N.E. (2016). BiGG Models: A platform for integrating, standardizing and sharing genome-scale models. Nucleic Acids Res. 44, D515-D522. https://doi.org/10.1093/nar/gkv1049.

14. Noronha, A., Modamio, J., Jarosz, Y., Guerard, E., Sompairac, N., Preciat, G., Daniëlsdóttir, A.D., Krecke, M., Merten, D., Haraldsdóttir, H.S., et al. (2019). The Virtual Metabolic Human database: integrating human and gut microbiome metabolism with nutrition and disease. Nucleic Acids Res. 47, D614-D624. https://doi.org/10.1093/nar/gky992.

15. Büchel, F., Rodriguez, N., Swainston, N., Wrzodek, C., Czauderna, T., Keller, R., Mittag, F., Schubert, M., Glont, M., Golebiewski, M., et al. (2013). Path2Models: large-scale generation of computational models from biochemical pathway maps. BMC Syst. Biol. 7, 116. https://doi.org/10.1186/1752-0509-7-116.

16. Machado, D., Andrejev, S., Tramontano, M., and Patil, K.R. (2018). Fast automated reconstruction of genome-scale metabolic models for microbial species and

communities. Nucleic Acids Res. 46, 7542-7553. https://doi.org/10.1093/nar/gky537.

17. Magnúsdóttir, S., Heinken, A., Kutt, L., Ravcheev, D.A., Bauer, E., Noronha, A., Greenhalgh, K., Jäger, C., Baginska, J., Wilmes, P., et al. (2017). Generation of genome-scale metabolic reconstructions for 773 members of the human gut microbiota. Nat. Biotechnol. 35, 81-89. https://doi.org/10.1038/nbt.3703.

18. Gu, C., Kim, G.B., Kim, W.J., Kim, H.U., and Lee, S.Y. (2019). Current status and applications of genome-scale metabolic models. Genome Biol. 20, 121. https://doi.org/10.1186/s13059-019-1730-3.

19. Feist, A.M., and Palsson, B.O. (2010). The biomass objective function. Curr. Opin. Microbiol. 13, 344-349. https://doi.org/10.1016/j.mib.2010.03.003.

20. Posfai, A., Taillefumier, T., and Wingreen, N.S. (2017). Metabolic Trade-Offs Promote Diversity in a Model Ecosystem. Phys. Rev. Lett. 118, 028103. https://doi.org/10.1103/PhysRevLett.118.028103.

21. Baldini, F., Heinken, A., Heirendt, L., Magnusdóttir, S., Fleming, R.M.T., and Thiele, I. (2019). The Microbiome Modeling Toolbox: from microbial interactions to personalized microbial communities. Bioinformatics 35, 2332-2334. https://doi.org/10.1093/bioinformatics/bty941.

22. Diener, C., Gibbons, S.M., and Resendis-Antonio, O. (2020). MICOM: Metagenome-Scale Modeling To Infer Metabolic Interactions in the Gut Microbiota. mSystems 5, e00606-19. https://doi.org/10.1128/mSystems.00606-19.

23. Budinich, M., Bourdon, J., Larhlimi, A., and Eveillard, D. (2017). A multi-objective constraint-based approach for modeling genome-scale microbial ecosystems. PLoS One 12, e0171744. https://doi.org/10.1371/journal.pone.0171744.

24. Heinken, A., and Thiele, I. (2015). Anoxic Conditions Promote Species-Specific Mutualism between Gut Microbes In Silico. Appl. Environ. Microbiol. 81, 4049-4061. https://doi.org/10.1128/AEM.00101-15.

25. Mataigne, V., Vannier, N., Vandenkoornhuyse, P., and Hacquard, S.
<figure><img src="image_16.png" /></figure>

(2021). *Microbial Systems Ecology to Understand Cross-Feeding in Microbiomes*. *Front. Microbiol.* 12, 780469.

26. Scott, K.P., Gratz, S.W., Sheridan, P.O., Flint, H.J., and Duncan, S.H. (2013). The influence of diet on the gut microbiota. *Pharmacol. Res.* 69, 52–60. https://doi.org/10.1016/j.phrs.2012.10.020.

27. Segers, M.E., and Lebeer, S. (2014). Towards a better understanding of Lactobacillus rhamnosus GG-host interactions. *Microb. Cell Fact.* 13, S7. https://doi.org/10.1186/1475-2859-13-S1-S7.

28. Fitzgerald, J.W. (1976). Sulfate ester formation and hydrolysis: a potentially important yet often ignored aspect of the sulfur cycle of aerobic soils. *Bacteriol. Rev.* 40, 698-721. https://doi.org/10.1128/br.40.3.698-721.1976.

29. Harada, T., and Spencer, B. (1960). Choline suphate in fungi. *Microbiology* 22, 520-527. https://doi.org/10.1099/00221287-22-2-520.

30. van der Veen, J.N., Kennelly, J.P., Wan, S., Vance, J.E., Vance, D.E., and Jacobs, R.L. (2017). The critical role of phosphatidylcholine and phosphatidylethanolamine metabolism in health and disease. *Biochim. Biophys. Acta. Biomembr.* 1859, 1558-1572. https://doi.org/10.1016/j.bbamem.2017.04.006.

31. Zeisel, S.H., and da Costa, K.-A. (2009). Choline: an essential nutrient for public health. *Nutr. Rev.* 67, 615-623. https://doi.org/10.1111/j.1753-4887.2009.00246.x.

32. Shetty, S.A., Kostopoulos, I., Geerlings, S.Y., Smidt, H., de Vos, W.M., and Belzer, C. (2022). Dynamic metabolic interactions and trophic roles of human gut microbes identified using a minimal microbiome exhibiting ecological properties. *ISME J.* 16, 2144–2159. https://doi.org/10.1038/s41396-022-01255-2.

33. Zhao, T., Yue, H., Peng, J., Nie, Y., Wu, L., Li, T., Niu, W., Li, C., Zhang, Z., Li, M., and Ding, K. (2023). Degradation of xylan by human gut Bacteroides xylanisolvens XB1A. *Carbohydr. Polym.* 315, 121005. https://doi.org/10.1016/j.carbpol.2023.121005.

34. Despres, J., Forano, E., Lepercq, P., Comtet-Marre, S., Jubelin, G., Yeoman, C.J., Miller, M.E.B., Fields, C.J., Terrapon, N., Le Bourvellec, C., et al. (2016). Unraveling the pectinolytic function of Bacteroides xylanisolvens using a RNA-seq approach and mutagenesis. *BMC Genom.* 17, 147. https://doi.org/10.1186/s12864-016-2472-1.

35. Reunanen, J., Kainulainen, V., Huuskonen, L., Ottman, N., Belzer, C., Huhtinen, H., de Vos, W.M., and Satokari, R. (2015). Akkermansia muciniphila Adheres to Enterocytes and Strengthens the Integrity of the Epithelial Cell Layer. *Appl. Environ. Microbiol.* 81, 3655–3662. https://doi.org/10.1128/AEM.04050-14.

36. Liu, M.-J., Yang, J.-Y., Yan, Z.-H., Hu, S., Li, J.-Q., Xu, Z.-X., and Jian, Y.-P. (2022). Recent findings in Akkermansia muciniphila-regulated metabolism and its role in intestinal diseases. *Clin. Nutr.* 41, 2333–2344. https://doi.org/10.1016/j.clnu.2022.08.029.

37. Flint, H.J. (2012). The impact of nutrition on the human microbiome. *Nutr. Rev.* 70, S10-S13. https://doi.org/10.1111/j.1753-4887.2012.00499.x.

38. Jiménez, N.E., Acuña, V., Cortés, M.P., Eveillard, D., and Maass, A.E. (2023). Unveiling abundance-dependent metabolic phenotypes of microbial communities. *mSystems* 8. e00492-23. https://doi.org/10.1128/msystems.00492-23.

39. Herrgård, M.J., Fong, S.S., and Palsson, B.Ø. (2006). Identification of Genome-Scale Metabolic Network Models Using Experimentally Measured Flux Profiles. *PLoS Comput. Biol.* 2, e72. https://doi.org/10.1371/journal.pcbi.0020072.

40. Marinos, G., Kaleta, C., and Waschina, S. (2020). Defining the nutritional input for genome-scale metabolic models: A roadmap. *PLoS One* 15, e0236890. https://doi.org/10.1371/journal.pone.0236890.

41. Arias, N., Arboleya, S., Allison, J., Kaliszewska, A., Higarza, S.G., Gueimonde, M., and Arias, J.L. (2020). The Relationship between Choline Bioavailability from Diet, Intestinal Microbiota Composition, and Its Modulation of Human Diseases. *Nutrients* 12, 2340. https://doi.org/10.3390/nu12082340.

42. Gatarek, P., and Kaluzna-Czaplinska, J. (2021). Trimethylamine N-oxide (TMAO) in human health. *EXCLI J.* 20, 301–319. https://doi.org/10.17179/excli2020-3239.

43. De Filippis, F., Pasolli, E., and Ercolini, D. (2020). The food-gut axis: lactic acid bacteria and their link to food, the gut microbiome and human health. *FEMS Microbiol. Rev.* 44, 454-489. https://doi.org/10.1093/femsre/fuaa015.

44. Hillman, E.T., Lu, H., Yao, T., and Nakatsu, C.H. (2017). Microbial Ecology along the Gastrointestinal Tract. *Microbes Environ.* 32, 300–313. https://doi.org/10.1264/jsme2. ME17017.

45. Zhou, X., Kandalai, S., Hossain, F., and Zheng, Q. (2022). Tumor microbiome metabolism: A game changer in cancer development and therapy. *Front. Oncol.* 12.

46. Rautmann, A.W., and de La Serre, C.B. (2021). Microbiota's Role in Diet-Driven Alterations in Food Intake: Satiety, Energy Balance, and Reward. *Nutrients* 13, 3067. https://doi.org/10.3390/nu13093067.

47. Kolodziejeczyk, A.A., Zheng, D., and Elinav, E. (2019). Diet-microbiota interactions and personalized nutrition. *Nat. Rev. Microbiol.* 17, 742-753. https://doi.org/10.1038/s41579-019-0256-8.

48. Rodríguez, J.M., Murphy, K., Stanton, C., Ross, R.P., Kober, O.I., Juge, N., Avershina, E., Rudi, K., Narbad, A., Jenmalm, M.C., et al. (2015). The composition of the gut microbiota throughout life, with an emphasis on early life. *Microb. Ecol. Health Dis.* 26, 26050. https://doi.org/10.3402/mehd.v26.26050.

49. Cregut, M., Durand, M.-J., and Thouand, G. (2014). The Diversity and Functions of Choline Sulphatases in Microorganisms. *Microb. Ecol.* 67, 350–357. https://doi.org/10.1007/s00248-013-0328-7.

50. Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., et al. (2011). *Scikit-learn: Machine Learning in Python. Mach. Learn. Res.* 12, 2825-2830. PYTHON.

51. Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., et al. (2020). SciPy 1.0: fundamental algorithms for scientific computing in Python. *Nat. Methods* 17, 261-272. https://doi.org/10.1038/s41592-019-0686-2.

52. Caswell, T.A., Droettboom, M., Lee, A., de Andrade, E.S., Hoffmann, T., Klymak, J., Hunter, J., Firing, E., Stansby, D., Varoquaux, N., et al. (2022). `matplotlib/matplotlib: REL: v3.5.2`. Zenodo. https://doi.org/10.5281/zenodo.6513224.

53. Waskom, M. (2021). seaborn: statistical data visualization. *J. Open Source Softw.* 6, 3021. https://doi.org/10.21105/joss.03021.

54. Ebrahim, A., Lerman, J.A., Palsson, B.O., and Hyduke, D.R. (2013). `COBRAPy: Constraints-Based Reconstruction and Analysis for Python. BMC Syst. Biol.* 7, 74. https://doi.org/10.1186/1752-0509-7-74.

55. Satopaa, V., Albrecht, J., Irwin, D., and Raghavan, B. (2011). Finding a "Knee" in a Haystack: Detecting Knee Points in System Behavior. In 2011 31st International Conference on Distributed Computing Systems Workshops (IEEE), pp. 166-171. https://doi.org/10.1109/ICDCSW.2011.20.

56. Heinken, A., Hertel, J., Acharya, G., Ravcheev, D.A., Nyga, M., Okpala, O.E., Hogan, M., Magnúsdóttir, S., Martinelli, F., Nap, B., et al. (2023). Genome-scale metabolic reconstruction of *7,302* human microorganisms for personalized medicine. *Nat. Biotechnol.* 41, 1320-1331. https://doi.org/10.1038/s41587-022-01628-0.

57. Sahoo, S., and Thiele, I. (2013). Predicting the impact of diet and enzymopathies on human small intestinal epithelial cells. *Hum. Mol. Genet.* 22, 2705-2722. https://doi.org/10.1093/hmg/ddt119.

58. Löhne, A., and Weißing, B. (2016). Equivalence between polyhedral projection, multiple objective linear programming and vector linear programming. *Math. Methods Oper. Res.* 84, 411–426. https://doi.org/10.1007/s00186-016-0554-0.

59. Ciripoi, D., Löhne, A., and Weißing, B. (2018). A vector linear programming approach for certain global optimization problems. *J. Glob. Optim.* 72, 347–372. https://doi.org/10.1007/s10898-018-0627-0.

60. Seemann, T. (2014). Prokka: rapid prokaryotic genome annotation. *Bioinformatics* 30, 2068–2069. https://doi.org/10.1093/bioinformatics/btu153.

61. Van Rossum, G., and Drake, F.L. (2009). Python 3 reference manual (`CreateSpace`, `Scotts Valley`).

62. Inc. PT. Collaborative data science [Internet]. Montreal, QC: Plotly Technologies Inc.; 2015. Available from: https://plot.l
## STAR★METHODS

### KEY RESOURCES TABLE

<table><thead><tr><th>REAGENT or RESOURCE</th><th>SOURCE</th><th>IDENTIFIER</th></tr></thead><tbody><tr><td colspan="3"><strong>Deposited data</strong></td></tr><tr><td>Small epithelial intestinal cell model</td><td></td><td>https://www.thielelab.eu/in-silico-models/thielelab.eu/in-silico-models</td></tr><tr><td>EMBL GEMs</td><td></td><td>https://github.com/cdanielmachado/embl_gems.github.com/cdanielmachado/embl/_gems</td></tr><tr><td>Virtual metabolic human (VMH)</td><td></td><td>https://www.vmh.life</td></tr><tr><td colspan="3"><strong>Software and algorithms</strong></td></tr><tr><td>Custom code</td><td>This paper</td><td>https://github.com/Anna-cell/MO_GEMs_Score (https://doi.org/10.5281/zenodo.11494989)</td></tr><tr><td>Mocbapy</td><td>Budinich et al.<sup>23</sup></td><td>https://gitlab.univ-nantes.fr/mbudinich/mocbapy</td></tr><tr><td>Benpy</td><td>Budinich et al.<sup>23</sup></td><td>https://gitlab.univ-nantes.fr/mbudinich/benpy</td></tr><tr><td>CPLEX 12.10</td><td>IBM</td><td>https://www.ibm.com/fr-fr/analytics/cplex-optimizer</td></tr><tr><td>Python 3.7.13</td><td>https://www.python.org</td><td>N/A</td></tr><tr><td>scikit learn 1.0.2</td><td>Pedregosa et al.<sup>50</sup></td><td>https://scikit-learn.org/stable/</td></tr><tr><td>Scipy 1.7.3</td><td>Virtanen et al.<sup>51</sup></td><td>https://scipy.org/</td></tr><tr><td> matplotlib 3.5.2</td><td>Caswell T.A et al.<sup>52</sup></td><td>https:// matplotlib.org/</td></tr><tr><td>seaborn 0.11.2</td><td>Waskom<sup>53</sup></td><td>https://seaborn.pydata.org/</td></tr><tr><td>plotly 5.10.0</td><td>Plotly Technologies Inc.</td><td>https://plotly.com/</td></tr><tr><td>CarveMe 1.5.1</td><td>Machado et al.<sup>16</sup></td><td>https://github.com/cdanielmachado/embl_gems</td></tr><tr><td>Cobrapy 0.22.1</td><td>Ebrahim et al.<sup>54</sup></td><td>https://opencobra.github.io/cobrapy/</td></tr><tr><td>Kneed 0.8.5</td><td>Satopää et al.<sup>55</sup></td><td>https://pypi.org/project/kneed/</td></tr></tbody></table>

## RESOURCE AVAILABILITY

### Lead contact

Further information and requests for resources should be directed to and will be fulfilled by the lead contact, Damien Eveillard (damien.eveillard@univ-nantes.fr) or the other corresponding author, Anna Lambert (anna.lambert@univ-nantes.fr).

### Materials availability

This study did not generate new materials.

### Data and code availability

* This paper analyzes existing, publicly available data. These accession numbers for the datasets are listed in the key resources table.

* All original code has been deposited at https://github.com/Anna-cell/MO_GEMs_Score (https://doi.org/10.5281/zenodo.11494989) and is publicly available as of the publication date.

* Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

## METHOD DETAILS

### Genome-scale metabolic models

The 331 bacteria models used in this study are CarveMe reconstructions extracted from the publicly available EMBL GEMs database.¹⁶ These bacteria were selected based on their description as gut bacteria in the Virtual Metabolic Human (VMH) database.¹⁴ The bacteria described in the VMH database's metadata as "Pathogen," "Opportunistic pathogen," or "Putative Pathogen" were joined in the "Pathogen" ecological type (n = 46). The ones described as "Probiotic" or "Probiotic potential" were joined in the "Probiotic" ecological type (n = 18). The "Commensal" (n = 123) and "uncharacterized" (n = 144) ecological types were conserved as is. Although the analysis described in this paper
have been applied on models reconstructed with CarveMe, the method is generic and is compatible with high-quality GEMs such as those available via AGORA.⁵⁶

The small intestine epithelial cell (sIEC) model relies on previously published results⁵⁷ and includes 1282 reactions and 844 metabolites. Its objective, referred to as the enterocyte maintenance, is an adaptation of the Recon 2 biomass reaction, in which lipids and vitamins are added and dNTPs are removed. This models the enterocyte metabolic and cellular maintenance, without replication. The namespace of the exchange reactions (i.e., the nomenclature of the model elements) was adapted to the EMBL models (i.e., BiGG's namespace) for model compatibility. The sIEC model includes two external compartments: the blood and the lumen. In constructing the ecosystem model, the reactions controlling the apparition and disappearance of metabolites in the blood were not considered exchange reactions, in contrast with those in the lumen.

## Ecosystem model: Pool compartment, diet

bacterial model, and the sIEC's model were joined to exchange metabolites through a **pool compartment** to build an ecosystem model. Mathematically, all models' stoichiometric matrices were diagonally assembled into a new one, in which a pool compartment was added. The original exchange reactions of the models became transport reactions (TR) from the organism's external compartments to the pool compartment. All exchange reactions in at least one of the models were duplicated to form the pool's exchange reactions (ER) group. TR were unconstrained to enable free transit of the metabolites between the pool and the organism's external compartment, and the media constraint was applied by restricting the bounds of ER.

This metabolic model of the ecosystem was built using a Python version of mocba²³ (i.e. mocbapy), and further constrained to fit the nutritional conditions described above on the pool compartment.

In this study, the **unconstrained diet** consists of the absence of constraint on the exchange reactions of the pool. The **Western diet** (WD) and **Protein diet** (PD) were extracted from Heinken et al.²⁴ WD is high in simple sugar (47%) and fat (35%) and is low in fiber (3%) and proteins (15%). PD is high in protein (30%) and balanced in fat (20%), simple sugar (25%), and fiber (25%). For WD and PD conditions, exchange reactions between the sIEC and the blood were constrained based on the Average American Diet (AAD).⁵⁷

When constraining the pool's content strictly to the diet information, with every other reaction blocked, the bacteria could not grow. Therefore, the exchange reactions not described in the diet were left unconstrained.

## Multi-objective linear programming: Pareto front

Multi-objective linear programming is a mathematical optimization technique that optimizes multiple conflicting objectives by finding optimal solutions representing the trade-offs between the different objectives. This trade-off is formally a Pareto front describing the optimal behaviors of the resulting metabolic ecosystem, in the sense that no increase in an objective can be done without affecting (i.e., decreasing) others. For identifying the Pareto front for an ecosystem composed of multiple metabolic models, we used the mocbapy python package. It translates the metabolic model of the ecosystem into a multi-objective linear problem and solves it using a Python adaptation of Bensolve⁵⁸,⁵⁹ to identify the Pareto front. This Pareto front is a set of extreme points in the objective space. This case study's objectives are the bacteria and the sIEC biomasses.

## Linear problem formalisms

A metabolic model comprises a stoichiometric matrix (S) representing metabolite-reaction relationships, fluxes (v) representing reaction rates, and bounds (and) defining flux constraints. A solution for a usual mono-objective FBA is obtained by solving the following linear problem:

$$ \begin{array}{ll} \text{maximize} & z = c^T v \\ v & \in \mathbb{R}^n \end{array} $$

subject to

$$ \begin{array}{l} S_v = 0 \\ l_i \le v_i \le u_i \quad i = 1, \dots, n \end{array} $$

where $c^T v$ is a linear combination of fluxes representing the objective function.

In a multi-objective linear problem, the stoichiometric matrix is organized in several compartments, being able to exchange metabolites through a pool compartment. There are as many optimized objective functions as organisms modeled in the ecosystem. As described inBDinich et al.,²³ the multi-objective linear problem solved in this instance can be defined as

$$ \begin{array}{l} \text{maximize} \\ v \in \mathbb{R}^n \end{array} \quad \begin{pmatrix} f_1 \\ \vdots \\ f_k \end{pmatrix} = \begin{pmatrix} c_1^T v \\ \vdots \\ c_k^T v \end{pmatrix} $$
subject to

$$ 
\begin{array}{l}
S^a v = 0 \\
I_i \le v_i \le u_i \quad i = 1, \dots, \bar{n}
\end{array}
$$

where $(f_1, ..., f_k)^T$ are the objective functions of the k organisms, $S^a$ is the ecosystem stoichiometric matrix, and $\bar{n}$ is the total number of reactions (i.e., the sum of reactions of each organism and exchange reactions from the pool compartment).

## Interaction score

The Pareto front of **bi-objective** problems is described in a two-dimensional space, each axis describing the possible values for one objective. To calculate an interaction score, the values of each dimension are normalized by the maximal value of their respective objective when the organism is alone. The maximal growth of the bacterium and the enterocyte alone are added to the Pareto front as (0, 1) and (1, 0), respectively.

The interaction score is inferred from the area under the normalized Pareto front (AUC<sub>P</sub>) curve, from which we subtract the area under the curve of the non-interaction front (AUC<sub>NI</sub>).

$$ S = \text{AUC}_{\text{P}} - \text{AUC}_{\text{NI}} $$

In a problem of **more than two objectives**, the same normalization is applied to the extreme points defining the Pareto front. The points of maximal growth when alone of each organism are added accordingly, as well as the origin point (a vector with zero values for all dimensions). The convex envelope of this set of points and its hypervolume are calculated using Scipy's ConvexHull function.⁵¹ The interaction score is the hypervolume, thus calculated from which is subtracted the non-interaction hypervolume.

When alone, all models were solved with CPLEX. Multi-objective problems were solved with benpy, a Python adaptation of bensolve.

## Interaction types

In **two dimensions**, the various forms taken by the Pareto front were discriminated into categories (illustrated in Figure S2) based on four conditions.

* Sign: Sign of the score.

* E+: A solution exists on the Pareto front where the biomass value of the enterocyte in the ecosystem is higher than that of the enterocyte alone.

* B+: A solution exists on the Pareto front where the biomass value of the bacterium in the ecosystem is higher than the biomass value alone.

* E+B+: A solution on the Pareto exists where the biomass values of both organisms are at their highest.

They were not studied in problems of higher dimensions.

<table><caption>Table : Interaction types inferred from Pareto appearance</caption><thead><tr><th>Sign</th><th>E+</th><th>B+</th><th>E+B+</th><th>Interaction type</th></tr></thead><tbody><tr><td>+</td><td>False</td><td>False</td><td>False</td><td>Neutralism</td></tr><tr><td>-</td><td>False</td><td>False</td><td>False</td><td>Competition</td></tr><tr><td>+</td><td>True</td><td>False</td><td>False</td><td>Favors enterocyte</td></tr><tr><td>+</td><td>False</td><td>True</td><td>False</td><td>Favors bacteria</td></tr><tr><td>+</td><td>True</td><td>True</td><td>False</td><td>Limited mutualism</td></tr><tr><td>+</td><td>True</td><td>True</td><td>True</td><td>Mutualism</td></tr></tbody></table>

## Sampling of the Pareto front

A homogeneous sampling of 100,000 solutions was performed along the Pareto front, describing the interaction between LGG and the enterocyte. Each sample was obtained by conducting a Flux Balance Analysis (FBA) with fixed biomasses of both organisms at the chosen Pareto point. The ecosystem biomass (sum of each biomass) was added to the resulting sampling, and a Spearman correlation matrix was generated. The exchange reactions with a correlation higher than 0.95 (secreted) or lower than -0.95 (uptaken) with the ecosystem biomass were considered drivers of the mutualistic interaction between LGG and the enterocyte (Figure 3B).

## Exchanged metabolites and obligatory reaction

Exchanged metabolites are metabolites for which the secretion by an organism is correlated (Spearman, corr > 0.5) to their absorption by the other organism. Obligatory reactions cannot be inactive based on a Flux Variability Analysis (FVA), which determines the upper and lower possible value for flux in a given condition.
<figure><img src="image_17.png" /></figure>

## Building a five-dimension ecosystem

Based on Shetty et al.,³² 16 bacteria strains were selected to represent a minimal human gut microbiome. First, the sequenced genome of the strains were annotated with prokka.⁶⁰ Then, CarveMe¹⁶ was used to generate metabolic models for each organism (See genome references in Table S3). For computational reasons, a multi-objective problem could be solved with up to 10 bacteria models in addition to the enterocyte (Figure S5). We observed that some bacteria were in total competition, meaning that if one bacteria was growing, the others had null growth. This led to the selection of four cohabiting bacteria: Akkermansia muciniphila ATCC BAA-835, Bacteroides xylanisolvens HMP 2_1_22, Faecalibacterium prausnitzii A2-165 and Ruminococcus bromii ATCC 27255. They were joined with the enterocyte in an ecosystem constrained with a WD, and the extreme solution points in the inferred Pareto were retrieved.

## Five-dimensional principal component analysis (PCA)

In the five-dimensional ecosystem, 34 extreme points were identified on the Pareto front using mocbapy. To integrate solutions embedded in the solution space, the model was converted from mocbapy format to cobapy. Then, the sampling function from cobapy was applied to retrieve 3000 random solutions. In each solution, only the objective values for each organism were kept for the PCA.

The data was processed and the model was built using scitlearn⁵⁰ tools (preprocessing.StandardScaler() and decomposition.PCA()). The explained variance of the five principal components were 30%, 21%, 18%, 17% and 14%, respectively. As the principal components three and four had non-negligible variance explanations, the knee point was found using the kneed python package.⁵⁵ The knee point is the point of maximum curvature in a function. Only the first two principal components were explored as the knee-point was equal to two.

## QUANTIFICATION AND STATISTICAL ANALYSIS

All statistical tests were realized using Python 3.7.13.⁶¹ Significance was determined from a p-value inferior to 0.05. The comparison of interaction score values for ecological types (Figure 2C) was assessed by pairs of ecological types with a Mann-Whitney U test using numpy.stats.mannwhitneyu(). The ecological types comprised n = 46 pathogens, n = 18 probiotics, n = 123 commensals and n = 144 uncharacterized. Every correlation was calculated with pandas.corr(method = "spearman"). A correlation was considered of interest if its absolute value was equal or over 0.95. To calculate the score, the AUC (in two-dimensional Pareto fronts) was calculated with scikit learn 1.0.2⁵⁰ (sklearn.metrics.auc) and the hypervolume (Pareto fronts of three or more dimensions) was calculated with numpy 1.7.3⁵¹ (scipy.spatial.ConvexHull.volume). In case on non-monotonic two-dimensional Pareto fronts, the AUC excluding the last point was first computed, followed by the AUC for the last two points. The second AUC was subtracted from the first to determine the Pareto front AUC. This is possible because in this context, a Pareto front can be non-monotonic only by its last point. FBA and FVA were performed using cobapy 0.22.1.⁵²

The figures were built using matplotlib 3.5.2,⁵² seaborn 0.11.2⁵³ and plotly 5.10.0.⁶²