Personalized whole‐body models integrate metabolism, physiology, and the gut microbiome

We developed a novel, iterative approach to assemble WBM reconstructions by leveraging an existing generic human metabolic reconstruction, Recon3D (Brunk et al, 2018), reconstruction algorithms (Vlassis et al, 2014), omics data (Wishart et al, 2013; Kim et al, 2014; Uhlen et al, 2015), and manual curation of more than 600 scientific literature articles and books (Fig 1A). The female WBM, termed Harvetta 1.0, was generated by starting with a meta‐reconstruction, composed of a set of 30 identical Recon3D models connected through anatomically consistent biofluid compartments (Fig 1B, Materials and Methods, Reconstruction details). Similarly, the male WBM, termed Harvey 1.0, started from a meta‐reconstruction with one copy of Recon3D for each of its 28 organs, tissues, and cell types (Table 1). Note that we will refer in the following to all cell types, tissues, and organs collectively as organs. We extensively reviewed organ‐specific literature for genes, reactions, and pathways known to be present in the different organs (Materials and Methods, Reconstruction details, Table EV1). This review effort yielded 9,258 organ‐specific reactions, which we added to a core reaction set, defining reactions that had to be present in the WBM reconstructions. We removed 912 organ‐specific reactions reported to be absent in the literature (Materials and Methods, Reconstruction details, Table EV1). Moreover, evidence for the presence of organ‐specific proteins for all organs except the small intestine was obtained from the human proteome map (Kim et al, 2014) and the human proteome atlas (Uhlen et al, 2015) (Materials and Methods, Reconstruction details, Tables EV2 and EV3). For four organs, metabolic reconstructions have been published, namely the red blood cell (Bordbar et al, 2011b), adipocyte (Bordbar et al, 2011a), small intestine (Sahoo & Thiele, 2013), and liver (Gille et al, 2010). Thus, we added the respective organ‐specific reactions to the core reaction set, after mapping the reaction identifiers to the Virtual Metabolic Human (VMH, www.vmh.life) reaction identifiers (Materials and Methods, Reconstruction details). The literature was also reviewed for the absence of cellular organelles in the different organs (Materials and Methods, Reconstruction details, Table EV4), which were removed accordingly from the relevant organs in the meta‐reconstructions. As a next step, we added transport reactions to the core reaction set, based on organ‐specific transporter expression data collected from the literature (Materials and Methods, Reconstruction details, Table EV5; Sahoo et al, 2014), and organ‐specific sink reactions for metabolites known to be stored in different organs (Materials and Methods, Reconstruction details, Tables EV6 and EV7). We used metabolomics data from 16 different resources to enforce the presence of metabolites detected in the different biofluid compartments (Fig 1, Materials and Methods, Reconstruction details, Table EV8), and added the corresponding reactions to the core reaction set. Furthermore, we used literature information to define metabolites that are known to cross, or not, the blood–brain barrier and either added them to the core reaction set or removed them from the meta‐reconstructions, respectively (Materials and Methods, Reconstruction details, Table EV9). Additionally, we included the dietary uptake reactions to the core reaction set for metabolites that have been identified in food (Materials and Methods, Reconstruction details, Table EV10). Finally, to enable the integration of the gut microbiome with the WBM reconstructions (see below), we added sink reactions to the core reaction set for metabolites known to be produced by human gut microbes (Table EV11). Each organ contains a biomass maintenance reaction representing the macromolecular precursors (e.g., amino acids) required for organ maintenance (Materials and Methods, Reconstruction details). To represent the energy required to maintain the body's cellular function and integrity, we added a whole‐body maintenance reaction to both meta‐reconstructions, in which each organ biomass maintenance reaction is weighted based on its respective organ weight in the reference man or woman (Snyder et al, 1975) (Materials and Methods, Reconstruction details).

The collected information, in the form of the core reaction set, and the tailored meta‐reconstructions were used as inputs for a model extraction algorithm (Vlassis et al, 2014) to generate sex‐specific draft WBM reconstructions (Fig 1A, Materials and Methods, Reconstruction details). This network extraction algorithm returns a compact subnetwork containing all core reactions and a minimal number of additional reactions necessary to render the subnetwork flux consistent, i.e., each network reaction can carry a non‐zero flux value. As the algorithm adds a minimal number of reactions to the subnetwork, it should be evaluated whether such additions can be supported or rejected based on biological evidence. A similar approach has been suggested for gap filling (Rolfsson et al, 2011). Consequently, we revisited the literature for evidence to support reaction inclusion and either expanded the organ‐specific core reaction set or removed organ‐specific reactions from the meta‐reconstructions as appropriate (Fig 1A). We then generated new draft WBM reconstructions using the model extraction algorithm with this updated input. This cycle was iterated over a hundred times, focusing on different organs and metabolic pathways each time, which was necessary due to the large size of the reconstructions (Materials and Methods, Reconstruction details). Throughout the reconstruction process, we followed an established quality control and quality assurance protocol (Thiele & Palsson, 2010) and performed core tests proposed by the systems biology community (preprint: Lieven et al, 2018).

The final sex‐specific WBM reconstructions account for 83,521 and 81,094 reactions for Harvetta and Harvey, respectively (Fig 1C, Tables EV12 and EV13). Harvetta contained more reactions than Harvey due to two more sex‐specific organs (Fig 1C and D). In total, 73,215 reactions were shared between both sexes, while 7,879 were unique to Harvey and 10,306 unique to Harvetta. Excluding sex organs, < 4% of all WBM reactions were sex‐specific, mostly involving alternative transport reactions (Dataset EV2: 3.6). Overall, these WBM reconstructions account for about 84% of the body weight. The remaining 16% correspond mostly to bones and connective tissue (Snyder et al, 1975; Fig 1C), both of which are as yet insufficiently represented with corresponding reactions in Recon3D. The resulting, sex‐specific WBM reconstructions comprehensively capture human whole‐body metabolism consistent with current knowledge.

A compendium of organ‐specific reconstructions was extracted from the WBM reconstructions (Figs 1B and 2, Tables 1 and EV14). All organ‐specific reconstructions can be obtained from the VMH at www.vmh.life. The technical quality of each organ was assessed using a comprehensive test suite, implemented in the COBRA Toolbox (Heirendt et al, 2019), that checks for adherence to quality standards developed for metabolic reconstructions (Thiele & Palsson, 2010; Brunk et al, 2018; preprint: Lieven et al, 2018). These tests included, e.g., testing for realistic ATP yields on various carbon sources under aerobic and anaerobic conditions or the ability to perform defined metabolic tasks (Dataset EV2: 3.7). A typical organ contained 2,947 ± 1,980 reactions, 2,076 ± 1,081 metabolites, and 1,298 ± 245 genes (average ± SD, Fig 2A). As one would expect, the organs with the most reactions were liver, kidney, and colon (Fig 1D). On average, 72% of reactions in each organ were gene‐associated, when excluding transport and exchange reactions (Fig 2A). About 10% of all metabolites were organ‐specific metabolites (Fig 2B), which may be used as biomarker metabolites for organ dysfunction. Notably, an additional 330 metabolites could only be found in two organs, potentially further expanding the set of potential organ dysfunction markers. Only 176 (10%) of the genes were core genes, which can be explained by the absence of the mitochondrially localized gene products from the core set, as the RBCs are lacking this organelle. When excluding the RBCs, an additional 142 genes were present in all remaining organs. The organ compendium represents a comprehensive set of manually curated, self‐consistent organ‐specific reconstructions that may be used for a range of biomedical applications (Aurich & Thiele, 2016; Nielsen, 2017).

The conversion of a reconstruction into a model involves the imposition of condition‐specific constraints (Orth et al, 2010). Each WBM reconstruction was constrained with 15 physiological parameters (Fig 3A) and metabolomics data (Wishart et al, 2013; Materials and Methods, Conversion from reconstruction to condition‐specific models, Dataset EV2: 3.2). For instance, metabolite transport across the capillaries of well‐perfused organs is bulk flow limited rather than diffusion‐limited (Feher, 2012). Using organ‐specific blood flow rates at rest (Price et al, 2003) and plasma metabolite concentrations (Wishart et al, 2013), we placed upper bounds on the corresponding organ‐specific metabolite transport reactions. Hence, by combining physiological and metabolomic constraints, we defined how much of each metabolite can be maximally taken up by an organ in the WBM, thereby addressing the challenge of appropriately constraining organ‐specific models. Additionally, we assumed that the kidney filters 20% of blood plasma along with all the metabolites, irrespective of their nature (Martini & Bartholomew, 2006). Accordingly, multiplying the glomerular filtration rate by healthy plasma metabolite concentration ranges (Wishart et al, 2013) allowed us to set lower and upper bounds for filtration of each metabolite by the kidney. Similarly, we assumed that cerebrospinal fluid (CSF) metabolites are unselectively returned into the bloodstream, permitting us to use healthy CSF metabolite concentration (Wishart et al, 2013) as constraints on the rate at which CSF metabolites enter the bloodstream. Using the daily urine excretion rate and urine metabolite concentrations (Wishart et al, 2013), we constrained the rate of urinary metabolite excretion (Fig 3A). We also applied constraints corresponding to an average European diet (Fig 3B; Elmadfa, 2012). All applied constraints are described in detail in the Dataset EV2: 3.2.1. In total, 12.5% of the WBM model reactions had a constraint placed on their bounds, leading to a substantially reduced steady‐state solution flux space (Appendix Fig S1A, Table EV12). We will refer to this novel paradigm in constraint‐based modeling as physiologically and stoichiometrically constrained modeling (PSCM) and provide a PSCM toolbox (opencobra.github.io/cobratoolbox) as well as a MATLAB (Mathworks, Inc.) Live Script (Dataset EV2, https://opencobra.github.io/cobratoolbox/stable/) enabling the reproducibility of all presented simulations.

To assess the predictive potential of the WBM models against current knowledge, we computed the metabolic essentiality of an organ by maximizing a whole‐body maintenance reaction, without the corresponding organ biomass contribution, under the aforementioned physiological and dietary constraints. We predicted an organ to be metabolically non‐essential if a non‐zero whole‐body maintenance reaction flux was possible in a WBM model deficient in that organ (Dataset EV2: 3.5). We defined an organ to be metabolically non‐essential in vivo if it can be surgically removed or replaced by an artificial organ (e.g., heart; Table EV15). Of the 32 organs, for 28 we could find corresponding information and predicted the metabolic essentiality correctly in all but one case (Fig 4A). It is important to note that the WBM models can only assess metabolic essentiality but do not capture other essential organ functions (e.g., those of the brain).

To further assess the predictive capacity of the WBM models, we simulated 57 inborn errors of metabolism (IEMs) by deleting the corresponding reaction(s) in all organs having the known defective gene(s) (Dataset EV2: 3.8, Materials and Methods, Statistical analyses). We then predicted and compared for each IEM the known biomarkers in either blood, urine, or CSF. Across the 57 IEMs, a total of 215 out of 252 (85.3%) biofluid specific biomarkers were qualitatively predicted correctly by Harvetta and 214 out of 252 (84.9%) by Harvey (Fig 4B, Table EV16). We compared this with similar predictions using the Recon3D model, which has only one extracellular compartment. As expected, the number of biomarkers was smaller than could be predicted (205 metabolites). Notably, the predictive potential of the Recon3D model was lower than that of the WBM models, with only 103 biomarker metabolites predicted correctly (50.2%, Fig 4B). Overall, the biomarker predictions with Harvetta and Harvey showed higher agreement percentages, increased Cohen's kappa metrics, and significantly better agreement percentages per disease (P = 0.0003 and P = 0.0001 for Harvetta and Harvey, respectively; Fig 4B). The better prediction capabilities of the WBM models were mainly caused by having fewer inconclusive or unchanged predictions, meaning that those metabolites were not predicted to be biomarkers. While for the Recon3D model, 45.3% (93 of 205) of the predictions failed to predict any changes, only 7.1% (18 of 252) fell into this category for Harvetta and 8.7% (22 of 252) for Harvey. Hence, the ability to account for physiological and metabolic constraints as well as organ‐level separation of metabolic pathways enhanced the prediction capacity considerably and presented a key advantage of using a WBM model over the generic Recon3D model.

In the WBM models, we investigated the constraints sufficient for the activity of two inter‐organ metabolic cycles (Frayn, 2010), when minimizing the Euclidean norm of the flux vector, with the rate of the whole‐body maintenance reaction set to one (Dataset EV2: 3.3.2.2). In the Cori cycle, glucose is secreted by the liver and taken up by the muscle, which produces lactate, which is then taken up by the liver to produce glucose (Fig 4C). Elevated postprandial blood glucose concentration, mediated by insulin levels, determines the rate at which glucose is taken up in vivo (Wasserman et al, 2011), but the regulatory functions of insulin are beyond the scope of the present WBM models. Consistently, only enforcing muscle glucose uptake was sufficient to achieve activity of the Cori cycle in the WBM models, at rest (Dataset EV2: 3.3.2.2). In the Cahill cycle, glucose is secreted by the liver and taken up by the muscle, which converts it to alanine, which is taken up by the liver to produce glucose (Fig 4C). The Cahill cycle is increased during exercise, when muscles degrade amino acids for energy needs, as a means to remove ammonium from the body via urea production in the liver (Frayn, 2010). Consistently, enforcing muscle glucose uptake, and enforcing either liver alanine uptake or increasing urinary urea excretion was sufficient to achieve activity of the Cahill cycle in the WBM models (Dataset EV2: 3.3.2.2).

To demonstrate that the aforementioned results are a consequence of including reactions and constraints based on organ‐level data, we then tested the degree of functional redundancy in the WBM models. Therefore, we removed random subsets of reactions and checking the feasibility of a unit flux through the whole‐body maintenance reaction. We found that there was an exponential drop in the number of feasible models as the fraction of reactions removed was increased (Appendix Fig S1B). Specifically, 60% of the WBM models (600/1000) were feasible if 0.1% of reactions (81/81,094) were randomly removed and none of the models was feasible when 1% of the reactions were randomly removed.

This reliable prediction of organ and gene essentiality, as well as the feasibility of inter‐organ metabolic cycles, and the infeasibility of whole‐body maintenance upon random reaction removal, demonstrates that the topology of the WBM models is a specific and sensitive representation of whole‐body metabolism but that constraints are required to enforce condition‐dependent metabolic activity.

The basal metabolic rate (BMR) varies between individuals, depending on their age, sex, weight, height, and body composition. Phenomenological methods (Harris & Benedict, 1918; Mifflin et al, 1990) have been developed for estimating BMRs, but they do not account for differences in metabolism. Here, we investigated whether the WBM models, which incorporate metabolism and physiological constraints (Fig 3), could be used to predict BMRs with higher accuracy than the phenomenological methods. To this end, we assumed that BMR is equal to the energy yield from ATP hydrolysis (64 kJ; Wackerhage et al, 1998) times the resting rate of steady‐state ATP consumption. The latter was predicted by setting the whole‐body maintenance reaction flux to one, minimizing the Euclidean norm to obtain a unique steady‐state basal metabolic flux vector (Heirendt et al, 2019) and summing the rate of ATP consumption of all reactions in all organs (Dataset EV2: 3.3). The whole‐body maintenance reaction represents the material and energy (ATP) required to maintain the non‐metabolic cellular functions of the body and was constructed based on the fractional weight contribution for each organ for a reference man or woman (Snyder et al, 1975) (Dataset EV2: 3.2). In the basal metabolic flux vector, ATP consumption by the whole‐body maintenance reaction is < 2% of the total ATP consumption. Without any fitting, or accounting for age, we predicted overall energy consumption rates of 1,344 kcal for Harvetta and 1,455 kcal for Harvey, which were consistent with the expected values for a reference man and woman, respectively.

Age, weight, and height from 13 healthy women with normal body mass index (BMI < 22.1 ± 2.4) (Prentice et al, 1986) (Fig 5A, Dataset EV2: 3.9, Table EV17) were used, together with established polynomial functions (Dataset EV2: 3.10.1.5; Young et al, 2009), to fix the physiological parameters (organ weights, blood flow rate, and cardiac output) of a set of personalized Harvetta models. We then tuned Harvetta to fit corresponding personalized BMR measurements by performing a sensitivity analysis (Dataset EV2: 3.9.1) thereby identifying parameters in Harvetta that could be adjusted to fit the BMR data better. As the fat‐free body mass was provided for each individual, we adjusted the muscle and adipocyte coefficients in the whole‐body maintenance reactions, instead of solely relying on the estimated values based on age, weight, and height (Young et al, 2009). We also increased the ATP coefficient in the muscle biomass maintenance reaction, starting from a previously estimated coefficient (Bordbar et al, 2011a), to better reflect the energetic requirements associated with protein turnover in the muscle. After optimizing the ATP coefficient, the WBM model was able to predict the measured BMR with high accuracy (explained variance = 79.7%, F(1, 11) = 43.14, P = 4.04e‐05, Fig 5B) and better than the Mifflin‐St Jeor equations (Fig 5C). However, as parameter tuning was involved, potentially leading to overfitting, we validated the WBM BMR prediction models with an independent data set not utilized in parameter estimation.

An independent, external data set, composed of age, weight, height, fat‐free body mass, and BMR data from six female and 16 male athletes (Loureiro et al, 2015) (Fig 5D) was used to assess the ability of the previously tuned WBM models to predict BMR (Dataset EV2: 3.9.2, Table EV17). Good correlation was obtained between predicted and measured BMR (explained variance = 61.3%, F(2, 19) = 15.06, P = 1.10e‐4) (Fig 5E). In comparison, the Mifflin‐St Jeor equations, a standard phenomenological method of estimating the BMR from age, sex, height, and weight (Mifflin et al, 1990), gave a lower correlation (R² = 0.49, P = 0.002) (Fig 5F). Furthermore, given the WBM models, the Mifflin‐St Jeor equations were not informative for the true BMR (linear regression b = 0.01, 95%‐CI:(−0.99, 1.00), P = 0.90). In contrast, the WBM prediction added predictive value to the Mifflin‐St Jeor equation (linear regression b = 2.24, 95%‐CI: (0.32; 4.16), P = 0.025). The last result implies that the WBM models are significantly superior to the Mifflin‐St Jeor equations in terms of prediction accuracy. Taken together, the results demonstrate that the WBM models provide good and replicable estimations of BMRs that improve on standard phenomenological methods. However, the analyses also highlight the need for personalized parameterization of the WBM models for a given simulation and that the use of further parameters is likely to improve the prediction fit. Equally, confounding factors, beyond age and sex, should be considered and validation in a larger cohort will be required.

Beyond the basal metabolic rate, activity‐induced energy expenditure consumes the remaining energy and nutrient resources provided by the diet (Biesalski & Grimm, 2005). With the aforementioned basal metabolic flux distribution (Dataset EV2: 3.3.2.3), the resting brain energy consumption rate was initially predicted to be 10.5 mole ATP/female/day (Dataset EV2: 3.3.2.3). Assuming a conversion of 1 mole of glucose into 31 moles of ATP, this value is lower than a total brain consumption rate of 20.65 mole ATP/day/person, based on the reported consumption rate of 120 g of glucose per day (Berg et al, 2002). To achieve a 20 mole ATP/female/day total brain consumption rate, we had to set a requirement for the brain to hydrolyze 3.5 mole ATP/female/day for electrophysiological processes (VMH ID: Brain_DM_atp_c_). Note that this demand for hydrolysis of ATP requires further hydrolysis of ATP by supporting metabolic pathways elsewhere in the brain. Similarly, the pumping of blood through the body by the heart requires the hydrolysis of ATP. However, literature reports differ in their estimated energy requirements with one source reporting 6,000 g ATP/day/person to be consumed (https://heartmdinstitute.com/heart-health/metabolic-cardiology-basics/), which corresponds to 11.8 mole ATP/day/person with the molecular weight of ATP being 507.18 g/mole (Dataset EV2: 3.2.2.7). A lower bound on a heart demand for hydrolysis of ATP (VMH ID: Heart_DM_atp_c_) was set to 6 mole ATP/day/person to account for the energy required for pumping blood through the body while allowing the remained to be consumed by supporting heart metabolic reactions. Despite the lack of comparable measurements of activity‐induced ATP consumption rates, we set these parameters on brain and heart ATP hydrolysis because each organ contributes to the overall energy consumption in a whole‐body model. After accounting for whole‐body maintenance and energy demand due to brain and heart activity, the basal metabolic flux distribution predicted that the remaining energy that the WBM models could extract from the dietary input was 515.7 kcal for Harvetta and 502.5 kcal for Harvey. These values correspond to approximately 20% of the average European dietary input (2,372 kcal) used for the WBM simulations. This result is consistent with a 15–30% allocation of dietary input toward physical activity (Biesalski & Grimm, 2005).

Overall, we demonstrated that WBM model predictions of basal metabolism and energy requirements are overall consistent with the literature and they provide an opportunity for integration of both physiological and molecular biological data to enable enhanced predictions when compared with purely phenomenological methods.

The human gut microbiota influences host organ metabolism, including the brain (“gut‐brain axis”), the liver (“gut‐liver‐axis”), and the colon (Clarke et al, 2014). However, our understanding of host–microbiome co‐metabolism remains limited. The WBM models can be used to predict how the microbiome may modulate human metabolism, thereby complementing existing in vitro and in vivo analysis approaches.

To constrain Harvetta and Harvey in a gut microbial context, we used published shotgun metagenomics and clinical data on 149 healthy individuals provided by the Human Microbiome Project (HMP) Consortium (Consortium, 2012a; Peterson et al, 2009) (Dataset EV2: 3.10, Fig 6A). First, we generated personalized germ‐free WBM models using clinical data on sex, weight, height, and heart rate of each individual. Second, we mapped the reported strain‐level abundances for each individual onto a resource of 773 gut microbial metabolic reconstructions (Magnusdottir et al, 2017; Dataset EV2: 3.10.1). The resulting 149 personalized microbial community models captured 131 ± 19 microbes per HMP individual, accounting for 91 ± 7% of the reported relative abundance (Appendix Fig S2). Third, the microbial community models were added to the large‐intestinal lumen of the corresponding personalized WBM models, resulting in personalized microbiome‐associated (my) WBM models. In the large‐intestinal lumen of these personalized myWBM models, gut microbes can take up dietary metabolites after they transit from the small intestinal lumen where small intestinal enterocytes may also absorb them. Metabolites produced by gut microbes can be taken up by colonocytes or excreted via the feces.

The human gut microbiome is a source of amino acids and neuroactive compounds (Clarke et al, 2014; Lin et al, 2017). Hence, we investigated whether the presence of the gut microbiota in personalized myWBM models could increase the maximal possible production rate of eight prominent neurotransmitters. We first maximized the neurotransmitter production in the brain in the personalized germ‐free WBM models and observed substantial inter‐individual variation (Fig 6B). This variation could be attributed to the changed bounds on the organ uptake reactions due to the personalization of the physiological parameters. We then used the personalized myWBM and observed substantial increase in the brain neurotransmitter production capability. No clear correlations between brain neurotransmitter synthesis potential and individual microbial species were observed (Fig 6C). The high inter‐individual variation in these metabolic objectives could neither be explained by the presence of a single species or genus (Fig 6B and C) nor by the provided meta‐data. Consequently, we hypothesized that the simulation results were the direct consequence of the host–gut microbial co‐metabolism and could not be derived from the metagenomic data alone. Using the 149 personalized microbiome models without the WBM, we calculated the neurotransmitters and neurotransmitter precursors predicted to be secreted by each strain (Appendix Fig S3). For instance, GABA was predicted to be produced by Akkermansia muciniphila and strains of the Bacteroides, Escherichia, and Shigella genera. Consistently, GABA‐producing capability has been shown for Bacteroides, Parabacteroides, and Escherichia species (Strandwitz et al, 2019). Glutamate, the precursor of GABA and itself a neurotransmitter, was predicted to be mainly produced by Eubacterium rectale, Faecalibacterium, Roseburia, and Bacteroides. Histamine was predicted to be produced mainly by Gammaproteobacteria species and Streptococcus sp., in agreement with reports that members of these taxa synthesize histamine (Clarke et al, 2014). Histidine, the precursor of histamine, was predicted to be produced mainly by commensals of the Bifidobacterium genus and of the Clostridiales order (e.g., Faecalibacterium, Roseburia). Alistipes and Bacteroides sp. were predicted to produce the majority of L‐tryptophan, a precursor of serotonin. Finally, L‐tyrosine, a precursor of dopamine, was predicted to be produced by a variety of species belonging to the Actinobacteria, Firmicutes, and Proteobacteria phyla. Taken together, the modeling results predicted that gut microbes have the molecular mechanisms that may influence host brain metabolism directly by synthesizing neurotransmitters and indirectly by providing the amino acid precursors to be further metabolized by human enzymes. The predicted maximum secretion potential, as well as the contributing strains, depended on each individual's microbiome composition. Thus, the WBM models can be used to generate novel mechanistic hypothesis as to how the gut microbiota might influence brain metabolism (Clarke et al, 2014; Lin et al, 2017).

Ethanol causes liver toxicity as its product, acetaldehyde, is cytotoxic and carcinogenic (Neuman et al, 2017). When maximizing the flux through the liver alcohol dehydrogenase (VMH ID: ALCD2if), the myWBM models predict an 11.4 ± 8.5 fold increase in maximal potential flux compared with corresponding personalized germ‐free WBM models (Fig 6E). This fold change in flux strongly correlated with species belonging to the Clostridia class and negatively correlated with representatives of the Bacteroidia class (Fig 6E). Consequently, the Bacteroidia/Clostridia ratio determined 92.1% of the variation in the predicted maximum alcohol dehydrogenase flux (F(2, 144) = 852.21, P = 1.56e‐80, see Table EV18). Species that affected the most the alcohol dehydrogenase flux included Clostridium scindens, Blautia hydrogenotrophica, and the well‐known pathobiont Clostridium difficile (Fig 6C). Endogenous ethanol production by the microbiota has been shown in non‐alcoholic steatohepatitis (NASH) patients (Zhu et al, 2013). Further support that microbial endogenous ethanol production can perturb liver function and cause non‐alcoholic fatty liver disease has been recently published (Yuan et al, 2019).

Microbially produced butyrate serves as the main energy source for the colonocyte (Donohoe et al, 2011). We computed the maximum flux through the butyrate‐CoA ligase in the colonocytes in the personalized myWBM models. The predicted maximum butyrate‐CoA ligase flux was on average 78.7 ± 64.7‐fold increased in the presence of the microbiome (Fig 6F) and strongly correlated with the abundance of known butyrate producers, such as Faecalibacterium prausnitzii, Butyrivibrio crossotus, and Subdoligranulum variabile (Louis & Flint, 2009). The Bacteroidia/Clostridia ratio determined 92.9% of the variance in the predicted maximum butyrate‐CoA ligase flux (F(2, 144) = 953.32, P = 8.8e‐84, Table EV18). Consistently, it is known that butyrate is mainly produced by microbes belonging to the Clostridia class (Louis & Flint, 2009).

The human gut microbiota can also influence drug metabolism (Spanogiannopoulos et al, 2016). Flux balance analysis was used to predict the maximal possible flux through the liver sulphotransferase reaction (VMH ID: PCSF) in each myWBM model (Fig 6D). myWBM models with low Clostridioides genus abundance were predicted to have lower maximum liver sulfonation flux (Fig 6D). The potential for liver sulfonation was inversely correlated with the abundance of the Bacteroidia class (Fig 6D), such that the Bacteroidia/Clostridia ratio contributed 53.6% of variance to maximal liver sulfonation rates (F(2, 144) = 85.37, P = 3.547e‐25, Table EV18). Note that, at higher Clostridioides genus abundances, substrate availability was limiting for p‐cresol production and subsequent sulfonation (Appendix Fig S4). Overall, consistent with our results, it is known that the production of p‐cresol, e.g., by Clostridium difficile, competes with drugs, such as acetaminophen, for sulfonation in the liver and interferes with drug detoxification (Spanogiannopoulos et al, 2016).

Taken together, using the personalized myWBM models, we demonstrate that this modeling framework can generate predictions of how microbial metabolism may modulate host metabolism. Importantly, some of these predictions have already been shown experimentally to be consistent with known host–microbe co‐metabolism. Thus, the myWBMs are suitable for the generation of novel, experimentally testable hypotheses on host‐microbe co‐metabolism, the development of microbe‐modulated biomarkers, and for the identification of individuals with increased risk of toxicity.

In this part, we describe the reconstruction approach and the information used to generate the WBM reconstructions, Harvey and Harvetta. The conversion of the reconstructions into condition‐specific WBM metabolic models is then described Statistical analyses details can be found at the end of the Materials and Methods section. Simulation‐relevant information and methods can be found in the Dataset EV2 (or here: https://opencobra.github.io/cobratoolbox/stable/).

To ensure reproducibility of the presented work, and also to reduce the number of supplemental tables, we provide a MATLAB (Mathworks, Inc) live script along with a PSCM toolbox, which is an extension to the COBRA Toolbox (Heirendt et al, 2019) (to be obtained from https://github.com/opencobra/cobratoolbox).

Please refer to the accompanying Matlab LiveScript and its html version. Dataset EV2. The PSCM Toolbox contains all information, data, scripts, and functions to repeat the simulations in this paper Dataset EV1.