How gut bacteria relate to metabolism changes in Parkinson's disease

To ensure that all modeling predictions were relevant to PD, we assessed only metabolites in the microbiome-personalized WBMs with previously known, replicated associations with PD. Therefore, we used a meta-analysis of 74 metabolomic studies on PD patients,³² which identified 190 metabolites to be associated with PD diagnosis in more than one study. Of these 190 metabolites, 116 were present in the blood compartment of the WBMs and belonged to 9 biochemical superclasses, 24 classes, and 34 subclasses (Table S4).³² Organic acids (46 metabolites) were the most abundant superclass, followed by lipids (26 metabolites). Within the organic acids, carboxylic acids, which include amino acids and peptides, were the largest biochemical class with 39 metabolites. The selected lipid metabolites contained sixteen fatty acyls, which include fatty acids, and six steroid metabolites, including five bile acids. The analyzed metabolites were further involved in 42 metabolic subsystems, including alanine and aspartate metabolism (11 metabolites) and fatty acid oxidation (10 metabolites) were the most common metabolic subsystems. The associated metabolic subsystems captured 26% of all 161 metabolic subsystems in the VMH (virtual metabolic human) database,³¹ indicating that the curated metabolites capture a diverse range of metabolic pathways. Not all selected metabolites were present in the gut microbiota lumen. Only 72 metabolites could be produced in microbial pathways (Table S4). The remaining 44 metabolites, however, could be indirectly influenced by their metabolic precursors and could thus still be influenced by gut microbiota. Taken together, we systematically selected 116 metabolites for further analysis in the WBM blood compartments, all of which had known replicated metabolomic associations with PD diagnosis. Although this set of metabolites likely represents a fraction of all PD-associated metabolites, it captures major parts of PD-associated metabolism, containing 61% of the known metabolites with replicated metabolomic PD associations, nine biochemical superclasses, and is involved in 26% of all metabolic subsystems in the VMH resource.

To identify which of these gut microbial influencers of the flux predictions could potentially explain the lower predicted blood fluxes in PD models, we extracted the potential microbial species flux contributions from the FBA predictions (Methods, Table S6), which were defined as the maximal change in predicted flux observed when a microbiome-personalized WBM was restored to its original state after the removal of a microbial species. All microbial species were then ranked based on their average metabolic contribution potentials to the predicted metabolic blood fluxes of the selected metabolites and pruned to find the smallest set of microbial species that could contribute to, on average, 95% of the total potential gut microbiome flux contributions (Methods, Table S6). The resulting pruned sets of flux-associated microbial species were much smaller than the original sets (Table 4), showing that most of the predicted gut microbiome influences were the result of a minority of the gut microbial species in the models.

Besides B. uniformis and B. vulgatus, our results also identified large metabolic contribution potentials for all six metabolites for F. prausnitzii, Eubacterium rectale, and Blautia wexlerae. Their potential metabolic contributions could again be attributed to their high relative abundances, which were the third highest (Table S1) of the mapped microbial species. In contrast to B. uniformis and B. vulgatus, these gut microbes had lower relative abundances in PD microbiomes compared to gut microbiomes of neurologically healthy individuals, meaning that their lower relative abundance in PD patients could have influenced the flux predictions. The correlation analysis further found moderate correlations with predicted fluxes for butyrate, myristic acid, and pantothenate (Figure 3(g)). The strongest correlations were found for F. prausnitzii with both butyrate and myristic acid blood fluxes. In contrast, weak correlations were found for L-leucine, leucylleucine, and nicotinic acid across the three microbial species (Figure 3(g)). These weaker correlations could be explained by the larger number of microbial species that contribute to L-leucine, leucylleucine, and nicotinic acid, diluting the relative metabolic influences on the predicted fluxes. Taken together, we found that lower abundances of F. prausnitzii, and to a lesser extent, E. rectale and B. wexlerae, potentially contribute to lower microbial production of blood metabolites in PD patients.

We also identified microbial contributors to the fluxes with higher relative abundances in PD patients. One of these contributors was Bifidobacterium longum, which had large metabolic contribution potentials for L-leucine, leucylleucine, pantothenate, and nicotinic acid (Figure 3(a,b,d,e)). When correlating the B. longum relative abundances with flux predictions for these metabolites, negative correlations were found (Figure 3(g)). This negative correlation indicates that a higher relative abundance of B. longum associated with lower gut microbiome contributions to L-leucine, leucylleucine, pantothenate, and nicotinic acid, and can be interpreted as B. longum potentially limiting other microbial producers of these metabolites. Another major microbial contributor to the predicted fluxes of the selected metabolites with higher relative abundances in PD microbiomes was Parabacteroides distasonis. In contrast to B. longum, P. distasonis had weakly positive correlations with its associated predicted metabolic fluxes (Figure 3(g)). This result could indicate that P. distasonis was less important than B. longum for the predicted reduction in the blood fluxes from controls to PD patients.

Other potential microbial influencers of the predicted fluxes were Escherichia coli and M. smithii. Despite these microbial species not being among the top microbial contributors to the predicted fluxes, they were top negative microbial correlates with the selected metabolites, except for butyrate (Figure 3(g)). In particular, high negative correlations were identified for M. smithii relative abundances and predicted fluxes for L-leucine and leucylleucine. When cross-matching this finding with previously found associations with PD diagnosis, we found that the relative abundances of both E. coli and M. smithii positively associate with PD status in the analyzed samples. This finding suggests that elevated relative abundances of E. coli and M. smithii in PD microbiomes could contribute to lower predicted levels of microbial metabolites in blood.

For butyrate production in blood, an optimal cluster of four microbial species was found consisting of B. uniformis, F. prausnitzii, B. wexlerae, and E. rectale. The combined relative abundances of these microbial species resulted in a Spearman correlation coefficient of 0.81 (Figure 4(a)), suggesting that these microbial species were the main microbial influencers of the predicted butyrate blood fluxes. Within this microbial cluster, F. prausnitzii, B. wexlerae, and E. rectale had lower relative abundances in the analyzed PD microbiomes, suggesting that these microbial species were likely the most important factors of the lower predicted butyrate blood fluxes in PD samples. This result further confirms B. uniformis as a major contributor to the predicted butyrate blood fluxes. Similarly to butyrate, a strong Spearman correlation was found for myristic acid blood fluxes and a microbial cluster consisting of B. uniformis, B. vulgatus, F. prausnitzii, and Prevotella copri (ρ = 0.77). As both F. prausnitzii and P. copri had lower relative abundances in PD microbiomes, these microbial species were likely major influencers of the lower predicted myristic acid fluxes. Next, for pantothenate, a cluster consisting of Bacteroides dorei, B. vulgatus, and F. prausnitzii was found. As B. dorei and B. vulgatus were not associated with PD microbiomes, the lower pantothenate fluxes were potentially influenced by F. prausnitzii. Finally, for nicotinic acid, the best correlating microbial combination contained Alistipes putredinis, Collinsella aerofaciens, Phocaeicola plebeius, and R. lactatiformans. Together, these microbial species negatively correlated with the predicted nicotinic acid fluxes (Figure 4(a)), indicating that a combined growth of these microbial species could contribute to lower nicotinic acid production in blood. Although the relative abundances of A. putredinis, C. aerofaciens, and P. plebeius were stable between PD microbiomes and controls, R. lactatiformans relative abundances were higher in PD microbiomes. This finding suggests that R. lactatiformans potentially contributed to the lower predicted nicotinic acid blood fluxes in PD patients, which aligns with the previously found association between R. lactatiformans relative abundances on predicted nicotinic acid blood fluxes (Figure 3(g)). Taken together, we identified combinations of gut microbial contributors to the predicted fluxes that best explained the flux predictions and found potential influencers of the predicted lower blood fluxes in PD patients. Overall, F. prausnitzii was most often found to explain the PD-associated predicted blood fluxes (Figure 4(b)), suggesting that this microbial species is a potential key microbe in gut microbiome influences on blood metabolism in PD.

Participant recruitment has been described before by Wallen et al.²⁹ Briefly, participants were enrolled in the city of Birmingham, Alabama and the surrounding region between October 2018 and March 2020. The eligibility criteria for PD cases have been a confirmed PD diagnosis by a specialist at the Movement Disorder Clinic at the University of Birmingham, Alabama, and informed consent of the PD patient. Healthy controls were recruited from spouses or friends of the PD patients. Eligibility criteria for these controls were the absence of a PD diagnosis, and a lack of the following disorders: Alzheimer's disease, dementia, REM sleep behavior disorder, multiple sclerosis, ataxia, dystonia, autism, epilepsy, stroke, amyotrophic lateral sclerosis, self-reported schizophrenia, and bipolar disorder. In addition, healthy controls were preferably aged 50 or older to ensure an age-matched cohort. Enrollment and sample collection was stopped in March 2020 to avoid confounding influences of COVID-19 and stress related to the COVID-19 pandemic. All sample processing steps, including DNA extraction, shotgun sequencing, quality and control steps, taxonomic profiling, and the collection of metadata, has been carried out blinded by Wallen et al.

The gut metagenomics reads and the associated metadata were retrieved from a publicly available archive: https://zenodo.org/records/7246185↗. The retrieved gut metagenomics were already preprocessed and taxonomically profiled and have been comprehensively described here.²⁹ Briefly, the metagenomics data were obtained from stool samples acquired from 490 PD patients and 234 neurologically healthy controls. The stool samples were sequenced by shotgun sequencing on an Illumina NovaSeq 6000, which included positive and negative controls. Human contamination has been removed by mapping the quality-controlled reads to the GRCh38.p13 human reference genome (https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000001405.39/↗). Taxonomic profiling has been performed with MetaPhlAn3,⁹¹ using the CHOCOPhlAn V30 metagenomic reference database (http://cmprod1.cibio.unitn.it/biobakery3/metaphlan_databases/↗). The participant metadata was previously obtained from questionnaires without interference from investigators (https://static-content.springer.com/esm/art%3A10.1038%2Fs41467-022-34667-x/MediaObjects/41467_2022_34667_MOESM1_ESM.pdf↗.)

Microbial species in the retrieved gut metagenomics were mapped onto genome-scale metabolic reconstructions in the combined AGORA2²⁴ and APOLLO²⁵ resources, using the MARS pipeline.⁹² After mapping, the metagenomic reads were normalized by calculating the relative read abundances per sample. To limit false positives, sample relative abundances below the threshold of 0.00001% were further removed, whereafter the read abundances were normalized again. The mapping coverages, i.e., the fraction of available metagenomic data in a sample after mapping, were calculated from the original and mapped relative abundances after removing the reads below the relative abundance threshold.

Genome-scale metabolic reconstructions for the mapped microbial species were retrieved from AGORA2 (version 2.01, https://www.vmh.life/files/reconstructions/AGORA2/↗.) and APOLLO²⁵ (version 1.0 https://doi.org/10.7910/DVN/PIZCBI↗). As the metabolic reconstructions in AGORA2 and APOLLO are strain-resolved, pan-species models were generated to match the mapped species in the metagenomic data. The pan-species models were generated using the createPanModels.m CobraToolbox¹⁸ function and contained all unique reactions, metabolites, and genes in the set of AGORA2 and APOLLO strains within each microbial species in the mapped metagenomics data.

Personalized metabolic models of gut microbiome communities were created using the mgPipe module in the Microbiome Modelling toolbox 2.0.⁹³ The creation of microbiome community models has been described elsewhere.⁹⁴ Briefly, pan-species models that mapped onto the gut metagenomics data in a sample were connected to a shared microbiota lumen, while transport reactions were added to enable metabolic exchanges between the pan-species models and the shared microbiota lumen. The mapped relative abundances were then integrated into the microbiome communities via the addition of a community biomass reaction (VMH reaction ID: communityBiomass), which lists each microbe's biomass reaction with the microbe's relative abundance as the stoichiometric coefficient. The community biomass reaction ensured proportional microbial growth for each integrated pan-species model per the mapped relative species abundances. Consequently, the community biomass reaction ensured that the metabolic contribution capacity of each microbial species to the host was proportional to their relative abundances.

We used generic sex-specific male and female WBMs²⁶ (version 1.04c) for the host part of the microbiome-personalized WBMs in this study. These organ-resolved WBMs captured the metabolism of 26 organs and 13 biofluid compartments, which were connected in an anatomically accurate manner. The male generic WBM included 83,395 reactions, 58,095 metabolites, and 105,936 constraints, while the female generic WBM had 85,892 reactions 60,537 metabolites, and 109,757 constraints. All parameters were set to their defaults.

The microbiome-personalized WBMs were created using the combineHarveyMicrobiota.m function, which joined the personalized microbiome community models with the large intestinal lumen compartment of the generic WBMs, as described previously.²⁶ This function also added transport reactions to enable metabolic exchanges between the microbiota and the host. Moreover, this function updated the coupling constraints for each microbe's reaction to their biomass reactions, to a value of 400, following previous WBM modeling studies.^26–28,95

Further parameterization of the WBMs was done by setting the upper and lower flux bounds of the Excretion_EX_microbiota_LI_biomass[fe] reaction to one mmol/day/person. This reaction excretes the microbiota_LI_biomass compound, which is produced by the communityBiomass reaction. Fixing the flux bounds of the Excretion_EX_microbiota_LI_biomass[fe] reaction ensured that all generated WBMs produced the same metabolic flux over their microbiome community, despite differences in community compositions. In addition, the upper and lower flux bounds of the Whole_body_objective_rxn were set to one mmol/day/person, following previous WBM studies.^26–28,96 The Whole_body_objective_rxn includes a linear combination of the whole-body organ biomass productions and associated coefficients for the proportional organ weights. This additional constraint ensured proportional organ biomass production between samples, corresponding to the organ weights.

The microbiome-personalized WBMs were further parameterized with dietary uptake reactions and reaction constraints from the Virtual Metabolic Human database (https://www.vmh.life↗.).³¹ Reactions and constraints were chosen for the average European diet, as sufficient information on dietary uptake was not available. The dietary uptake formulation represented the metabolite-level average intake of a one-day meal plan for a 70 kg adult³¹ and was added to the models using the setDietConstraints.m CobraToolbox function. The added diet contained 192 metabolites (Table S3), including trace amounts, i.e., dietary uptakes below 0.1 mmol/day/person, of 92 metabolites.

The metabolic fluxes in blood were predicted by performing FBA, which translates a genome-scale metabolic reconstruction into a numerical model and interrogates the model using a linear programming solver according to an objective function.³³ The objective function describes a linear optimization problem for finding extreme points in a constrained solution space of metabolic states as follows:maxvcT⋅v(1)s.t,S⋅v=0,lb≤v≤ub

The linear programming solver then finds a reaction flux vector v containing a maximized metabolic flux through an objective reaction, which is defined by a single non-zero coefficient in c. The maximized flux value must lie in the feasible solution space defined by constraints, such as the steady state constraint. Steady state is defined by the system of linear equations, S\cdotv=0, and enforces unchanging metabolite concentrations over time by ensuring that the amount of produced metabolites (products) equals the amount of consumed metabolites in each reaction (substrates). The metabolic products and substrates in each reaction are encoded by non-zero coefficients in the stoichiometric matrix, S. Substrates are represented by a negative coefficient, while reaction products are indicated by positive coefficients. In addition to the steady state constraint, each reaction in v was constrained by reaction-specific upper (ub) and lower (lb) bounds, derived from literature sources such as thermodynamic data and experimental measurements. If no reaction information was available, reactions were left unconstrained by setting lb and ub to arbitrary values of −1,000,000 and 1,000,000 mmol/person/day (see also Thiele et al.²⁶ for details). For the objective reactions, we introduced demand reactions to the WBM blood compartments (VMH name: [bc]), for each of the analyzed blood metabolites. Demand reactions are unbalanced left-sided reactions (e.g., DM_leu_L[bc]: 1 leu_L[bc] →∅ for L-leucine in blood), commonly used for predicting the maximal potential metabolite productions in a system. Demand reactions have been used in WBMs before to predict blood biomarkers of single gene defects⁹⁷ and gut microbiome flux influences associated with risk markers of Alzheimer's disease.²⁸ Demand reactions were separately introduced for each FBA using the addDemandReaction.m CobraToolbox function.¹⁸ FBA was performed for each model by running fba = optimizeWBModel(model) with the model.osenseStr = 'max' parameter to ensure maximization of the flux through the objective reaction. All default parameters were used for the linear optimizations. After performing FBA, the maximized fluxes were extracted from the fba.f solution field and rounded to six decimals to remove potential numerical artifacts that can appear in small numbers.

Outlier samples in the predicted fluxes were identified from a robust PCA on the z-transformed fluxes, performed with the alternating least squares algorithm. Potential outliers were identified by taking the Euclidean norm of the principal component scores for each sample, weighted by the explained variance of each principal component. This operation finds the average distance of a sample in a principal component to the center of the flux data. Values close to zero indicate samples for which average fluxes were predicted, while high values indicate potential outlier samples across the predicted fluxes. Together with a visual inspection of the first three PCA components, the ranked outlier scores were used to manually select outlier samples to be removed.

The Fisher exact test was used to quantify sample differences in PD patients and controls for the binary variables in the metadata. The fishertest.m Matlab function was used to calculate the Fisher-exact p-values, the associated odds ratios, and the 95% confidence intervals. For the numerical metadata, including age in years, BMI in kg/m², and total read counts, differences between PD patients and controls were tested by performing the Wilcoxon rank sum test, using the ranksum.m Matlab function.

Logistic regressions were performed to associate the predicted blood flux values (predictor) with disease status (response, PD/Control). The flux predictions were log2 transformed and normalized via Z-transformation for the regressions. The following confounders were included in the regressions: age in years at stool collection (Z-transformed), sex (male/female), alcohol intake (yes/no), laxative usage (yes/no), antihistamine usage (yes/no), mood medication usage (yes/no), pain medication usage (yes/no), sleep-aid medication (yes/no). Finally, the total sequence count (log10-transformed and Z-transformed) was added to account for technical biases in sequencing and read processing.

The confounder selection was informed by a previous confounder analysis on the same gut metagenomics by Wallen et al.²⁹ Their confounder analysis included a multivariable association study on both the available metadata and the processed gut metagenomics using the MaAsLin2 pipeline.⁹⁸ We also followed the decision of Wallen et al. to not include covariates for gastrointestinal issues, such as constipation, digestion issues, and rapid weight loss or weight gain, as these features are intrinsic to the PD disease phenotype and would mask PD status if included. As previously noted by Wallen et al., older age and male sex are also intrinsic risk factors for PD. Nevertheless, these household covariates were included to estimate general associations with PD status that are also representative of women and younger PD patients. The influence of age and sex on the flux predictions was evaluated by performing a second set of logistic regressions on the predicted fluxes associated with PD status. Age-related influences were investigated by adding an interaction term for flux and age to the logistic regressions as predictors of interest. Sex differences were subsequently analyzed using logistic regressions with an interaction term for flux and sex as the predictor of interest. The multiple logistic regressions were performed with the fitglm.m Matlab function, using all standard parameters. 95% confidence intervals were obtained using the Wald method by running the coefCI.m Matlab function. Finally, the regression p-values were adjusted to account for multiple testing using the Benjamini–Hochberg procedure,⁹⁹ by running the mafdr.m Matlab function with the 'BHFDR' option set to true.

The predicted PD-associated blood fluxes were systematically validated by qualitatively comparing the direction of the predicted associations with previous literature. Validation studies were selected by consulting only studies that were included in the same meta-analysis³² used for selecting the analyzed metabolites with replicated PD associations. All included studies for modeling validation passed quality assessments by having QUADOMICS¹⁰⁰ scores of eight or higher and a Newcastle-Ottawa Scale¹⁰¹ score of seven or higher.

Microbial influencers of the predicted fluxes were identified by retrieving shadow prices in the FBA solutions. Shadow prices are a routinely used metric for identifying the value of a metabolite toward an objective function. More precisely, shadow prices for a metabolite i, πi, represent the inverse ratio of an incremental relaxation of the steady state constraint for a metabolite i, δbi, and the resulting maximal change in the maximized flux through the objective reaction, δZ,(2)πi=−δZiδbi

To identify the microbial species that could have contributed to the flux predictions, we extracted the shadow prices of their microbial biomass metabolites. Microbial biomass metabolites have the useful characteristic that they are coupled through all associated reactions within the microbial pan model via flux capacity (coupling) constraints,⁹⁵ which enforce a minimal flux through the microbial reactions of 400 times the flux through the biomass reaction. A consequence of these biomass reaction couplings is that if a microbial species produces any metabolite in the pathway of the predicted metabolites in blood, a non-zero shadow price is found. Therefore, the shadow prices of the microbial biomass metabolites indicate which microbial species can contribute to the predicted fluxes. Practically, the shadow prices were automatically computed when performing FBA and were extracted from fba.y after performing fba = optimizeWbmodel(model). Specifically, the shadow prices of the microbial pan biomass metabolites were obtained, e.g., panBacteroides_uniformis_biomass[c]. To get the number of microbial influencers for each optimized reaction and each microbial species, we counted the number of microbial species with non-zero shadow prices for their biomass metabolite, e.g., panBacteroides_uniformis_biomass[c]. Here, a non-zero pan-biomass shadow price means that increased or decreased biomass production, i.e., increased flux through the biomass reaction of a microbial species, would result in an altered flux prediction. A numerical tolerance for the shadow prices of 1e-6 was used to identify shadow price values of zero. Finally, the total number of microbial influencers for the investigated reactions was calculated by taking the union of microbial species with non-zero shadow prices across the cohort. Taken together, the shadow price values of pan biomass metabolites were used to identify which microbial species could influence the predicted fluxes through the contribution of microbial metabolites, including precursors of the analyzed metabolites in blood.

The microbial species contribution potential was defined as the greatest possible change in the predicted flux between an analyzed microbiome-personalized WBM (mWBMfull) and a hypothetical reduced microbiome-personalized WBM, where the microbe of interest was removed (mWBMreduced). Crucially, we assumed that the microbial abundances in mWBMreduced are identical to mWBMfull, except for the removal of the microbial species of interest. This assumption enabled the prediction of how much each microbial species in the WBMs could contribute to the predicted blood fluxes in an ideal scenario. The simplest method for calculating the microbial species flux contribution potentials would be to iteratively predict the change in blood fluxes after removing each microbial species in all WBMs. However, as this approach would be prohibitively expensive to compute, we used an alternative computationally efficient approach that extracted the microbial species contribution potentials from the already calculated FBA solutions, thus eliminating the need for additional FBA predictions. As previously described, the shadow prices of microbial biomass metabolites in WBMs reflect the ratio of the maximal change in predicted flux, δZ, with a change in the availability of microbial biomass compound i, δbi, in the WBM. Each microbial biomass compound in the WBMs is produced through a microbial biomass reaction. Thus, a breakage of the steady state of biomass metabolite bi necessitates a corresponding (steady state breaking) change in flux through the microbial biomass reaction. Conversely, a breakage of steady state by a change in flux through a microbial biomass reaction will result in increased availability of the microbial biomass metabolite. Note that the availability of the microbial biomass only increased in this scenario because the artificial removal of microbial biomass metabolites in the communityBiomass reaction was constant, which was the result of enforcing the excretion of the output of the communityBiomass reaction, i.e., the total gut microbiome biomass production, to one mmol/day/person. Another consequence of fixing the production and excretion of the communityBiomass compound was that the steady state flux values of each microbial biomass reaction were fixed, as the defined constant removal of microbial biomass metabolites necessitates an exact opposite production rate of microbial biomass metabolites. The consequence of this linear relationship between changing biomass availability and corresponding changing fluxes through a microbial biomass reaction was that the maximal change in blood fluxes could be calculated from the associated shadow price value, πbiomass,i, and the change in biomass reaction flux. δvbiomass,i:(3)δZi=−πbiomass,i⋅δvbiomass,i

Note that this formula is a reformulation of the modified shadow price formula with the change in biomass reaction flux replacing the change in metabolite availability,, in the original formula. Another change was the multiplication of the shadow price value by −1, to reflect that a steady state-breaking increase in microbial metabolic contributions to the flux predictions resulted in increased blood flux productions and vice versa. Now, the maximal change in blood fluxes when removing and reintroducing a microbial species of interest could be simulated by setting the microbial biomass reaction flux ofto zero and the biomass reaction flux ofto the microbial biomass reaction flux value in the corresponding FBA solution. As the chosen microbial biomass reaction flux inwas zero, the potential microbial flux contributions could be calculated by multiplying the shadow price of a microbial biomass metabolite by the flux through the microbial biomass production reaction. Taken together, by multiplying the shadow prices of biomass metabolites with their corresponding biomass reaction fluxes, microbial flux contribution potentials were predicted in a computationally efficient manner. δ b i mWB M reduced mWB M full mWB M reduced

To identify potentially relevant microbial flux contributors, the bootstrapped mean average flux contribution was calculated for all microbial flux contributions to the six selected metabolites. Bootstrapping was performed with replacement for 50,000 samples for all potential microbial flux contributions and for each of the selected metabolites. Then, the gut microbial species were ranked on their mean potential flux contributions to the blood production of each metabolite and microbial subsets were obtained by filtering on the top ranked microbial species that together had a potential flux contribution of 95% of the sum of all flux-associated gut microbial species for each of the six selected metabolites. The pruned lists of potential microbial contributors to the predicted fluxes were then correlated to their respective flux predictions by performing Spearman correlations. Potential microbial clusters that best explained the predicted fluxes were found on the pruned lists of flux-associated microbial species using a "brute-force" method of i) calculating the summed relative abundances of all possible pairs, triplets, quartets, and quintets of flux-associated microbial species, and 2) performing Spearman correlations for all microbial combinations against the associated flux predictions. Optimal sizes of microbial combinations were then selected by finding the largest top microbial combination that improved the correlation strength with the associated fluxes by. | | ρ > 0.05

Modeling work, including the construction and analysis of the microbiome WBMs, was done in MATLAB 2020b and MATLAB 2024b (Mathworks^TM). In addition, the parallel computing toolbox, statistics and machine learning toolbox, and the bioinformatics toolbox were utilized to generate and analyze microbiome WBMs. The generation and analysis of the microbiome-personalized WBMs were performed using functions from the COBRA toolbox v3,¹⁸ including the microbiome modeling⁹³ and PSCM⁹³ toolboxes. The IBM ILOG CPLEX 12.10 linear solver (IBM Inc.) was used for performing FBA on the investigated metabolites. All visualizations were carried out in R 4.4.2, using the tidyverse package.¹⁰²