Understanding the host-microbe interactions using metabolic modeling

Abstract

The human gut harbors an enormous number of symbiotic microbes, which is vital for human health. However, interactions within the complex microbiota community and between the microbiota and its host are challenging to elucidate, limiting development in the treatment for a variety of diseases associated with microbiota dysbiosis. Using in silico simulation methods based on flux balance analysis, those interactions can be better investigated. Flux balance analysis uses an annotated genome-scale reconstruction of a metabolic network to determine the distribution of metabolic fluxes that represent the complete metabolism of a bacterium in a certain metabolic environment such as the gut. Simulation of a set of bacterial species in a shared metabolic environment can enable the study of the effect of numerous perturbations, such as dietary changes or addition of a probiotic species in a personalized manner. This review aims to introduce to experimental biologists the possible applications of flux balance analysis in the host-microbiota interaction field and discusses its potential use to improve human health.

Video abstract

Application of bacterial metabolic networks to study metabolic interactions

The gut microbiota is the community of microorganisms residing in the gut and include commensal, symbiotic, and pathogenic bacteria. Under normal circumstances, the gut microbiota and its host are in symbiosis [1]. Disruption of the symbiosis is detrimental for host health and can result in disease including gastrointestinal disorders such as inflammatory bowel disease [2], metabolic disorders such as diabetes mellitus [3], and mental disorders such as autism spectrum disorder [4], and major depressive disorder [5]. To understand the symbiotic relationship, the different members of the gut microbiota, and the way they communicate with each other and with the host need to be known. The gut microbiota communicates via the production of metabolites [6]. Therefore, it is key in the field of host-microbe interactions to identify which microbial members are present and what their metabolic output is. However, this does not fully elucidate the dynamic interactions within the microbiota, and between the host and the microbiota, since the metabolic output of microorganisms is dependent on their surroundings [7]. Therefore, the metabolic output and, in turn, the interactions between the host and the microbiota is different among individuals [8], making successful treatment of the aforementioned disorders more challenging. Although experimental approaches are crucial to the progress of the microbiota field, they are not able to fully capture the mechanisms, interactions, and behavior due to the huge complexity of the gut environment. These limitations have led to the development of a complimentary approach to completely understand the relationship between the host and its microbes; bacterial metabolic networks [9]. In this approach, bacterial interactions can be visualized in the form of a metabolic network. The metabolites comprise the nodes of the metabolic network. Biological processes such as conversions, uptake, and secretion are represented by the edges. By placing the set of metabolic reactions of a single bacterium into a compartment, one could separate the metabolic reactions of one bacterium from another bacterium or from host cells in the metabolic network. Placing each cell compartment into a shared compartment can simulate how the different cells metabolically interact with each other [10] (Fig. 1a).

Fig. 1

Application of flux balance analysis to simulate a metabolic interaction among multiple bacteria. a A metabolic network of a bacterial community consisting of Faecalibacterium prausntizii and Bifidobacterium adolescentis (adapted after El-Semman et al. [20]). Fluxes of the exchange reactions are represented by arrows. Solid black arrows indicate uptake, and secretion reactions of the bacteria, dashed black arrows indicate the flow of metabolites in or out of the system, and dashed grey arrows indicate the formation of new biomass, where metabolites cannot be secreted by the bacteria anymore, thus leaving the system. b Stoichiometric matrix of the exchange reactions depicted in panel a. BA depicts B. adolescentis and FP depicts F. prausnitzii. c Visual representation of the concept of flux balance analysis; (i) A solution space of the flux distribution in a system; (ii) The allowable solution space is the result of addition of constraints depicted in Eqs. 1 and 2, which equations limit the available flux distributions; (iii) Addition of an objective function as depicted in Eq. 3 determines the optimal solution for the flux distribution in the system. Adapted after Orth et al. [3].

Metabolic networks can be used to predict the effect of alterations in the metabolic network in silico. A constraint-based reconstruction and analysis (COBRA) approach is often used to simulate the operation of a metabolic network under different external nutrient conditions. A COBRA method suitable for investigating the metabolism of the microbiota is flux balance analysis (FBA) [11]. A flux is the rate of turnover of a metabolite through a metabolic pathway. To perform FBA, all the fluxes in the network should be represented by a set of linear equations. The equations are placed in a stoichiometric matrix, which consists of the substrates, products, and directionality of the reactions (Fig. 1b). Next, FBA uses constraints to limit the flow of metabolites through the network and calculates the distribution of metabolic fluxes in the metabolic network for a given objective function (OF), resulting in an optimal distribution of fluxes (Fig. 1c). For example, an OF that maximizes the production of the short chain fatty acid (SCFA) acetate will result in a different flux distribution compared to an OF that maximizes butyrate production. Other examples of OFs include minimizing the production of a particular metabolite, maximizing cell growth or for a community of bacteria maximizing growth of a single species or maximizing the growth of the whole community [11].

FBA works at steady-state, i.e., the amount of the metabolite produced is equal to the amount of the metabolite consumed. Accordingly, the set of linear equations is formulated as:

$$ \mathrm{S}\kern0.5em \ast \kern0.5em \mathrm{v}\kern0.5em =\kern0.5em 0 $$

(1)

Where (S) represents a stoichiometric matrix, and (v) represents the flux distribution. The formula describing a biomass function contains all the metabolites in the system that are required to build a new cell. In other words, a biomass function simulates cell growth. Importantly, since Eq. 1 depicts a set of linear equations and there are typically more reactions than compounds, there is more than one flux distribution possible. Constraints are imposed on a flux as upper- and lower-bounds to limit the maximum and minimum values that each flux can take. Constraints could reflect media conditions, where uptake and secretions rates are limited, or reaction velocities of internal enzymes obtained from experimental data [12, 13]. Accordingly, each flux with constraints is formulated as:

$$ {\mathrm{V}}_{\mathrm{i},\min}\kern0.5em <\kern0.5em {\mathrm{V}}_{\mathrm{i}}\kern0.5em <\kern0.5em {\mathrm{V}}_{\mathrm{i},\max } $$

(2)

Similarly, the OF is formulated as follows:

$$ \mathrm{Z}\kern0.5em =\kern0.5em {\mathrm{c}}^{\mathrm{T}}\ast \kern0.5em \mathrm{v} $$

(3)

where (Z) is the solution of the OF, (C) is the vector of weights, indicating how much each reaction contributes to the OF and (T) represents the matrix transpose. For example, when a single flux is maximized or minimized, C is a vector of 0s with a single 1 [14, 15]. However, formulating an OF can be challenging and is fully dependent on the research question [16, 17].

FBA is a versatile tool to employ for many purposes. By adjusting the upper- and lower- bounds of metabolites, growth on different media [18] or, in the case of the gut microbiota, changes in diet, can be simulated. Similarly, by setting the flux of a certain metabolite to zero, a gene knock-out or the absence of a member of the microbiota can be simulated [11, 15, 19]. In this way, we can estimate the viability of a microbial community under different conditions as well as the effect of adding new species to a bacterial community on host health [16]. For example: FBA used to study the effect of lactate production on Bifidobacterium adolescentis shows a reduction in the production of formate, ethanol and acetate as well as a reduction in biomass formation if lactate production is manually increased. The OF in this example is maximizing biomass production (Fig. 2) showing that the flux distribution changes with altering the environment [20]. To investigate interactions in a microbial community, the authors added another bacterium, Faecalibacterium prausnitzii, in the same metabolic environment as B. adolescentis. F. prausnitzii needs acetate to grow well on glucose and to produce butyrate [21]. If acetate is not supplied, F. prausnitzii will use the acetate produced by B. adolescentis to grow and consequently produce butyrate (Fig. 3). Since the biomass reaction simulates bacterial growth, altering the flux through the biomass reaction of one bacterium, while keeping the total biomass of the system constant, can simulate changes in the bacterial composition in the gut. Altering the amount of flux through a biomass reaction will result in alteration of the flux distribution of the whole system, which, in turn, causes an increase in butyrate production when the flux through the biomass reaction of F. prausntizii is increased in the system compared to a situation with more flux through the biomass reaction of B. adolescentis. The OF in this example is minimization of glucose consumption for both bacteria (Fig. 3) [20]. This example shows that FBA can be used to investigate interactions between bacteria. In a similar approach, the introduction of a compartment that represents a host cell and connecting it to a shared metabolic compartment that represents the intestinal lumen can be used to study host-microbiota interactions (Fig. 4). Using this approach, Heinken et al. combined metabolic compartments of Bacteroides thetaiotamicron and a generalized mouse cell [22] in a metabolic network. The authors investigated metabolic dependencies between the host and its microbe by simultaneous optimization of the growth rates of both the host cell and the microbe. The authors showed that the presence of B. thetaiotamicron could influence the growth of the generalized mouse cell by supplying the host cell with essential and non-essential amino acids. Furthermore, the authors simulated gene knockouts of B. thetaiotamicron and the mouse cell by setting the corresponding flux to 0. By optimizing the biomass function of the cell with the knockout the authors showed that B. thetaiotamicron is able to rescue lethal knockouts in the host cell and vice versa [23]. Although this metabolic network shows only an interaction between one bacterium and one host cell type, which is far from accurately depicting the community in the gut, expanding this network by introducing more bacterial and host compartments can give new insights into systemic effects of the microbiota on the host and vice versa. Metabolic networks depicting human metabolism, which includes host interactions with the microbiota can be used to investigate causality in microbiota-related disease and, in turn, can generate new hypotheses for the treatment of microbiota related diseases. To date, the most comprehensive model of human metabolism is Recon3D [24, 25].

Fig. 2

Representation of the use of FBA in metabolic modeling of an organism. The OF in panels a and b is maximizing the production of the biomass. The thickness of the arrows indicate the amount of flux, where a thicker arrow indicates a higher flux. When the production of lactate is manually altered in panel b, the flux distribution changes, whereby the flux of acetate, formate, ethanol, and biomass is lowered compared to panel a. Adapted after El-Semman et al. [20]

Fig. 3

Representation of the use of FBA in the modeling of a bacterial community consisting of two gut bacteria: F. prausnitzii and B. adolescentis. Fluxes are represented by arrows. Solid black arrows indicate uptake and secretion reactions of the bacteria, dashed black arrows indicate the flow of metabolites in and out of the system and dashed grey arrows indicate the formation of new biomass, where metabolites can no longer be secreted by the bacteria, thus leaving the system. The amount of flux is represented by the thickness of the arrows, a higher flux is a thicker arrow. The ratio of the produced biomass is used as a measure for the number of F. prausnitzii and B. adolescentis in the system, whereby the total biomass is kept constant. A higher flux through the biomass reaction represents more bacteria of that species in the system. The OF is the minimization of glucose uptake for both bacteria. Adapted after El-Semman et al. [20]

Fig. 4

Representation of the use of FBA to study the complex host-microbe interactions. A bacterial metabolic compartment is placed in a compartment which is connected to a metabolic compartment representing a host cell. The host cell compartment is connected to another compartment representing the bloodstream of the host. Arrows between the different compartments indicate exchange reactions. Solid arrows represent influx of metabolites into the system, which represents metabolites originating from the diet and/or metabolites present in the bloodstream. Dashed arrows represent efflux of metabolites out of the system symbolizing metabolites excreted in feces and/or translocating in the bloodstream

The examples of metabolic networks depicted in Figs. 2, 3, and 4 show the usefulness of computational methods in understanding the complex interactions within the gut microbiota. Fortunately, the principles of using FBA in small metabolic networks are the same as in large metabolic networks. Thus, a thorough understanding of these principles will help making predictions beyond what experimental biologists can expect from logical reasoning to eventually allow the discovery of novel interventions to improve human health.

Since FBA works at steady-state, it gives an optimal distribution of fluxes for a given OF in small metabolic networks. In large metabolic networks multiple solutions are possible. To find these alternative solutions flux variance analysis can be used, which uses FBA to maximize and minimize each flux in the system [15]. Finding optimal solutions by simulating a metabolic network once does not capture dynamic changes in metabolite levels and dynamic interactions between bacteria. Dynamic changes are an important factor when studying bacterial communities, because cells are dividing and dying. Additionally, in the human gut there is movement of cells from the upper part of the intestine to the lower part due to gut motility and metabolites levels change dynamically during the day due to the cycles of food intake [26,27,28]. Therefore, capturing the dynamic changes in metabolic networks is necessary to study the gut microbiota. To capture dynamic changes, dynamic FBA (DFBA) has been developed. DFBA uses ordinary differential equations (ODE) to couple FBA to a kinetic model [29]. This can be done using three approaches: static optimization approach (SOA), dynamic optimization approach (DOA), or direct approach (DA) [30]. The most widely used approach is SOA. Essentially, SOA makes a series of snapshots using FBA. The starting conditions of each snapshot is determined by the outcome of the previous snapshot. SOA requires small time steps between the snapshots to accurately capture the dynamic changes, making it computationally expensive [31]. DOA obtains time profiles of fluxes and metabolite levels by optimizing the entire time period of the simulation. The dynamic optimization problem is transformed into a non-linear programming problem which is solved once [31, 32]. In contrast to DOA, DA uses a linear program solver at the right side of the ODE. Similar to SOA, DA is also computationally expensive since the linear program needs to be solved each time the right side of the ODE is evaluated [30]. An example of the application of the SOA was shown by Mahadevan et al. to investigate the restructuring of the metabolic network of Escherichia coli during a diauxic shift. The authors simulated a batch culture of 10 h, which was divided into 10,000-time intervals. The concentration of the metabolites at the start of each interval was directly calculated from the previous interval. This approach showed that during the first 4.6 h the oxygen, and glucose uptake rates of the cells were the limiting constraints on biomass formation. During the second phase between 4.6 h and 6.9 h, the limiting constraint on biomass formation was the concentration of oxygen in the environment. During the last phase, where acetate is utilized from 6.9 h to 10 h the mass transfer coefficient was the limiting constraint on biomass formation. The computational values were similar to experimental values confirming the usefulness of DFBA to determine the limiting factors of growth in bacteria [32]. DFBA has been used to study communities of microorganisms including a community of E. coli and Saccharomyces cerevisiae [33]. However, DFBA can only be used in small communities of bacteria because adding more species in a network increases the amount of reactions dramatically, which consequently increases the time and costs for simulations. DFBA in microbial research is extensively reviewed elsewhere [30, 34]. To simulate microbial communities with FBA, not only the uptake and secretion reactions are necessary, but all the molecular capabilities of each bacterium in the metabolic network should be known. Thus, a metabolic network of each individual bacterium is needed. This is done by generating metabolic models of each bacterium directly from their sequenced genome.

Constructing metabolic networks from annotated genomes

To generate and simulate a metabolic network, the metabolic capacity of the bacteria in the network should be known and translated into a model that can be used by FBA. Genomic-scale metabolic models (GEMs), also known as genomic-scale metabolic reconstructions (GENREs), consist of all the metabolic capacities of a bacterium [35], and are automatically constructed from annotated genomes [36]. Several tools are available to construct GEMs, which are reviewed elsewhere [37]. However, not all genes of an organism are active during each growth phase, or in every environment [38, 39]; hence, the constraints of these automatically generated GEMs should be manually refined [40]. Recently, AGORA, a semi-automated GEM database of 818 members of the gut microbiota became available [41]. However, the AGORA GEMs differ from other available GEM databases such as BiGG, KBase, and CarveMe [42,43,44]. Taxonomically, the AGORA database contains a wider variety of organisms and the GEMs are more generally constructed, making it possible to use these GEMs in different ways. Furthermore, the AGORA database is constructed to simulate the gut microbiota community in general. Therefore, the AGORA GEMs have, for example, more carbon-uptake reactions compared to the BiGG database, which has a more limited number of bacteria and is generally used to simulate metabolic capabilities of a single bacterium and predicts possible changes in the metabolic capacity of this bacterium in case it harbors inoperative genes [45, 46]. For example, the BiGG database contains numerous GEMs of E. coli, which is a heavily researched organism, thus the metabolic reactions and constraints can be accurately determined from literature. On the contrary, the AGORA database contains numerous underinvestigated organisms, which results in more general GEMs, highlighting the need for further investigation, update of the metabolism of gut microbes, and refining of the GEMs accordingly [46]. The AGORA database should thus be seen as repository for general GEMs of gut bacteria, which should be made into condition-specific GEMs before usage by adding proper constraints and formulating OFs depending on the research question of the user. This needs to be done before using GEMs obtained from databases containing condition-specific GEMs as well [47]. For example, Pryor et al. showed an in vivo increase in the lifespan of Caenorhabditis elegans due to bacterial-produced agmatine, whereby agmatine production is induced by metformin a commonly used drug in type 2 diabetes patients. The authors used GEMs from the AGORA database to translate the observed in vivo results found in C. elegans to humans. The authors constructed a metabolic network based on 16S sequencing data from metformin-treated and untreated type 2 diabetic patients, respectively. The authors successfully showed a molecular mechanism between the host and its microbes via their findings of a higher production of agmatine in the metformin-treated patients. However, before the authors could use the GEMs from the AGORA database, they manually curated the GEMs and included agmatine uptake and secretion reactions. Moreover, the authors had to adjust constraints of the AGORA GEMs, which were originally constructed using a western diet [48] to the dietary information available for the patients in the cohort used for the construction of metabolic communities. Taken together, the metformin study highlights the need the carefully optimize the AGORA models specific for the research question of the user [19, 48,49,50].

Nevertheless, community efforts are being taken to standardize GEMs and assure the quality of the models [45]. To optimize GEMs of bacteria, Kuang et al. combined FBA with an untargeted mass spectrometry-based approach to identify metabolites produced by Citrobacter sedlakii, a non-pathogenic bacterium found in human stool. The authors investigated the metabolic output of C. sedlakii by taking samples at different growth stages and analyzing the extracts using two liquid chromatography mass spectrometry (LCMS) approaches, reverse phase (RP) and hydrophilic interaction liquid chromatography (HILIC). The obtained data was compared to a predicted list of metabolites generated via FBA using a tool called MS_FBA. The comparison showed that the metabolic output of the GEM of C. sedlakii did not cover all the metabolites measured with LCMS, highlighting the need to determine the accuracy of bacterial GEMs and more precise genome annotation [51]. Lastly, when comparing flux distribution between different sets of microbial profiles, for example between patients and healthy controls, the metabolic environment of the microbial community should be known. This can be partly inferred from the diet of the patients, but never to a full extent. Another approach to infer the metabolic environment of the microbial community is to infer the metabolic environment from abundance distribution in individual samples. The Metabolic Analysis of Metagenomes using FBA and Optimization (MAMBO) approach combines 16S sequencing data from individual fecal samples with GEMs obtained from reference genomes of the sequenced bacteria to infer the metabolome in the sequenced fecal sample. The rationale behind this approach is that the metabolic environment shapes the abundance profile in a given sample, because the metabolic potential of the bacteria indicate which species will thrive in the metabolic environment [52]. To answer this question, GEMs were constructed from annotated reference genomes made available in the human microbiome project [53] using the modelSEED pipeline [54]. In addition, 372 GEMs were obtained from the AGORA database. The constructed GEMs were added in the same metabolic environment in a similar way as seen in Fig 3. The starting concentration of the metabolites in the metabolic environment was defined at random and FBA was performed, where the OF is defined as: the biomass functions of all the GEMs should be as close to the obtained abundance distribution. In a next step the metabolites in the metabolic environment were slightly altered and FBA with the same OF was performed again. If the Pearson correlation between the biomass functions and the abundance profile is higher compared to the previous metabolic environment the metabolite change is accepted and the process is repeated until the Pearson correlation does not improve. Comparing the biomass distribution in the FBA results with the abundance profile obtained from the fecal sample shows which metabolic environment fits the abundance distribution best [55]. Collectively, the application of a well-defined metabolic environment, microbial composition and metabolic capabilities, the relationship between host and microbe can be studied.

Application of GEMs in modeling the gut microbiota

To understand the effect of the microbiota on the host, microbial metabolic networks can be extrapolated to include metabolic networks of the host. This enables predicting the effects of gut microbiota on the host and suggesting possible interventions to promote host health. Diener et al. constructed microbial metabolic networks from the metagenome data of a cohort of Swedish people with and without diabetes mellitus and used FBA to analyze SCFA production. In this study, 16S abundance data was integrated with metabolic fluxes to construct personalized predictions of metabolic output and interventions such as dietary changes and medical treatment. The authors used GEMs from the AGORA database, however the AGORA database did not cover all the species found in the Swedish cohort. Instead, the authors performed FBA on the genus and species level, whereby AGORA models were pooled together into higher phylogenetic ranks. This resulted in networks containing between 12 and 30 GEMs at the genus level and 23 and 81 GEMs at the species level. The GEMs were placed in a shared metabolic environment for each sample in the cohort representing the gut lumen. The relative biomass for each GEMs was estimated from the relative read distribution in each sequenced sample. The authors found a minimal overlap of resource utilization between microbes in different niches, suggesting an upper bound on alpha diversity in the gut. Next to ecological insights, the authors concluded from their model that SCFA production is highly specific per individual. Nonetheless production of butyrate and propionate was reduced in diabetic subjects compared to healthy controls and the overall SCFA production profile could be restored upon metformin treatment [19]. The authors did not measure metabolite levels in vitro, but the outcome is in line with other experimental-based research [56]. Another study expanded GEMs found in the AGORA database with reactions for bile acid metabolism. Here, paired bacteria were placed in a shared metabolic environment to compare their ability to metabolize bile acid with that of a single bacterium. An average European diet supplemented with the bile acids taurocholate, glycocholate, taurochenodeoxycholate, and glycochenodeoxycholate was used as modeling constraints and maximizing the exchange of the bile acids was used as the OF. The authors showed that metabolism of bile acids is higher in a microbial community compared to a single bacterium because individual bacteria cannot metabolize each bile acid. In a similar manner, the authors integrated publicly available microbial abundance data of healthy individuals and patients with inflammatory bowel disorders (IBD) into metabolic network to study bile acids metabolism in IBD. The study showed that the metabolism of bile acids was lower in patients with IBD [49], which is in line with in vivo results [57]. However, the authors did not perform interventions such as dietary change or introduction of other bacteria in the metabolic network to predict possible treatment of IBD patients. Dietary interventions were investigated in another study that used 16S sequencing data to construct metabolic networks of 28 Crohn’s disease (CD) patients, and 26 healthy controls. The study placed GEMs obtained from the AGORA database and used BacArena, a tool used for modeling bacterial communities [58], to place the GEMs in a grid environment for each individual (Fig. 5). The bacterial composition obtained from the 16S sequencing data was used to determine the number of GEMs for each species placed in the grid environment at the start of the simulation. The microbial biomass was used as the OF and all possible metabolites that could be taken up by the GEMs were added in high concentrations to the environment. Bacterial growth was simulated in 24 steps of 1 hour. The last time-point was used to compare the microbial abundance and metabolite production between CD patients and healthy controls. Interestingly, the outcome showed a higher production of SCFAs in the control group. Next, the authors identified metabolites, which could have resulted in increased SCFAs production in each of the constructed metabolic networks of the CD patients. By adding more of the identified metabolites in the environment in a personalized manner, the authors showed that altering the diet in silico can have different results in each individual based on their microbiota composition [59], which is in line with experimental results [60].

Fig. 5

Representation of the use of empty compartments in studying the gut microbiota. The metabolic compartments indicated with grey squares are organized in a two-dimensional grid. Each bacterium has its own metabolic compartment. Exchange of metabolites takes place between the bacterium and its own metabolic compartment. Metabolic exchange can also take place between adjacent metabolic compartments. The metabolic compartments without any bacterium can still exchange metabolites with the adjacent compartments. In this way, metabolic gradients occur. Introduction of a time step to a grid of compartments gives the opportunity to include movement and division of bacteria. a Bacterial community at the start of a simulation, movement is indicated with black, curved arrows and division is indicated with grey arrows. b After a time step, some cells move and others divide, resulting in a different distribution of bacteria in the grid

Currently, several tools are available for modeling bacterial communities, such as OptCom, BacArena, MiMoSa, COMETS, FLYCOP and MICOM [19, 58, 61,62,63,64,65]. To fully capture the effect of the microbial community on the host, metabolic networks of the gut community need to be extended with GEMs of host cells. A number of GEMs for human tissue such as the liver, blood vessels, and gut epithelial cells are already available to expand the existing metabolic networks of the gut microbiota to include human cells [66]. Combining metabolic networks based on 16S sequencing data with a representative GEM of human metabolism may be used to investigate the effect of the microbiota metabolic output on host health. A representative model of human metabolism is Recon3D, which contains over 13,000 metabolic reactions and can be used to integrate host and microbiota metabolism [24, 25]. However, combining host GEMs and bacterial GEMs can be challenging due to the formulation of the OF and spatial organization. For example, the intestine is spatially organized where intestinal cells surround the microbiota. Accordingly, gut microbiota metabolites have the highest concentration in the lumen and will be first available for the microbiota and not the host cells [67, 68]. When simulating interactions between gut microbiota and the host, metabolic networks should take these gradients into account. Similarly, metabolic and pH gradients exist along the length of intestinal tract [69, 70] and the microbiota composition in the small intestine differs from that in the colon. Furthermore, bacterial cells travel from the upper part of the intestine to the colon to be ultimately shed in the feces [71]. One way to simulate these gradients and spatial organization is through the addition of empty compartments between the compartments of the host cells and the bacterial compartments in a metabolic model. The compartments are organized in a two-dimensional grid and each bacteria has its own compartment from which it can exchange nutrients. The compartments are connected through exchange reactions, thus bacteria can interact with each other. By adding compartments without bacteria in between compartments filled with bacteria, a new bacterium can fill the empty compartment, which can simulate movement or replication of the bacteria (Fig. 5). Moreover, the empty compartments simulate gradients, since not all the nutrients entering the empty compartment will move to the next compartment, which is filled with a bacterium [58]. Another way to simulate gradients and movement of bacteria is by employing differential equations. Van Hoek et al. used a metabolic model in which spatial organization is included to investigate cross feeding. The authors expanded a GEM of Lactobacillus plantarum with metabolic reactions commonly found in the gut such as butyrate and propionate fermentation. Next, the authors placed multiple expanded GEMs of L. plantarum randomly in a tube-like grid environment and used the SOA to simulate the network in 80,000 time steps. For every cell in the grid the OF is maximizing the rate of ATP production, which is an alternative to biomass formation. Periodically, glucose was added to the system at the proximal side of the tube and concentrations of metabolites in each compartment were shifted to the right until they leave the system at the distal side of the tube. Together, this represents metabolic gradients that occur in the gut. The authors then included growth, removal, movement and evolution of single organisms in the network and showed that diversification can be an emergent property of cross feeding among microbial communities. When diarrhea was simulated by increasing the flux speed through the system, the microbial diversity was destroyed [72], which is in line with in vivo observations [73]. Chan et al. used a similar approach, to show that aerobic and anaerobic bacteria separate into different niches based on the oxygen gradient [74]. These examples show the importance of the gut environment and these factors should be taken into account. Persi et al. integrated pH dependent activity of enzymes obtained from experimental data in the BRENDA database in a GEM of cancer cells. The authors showed that cancer cells proliferate differently in silico based on the pH [75]. However, integrating pH dependent activity in a GEM is heavily dependent on experimental work done in vitro. Since most enzymes from gut bacteria are under investigated, integrating pH dependent activity in GEMs used for metabolic modeling is currently not feasible.

Another problem in combining bacterial GEMs and host GEMs is formulating an appropriate OF. In the majority of metabolic networks, the total biomass growth, or metabolic output is optimized. In the gut, the host benefits the most from a balance between the metabolic output of all bacteria, which requires a balance among the distribution of the different microbial species. To address this problem OptCom was developed. This approach uses two layers of OFs. The first layer maximizes the biomass formation of each individual species. The second layer maximizes the growth of the whole community, resulting in a more realistic growth distribution in a bacterial community [61]. The extension d-OptCom can be used in DFBA [76]. However, OptCom cannot easily be used with communities consisting of a large number of bacteria due to increasing computational time needed with increasing community complexity. Therefore, community and systems level interactive optimization (CASINO) was developed. The CASINO framework uses two layers of OFs, but differs from OptCom by optimizing both layers of OFs iteratively [77]. A problem with OptCom and CASINO is the description of the biomass function. As mentioned above, the host benefits the most from a balance in the metabolic output of all bacteria. Thus, for a healthy microbiota the community needs to function at steady-state. However each individual bacterium grows with a specific growth rate [78]. OptCom and CASINO do take bacterial growth into account, but not the community steady state. Therefore, optimizing biomass formation results in the domination of the fastest growing bacterium in the system, resulting in a flux distribution that does not represent the flux distribution in the gut community. To overcome this problem, SteadyCom was developed which takes the community steady-state into account [79].

Simulation of dysbiosis and treatment using metabolic networks

While the composition of the gut microbiota fluctuates over time due to, among others, the diet, the overall composition is more or less stable [80] and is able to recover from short-term perturbations such as short-term antibiotic administration, periods of starvation, and radical changes in diet [28, 81, 82]. However, long term perturbations such as prolonged antibiotic use [82], but also changing diet can cause a shift in microbiota composition, and in turn, a shift in its metabolic products [26, 83]. Alterations in microbiota composition, also known as dysbiosis, can have negative consequences on host health. Dysbiosis has been reported in patients with Alzheimer’s disease [84], Huntington’s disease [85], and Parkinson’s disease (PD) [86]. The microbiome modeling toolbox can be used to construct a metabolic network of the microbiota from the relative abundance data obtained from 16S sequencing studies of fecal samples and can be used to investigate dysbiosis [9]. Baldini et al. used the microbiome modeling toolbox to investigate differences in metabolic output between the microbiota of PD patients and healthy controls. The authors placed GEMs found in the AGORA database in a shared metabolic environment to construct a metabolic network. Differences in bacterial composition between each person were achieved by adjusting the coefficient for the biomass function in the stoichiometric matrix. Maximizing the total output of the microbial community for a produced metabolite was used as the OF. In total, 129 metabolites were investigated for each individual. The average European diet was used as input to simulate changes in the metabolic output. The study reported 9 metabolites to have the potential to be significantly altered in PD patients compared to healthy controls including methionine and cysteinylglycine, which are part of sulfur metabolism. Furthermore the authors showed that a higher presence of Akkermansia muciniphila, and Bilophila wadsworthia in PD patients, and identified a new research target in PD research with the use of metabolic networks [87]. Hertel et al. confirmed these results by using the same approach as Baldini et al. [87] to construct and simulate personalized community networks of 31 early stage, drug naïve PD patients and 28 age matched controls from data obtained from Bedarf et al. [88]. Simulation with FBA concluded that four microbial reactions involved in homoserine metabolism are altered in PD patients, which is consistent with measured levels in plasma of PD patients and homoserine, the precursor of methionine, is part of sulfur metabolism. Furthermore, the authors showed an increase in two bacterial species both involved in sulfur metabolism, B. wadsworthia and A. muciniphila [89]. Both species have been reported in PD progression [90]. The above-mentioned studies show the potential for metabolic networks and FBA in investigating causal relationships in understanding diseases. Additionally, dysbiosis can play a role in the treatment of neurodegenerative diseases. For example, PD patients receive levodopa (L-DOPA) as treatment. However, the dosage varies widely among patients [91]. Members of the gut microbiota can convert L-DOPA. Thus, having more L-DOPA converting bacteria in the microbiota results in a higher dosage of L-DOPA among PD patients [92, 93]. On the contrary, medication used to treat brain-related diseases can have negative consequences by inducing dysbiosis in the gut. For example fluoxetine, a drug used as an antidepressant [94], causes sporulation in members of the gut microbiota. Thus, changing the metabolic output of the microbiota, which may impact host health [95]. The above-mentioned studies show an association between medication, and dysbiosis but lack the causality of dysbiosis and options to treat the dysbiosis to improve drug effectiveness. Constructing metabolic networks of a dysbiosis can give insight in the interactions of the dysbiosis. Next, effectiveness of treatment can be tested in silico on the microbiota of each patient individually, which will improve the process of choosing the best treatment, when prescribing medication.

Currently, dysbiosis can be treated in several ways such as a fecal matter transplantation (FMT), antibiotic use and using pro-, pre-, or psychobiotics. FMT has been used successfully in recurrent Clostridium difficle infection [96], but not as successful in treating IBD [97]. Interestingly, a study investigating FMT in IBD showed that patients treated with fecal material from one particular donor showed more response compared to patients treated with fecal matter from other donors [98]. This suggests that for successful use of FMT the microbial species and metabolites responsible for the beneficial effects of FMT need to be identified. However, to identify the differences between donors and understand why some donors are better than others, not only the microbiota and metabolite composition is needed but also the interaction among the bacteria as well as between the bacteria and the host, indicating the need for the construction of metabolic networks of the donors.

Another possible treatment for dysbiosis is the use of probiotics. Probiotics are defined as live microorganisms that, when administered in adequate amounts, confer a health benefit on the host [99]. Numerous studies have shown positive effects of single strains or mixture of probiotics on host health in animals or humans [100,101,102]. In humans, the effect of probiotic administration is mostly studied in clinical trials, whereby the effect on health is measured [103]. However, these study often show contradictory or unexpected results [104,105,106]. For example, Suez et al., reported that administering probiotics after antibiotic treatment slowed the recovery of the microbiota [81]. Whereas other studies show that administering probiotics after antibiotic treatment does not impact the recovery of the microbiota [107, 108]. These contradicting results can be explained by differences in the probiotic strains used, types of antibiotics, dosage of the treatment, but also diet, medical history, initial microbiota composition and genetics of the patient [104, 109]. To address these problems there is a need to decipher the molecular mechanisms underlying the observed effects of probiotics. Furthermore, the viability of the probiotic in vivo, and the interaction within the microbiota and the metabolic output need to be investigated. In this respect, metabolic networks provide a useful tool to investigate the effect of probiotics on the gut microbiota. When a probiotic strain is added as a compartment in a metabolic network of a gut community, the metabolic flux distribution of the network will change. In other words, certain levels of bacterial metabolites produced in the network will be increased or reduced. In turn, this might have an effect on the microbiota composition presented in the network as shown in Fig. 6. From these changes observed in silico, the effectiveness of probiotic treatment can be predicted [110]. However, to accurately predict the effect of adding a new species to a community, the metabolic behavior of the new species should be known at the strain level and the metabolic capability under different environments should be known. Bifidobacterium species are widely used as probiotics [111]. As mentioned earlier, the GEMs of Bifidobacteria in the AGORA database are not condition-specific and do not give the metabolic behavior in detail. For example the GEMs of Bifidobacteria in the AGORA database do not show growth on starch, whereas most Bifidobacteria can metabolize starch [112]. Therefore, careful curation of the GEMs in the AGORA database is needed before using a FBA approach in probiotic research. Devika et al. refined GEMs of 36 strains of Bifidobacteria including probiotic strains used in commercialized products. The authors compared the metabolic capabilities of the GEMs under 30 different environmental conditions, with maximizing the biomass function as the OF. Based on metabolic capabilities in the different environments the Bifidobacteria could be divided into three groups. Based on the metabolic capabilities of the GEMs the authors hypothesized that the protective effect of the probiotic candidate Bifidobacterium thermophilum RBL67 on Salmonella and Listeria species comes from the production of SCFAs. Furthermore the authors hypothesize that Bifidobacterium gallicum DSM20093 and Bifidobacterium kashiwanohense DSM21854 can help relieve constipation via production of acetate [50], showing that FBA can be a useful tool for identifying new probiotic species.

Fig. 6

Representation of the application of FBA to study the effect of adding probiotic species in a bacterial community consisting of two gut bacteria, F. prausnitzii and B. adolescentis. Fluxes are represented by arrows. Solid black arrows indicate uptake and secretion reactions of the bacteria, dashed black arrows indicate the flow of metabolites in and out of the system and dashed grey arrows indicate the formation of new biomass, where metabolites can no longer be secreted by the bacteria, thus leaving the system. The amount of flux is represented by the thickness of the arrows, a higher flux is a thicker arrow. a The community depicted in Fig. 3Adapted after El-Semman et al. [20]. b Addition of another species to the shared metabolic environment will change the flux distribution of the whole system. The added species might produce a metabolite, which the other species can use, resulting in new products. Similarly, the abundance distribution of the community may change in response to the addition of a new bacterium. A change in abundance distribution of the community is depicted an altered flux through each biomass function

Future perspectives

Metabolic networks and FBA are promising tools to elucidate the precise interactions among bacteria, between bacteria and the host, and can be used to test the effectiveness of probiotic treatment before clinical trials Furthermore, metabolic networks can help point at knowledge gaps. For example, if there is a discrepancy between the outcome of model simulations, and experimental work, a more thorough understanding of the metabolic capabilities of the organism of interest is needed, which leads to new hypotheses and more experiments. Next, metabolic networks can elucidate causality in the microbiota research, because unlike in experimental work, numerous alterations can be implemented in metabolic networks [113]. Currently, metabolic networks of the gut microbiota are mainly based on databases from samples collected from western individuals [114]. To make sure that metabolic networks can be used globally, more global sampling is needed. Furthermore, to ensure the accuracy of the predicted output of the metabolic networks, sufficient information of concentrations of microbial metabolites in the gut, microbial composition, metabolic potential, and interactions is required to be implemented in the metabolic network [50, 115]. Furthermore, metabolic networks are mathematical descriptions of reality. To make sure the metabolic networks depict reality accurately, the outcome of simulating metabolic networks should be experimentally validated. However, the majority of the gut bacteria are not yet culturable [116]. Therefore, GEMs of those bacteria cannot be tested for accuracy. Furthermore, investigating interactions between bacteria can be challenging. Nevertheless, research is performed to identify, culture and characterize new species from the gut microbiota [117]. Currently, new in vitro methods to experimentally validate metabolic networks become available such as SHIME, HuMiX and others [113, 118, 119]. Medlock et al. combined metabolic networks with co-culturing to infer cross feeding between members of the Altered Schaedler Flora (ASF) and validated that this metabolic interaction leads to a growth benefit for the bacteria [120]. This review focused on the use of FBA in modeling the gut microbiota. However, more modeling approaches are being used to understand and improve gut health. An example of the use of a computational modeling approach to benefit human health is a study performed by Zeevi et al. where a machine learning algorithm advised a personalized diet based on, blood glucose levels and microbiota composition, to lower the glycemic response after a meal. Indeed, when participants followed the advised diet, a lower glycemic response was observed [121], indicating that algorithms can predict the influence of diet on the host if data on microbial composition and host parameters are gathered. However, since this approach was based on machine learning, not on FBA, the exact metabolic interactions were not elucidated in the study. Computational fluid dynamics is also used to study the digestion process [122]. Similarly, mathematical descriptions based on experimental data obtained from fermentation experiments are applied to predict the effect of a dietary change [123], and regression based network approaches were shown to be able to identify keystone species in the gut microbiota [124]. Combining knowledge obtained from different computational and experimental methods will help to better understand the intricate relationship between the microbiota and the host. Over time, this understanding can be used to design targeted, and personalized approaches to alter the microbiota in a way that it benefits the health of the host.

Conclusions

Interactions within the complex microbiota community and between the microbiota and the host are challenging to elucidate but understanding these interactions is vital for successful intervention in the microbiota. This review shows that the use of computational methods based on flux balance analysis provides novel understanding of the microbial interactions, helps formulating new testable hypotheses, and enables the investigation of the effectiveness of probiotic administration in a personalized manner.