Metabolic coupling of Mycobacterium tuberculosis to its host is foundational to its pathogenesis. Computational genome‐scale metabolic models have shown utility in integrating ‐omic as well as physiologic data for systemic, mechanistic analysis of metabolism. To date, integrative analysis of host–pathogen interactions using in silico mass‐balanced, genome‐scale models has not been performed. We, therefore, constructed a cell‐specific alveolar macrophage model, iAB‐AMØ‐1410, from the global human metabolic reconstruction, Recon 1. The model successfully predicted experimentally verified ATP and nitric oxide production rates in macrophages. This model was then integrated with an M. tuberculosis H37Rv model, iNJ661, to build an integrated host–pathogen genome‐scale reconstruction, iAB‐AMØ‐1410‐Mt‐661. The integrated host–pathogen network enables simulation of the metabolic changes during infection. The resulting reaction activity and gene essentiality targets of the integrated model represent an altered infectious state. High‐throughput data from infected macrophages were mapped onto the host–pathogen network and were able to describe three distinct pathological states. Integrated host–pathogen reconstructions thus form a foundation upon which understanding the biology and pathophysiology of infections can be developed.
Mycobacterium tuberculosis (M. tb) is an insidious and highly persistent pathogen that affects one‐third of the world's population (WHO, 2009). Metabolism is foundational to M. tb's infection ability and the ensuing host–pathogen interactions. In addition, M. tb has a heterogeneous clinical presentation and can infect virtually every tissue. Depending on the location of the infection, different metabolic pathways are active and inactive in both the host and pathogen cells. In this study, we sought to model the host–pathogen interactions of the human alveolar macrophage and M. tb as well as detail the metabolic differences in specific infection types using genome‐scale metabolic reconstructions (Figure 4A).
Genome‐scale metabolic reconstructions are knowledge bases of all known metabolic reactions of a given organism. Reconstructions have been shown to elucidate the mechanistic genotype‐to‐phenotype relationship through the integration of high‐throughput and physiological data (Oberhardt et al, 2009). Genome‐scale reconstructions are converted into mathematical models under the constraints‐based reconstruction and analysis (COBRA) platform (Becker et al, 2007). COBRA models use network stoichiometry and steady‐state mass balances to define a solution space of potential flux states that a network can take. Thus, the COBRA approach does not require kinetic parameters.
Recently, the global human metabolic network, Recon 1, has been reconstructed (Duarte et al, 2007). To understand the metabolic host–pathogen integrations of M. tb with its human host, we first tailored the global human metabolic network into a cell‐specific metabolic reconstruction of the human alveolar macrophage. This was carried out using established computational algorithms (Becker and Palsson, 2008; Shlomi et al, 2008) and manual curation to confirm the included and excluded reactions. The human alveolar macrophage reconstruction, iAB‐AMØ‐1410, accounts for 1410 genes, 3012 intracellular reactions, and 2572 metabolites (Figure 4C). iAB‐AMØ‐1410 was able to accurately predict maximum ATP and NO production rates obtained from experimental data (Griscavage et al, 1993; Newsholme et al, 1999).
The second step to studying host–pathogen interactions was integration of the human alveolar macrophage reconstruction with an existing genome‐scale metabolic model of M. tb, iNJ661 (Jamshidi and Palsson, 2007). Interfacial constraints were set to create a phagosomal environment that was hypoxic, nitrosative, rich in fatty acids, and poor in carbohydrates. From the onset, it was apparent that some oxygen (<15% of in vitro uptake) was required for proper simulations. In addition, algorithmic tailoring of the M. tb biomass objective function was performed to better represent an infectious state. The integrated host–pathogen metabolic reconstruction was dubbed iAB‐AMØ‐1410‐Mt‐661.
Analysis of the integrated host–pathogen metabolic reconstruction resulted in three main findings. First, by setting interfacial constraints and tailoring the biomass objective function, the solution space better represents an infectious state. Without adding artificial constraints to the host portion of the integrated model, the iAB‐AMØ‐1410 solution space is greatly reduced (Figure 4B). Macrophage glycolysis and nitric oxide production are up‐regulated and macrophage ATP production, nucleotide synthesis, and amino‐acid metabolism are suppressed. In addition, M. tb glycolysis is suppressed and isocitrate lyase is up‐regulated for generation of acetyl‐CoA. Fatty acid oxidation pathways and production of mycolic acids are increased, while production of nucleotides, peptidoglycans, and phenolic glycolipids are reduced. The modified solution space of the alveolar macrophage and M. tb better represents the infectious state.
Second, the host‐pathogen model more accurately predicts M. tb gene deletion tests than the current in vitro model, iNJ661. The host‐pathogen model predicted 11 essential genes and 37 unessential genes differently than iNJ661. A total of 22 of the differentially predicted genes have been experimentally characterized (Sassetti and Rubin, 2003; Sohaskey, 2008). The host‐pathogen model correctly predicted 18 of the 22 genes. Thus, iAB‐AMØ‐1410‐Mt‐661 is a more accurate platform for studying infectious states of M. tb.
Finally, we sought to determine metabolic differences in both the macrophage and M. tb between three different types of infection: latent, pulmonary, and meningeal. Transcription profiling data of the macrophage for the three infections (Thuong et al, 2008) were integrated in the context of the host–pathogen network to elucidate the reaction activity of the three infections. There was wide heterogeneity in the three infection states; some of these differences are highlighted. Macrophage hyaluronan synthase and export were only active in the pulmonary infection. This is potentially interesting from a pharmaceutical viewpoint as hyaluronan has been implicated as a potential carbon source for extracellular M. tb (Hirayama et al, 2009). In addition, we detected metabolic activity differences in M. tb pathways that have been previously discussed as potential drug targets (Eoh et al, 2007; Boshoff et al, 2008). Polyprenyl metabolic reactions were only active in the latent state infection, while de novo synthesis of nicotinamide cofactors was only active in latent and meningeal M. tb infections.
Host‐pathogen modeling represents a novel approach for studying metabolic interactions during infection. iAB‐AMØ‐1410‐Mt‐661 is a more accurate platform for understanding the biology and pathophysiology of M. tb infection. Most importantly, genome‐scale metabolic reconstructions can act as scaffolds for integrating high‐throughput data. Particularly, in this study we were able to discern reaction activity differences between different infection types.
A human alveolar macrophage genome‐scale metabolic reconstruction was reconstructed from tailoring a global human metabolic network, Recon 1, by using computational algorithms and manual curation.
A genome‐scale host–pathogen network of the human alveolar macrophage and Mycobacterium tuberculosis is presented. This involved integrating two genome‐scale network reconstructions.
The reaction activity and gene essentiality predictions of the host–pathogen model represent a more accurate depiction of infection.
Integration of high‐throughput data into a host‐pathogen model followed by systems analysis was performed in order to elucidate major metabolic differences under different types of M. tuberculosis infection.
Mycobacterium tuberculosis (M. tb) is a highly persistent pathogen that primarily affects the third world. About one‐third of the world's population is infected, with 9.27 million new cases and 1.76 million deaths in 2007 (WHO, 2009). Furthermore, the development of multidrug‐resistant tuberculosis and extensively drug‐resistant tuberculosis, which are infections that cannot be treated with first‐line and second‐line drugs, respectively, continue to keep M. tb a concern in developed countries.
A key aspect of M. tb's pathogenicity is its ability to dramatically change its metabolism in different states. After infecting an alveolar macrophage, M. tb shifts into an infectious state, stopping biomass accumulation. M. tb also accumulates mycolic and fatty acids on its cell wall in order to survive a very hostile phagosome environment that is nutrient poor, hypoxic, nitrosative, and oxidative (Schnappinger et al, 2003). The mycolic and fatty acids are essential to M. tb's resistance to drug therapies (Barkan et al, 2009). Under a latent state, the pathogen cannot be spread and the patient is not in danger. However, latent tuberculosis can activate in the lungs (75% of cases) or other parts of the body.
The lack of proper experimental M. tb models that provide a mechanistic understanding of in vivo conditions hinders a better understanding of the infection process (Boshoff and Barry, 2005). It is not only difficult to work with M. tb because of its slow growth rate, but most in vitro models are inaccurate simulations of in vivo conditions. In vivo animal models better characterize the disease, but there is less control over experimental conditions.
Genome‐scale metabolic reconstructions are useful for increasing the understanding of the genotype/phenotype relationship in organisms (Oberhardt et al, 2009). These networks are built in a bottom–up manner in which the nodes (metabolites) are connected by links (biochemical transformations and reactions) as defined by genetic and physiological data (Palsson, 2004; Reed et al, 2006). There are now genome‐scale reconstructions for numerous prokaryotes (Thiele et al, 2005; Feist et al, 2007), including two such networks for M. tb (Beste et al, 2007; Jamshidi and Palsson, 2007) and eukaryotes (Duarte et al, 2004; Sheikh et al, 2005). In addition, the global human metabolic network, Recon 1, has been reconstructed (Duarte et al, 2007) containing all the annotated biochemical reactions for human cells. Recon 1 has sparked interest in building specific networks for different human tissues and cells.
As the number of genome‐scale reconstructions increases, there has been interest in building networks that characterize the metabolic interaction between multiple organisms. In this study, we looked to increase the understanding of the metabolic changes in both the host (alveolar macrophage) and the pathogen (M. tb) during infection through the use of a host–pathogen genome‐scale reconstruction. We did this by first building a manually curated cell‐specific human network for the alveolar macrophage. We then integrated the host model (HM) with the pathogen model (PM) (iNJ661) to form the host–pathogen model (HPM). We used established constraint‐based analysis methods (Becker et al, 2007) and published gene expression data for M. tb infections to further our understanding of macrophage and M. tb metabolic functions.
Interrogation of the interactions between an HM and PM required completing four steps: (1) construction of the alveolar macrophage HM, (2) adaptation of the M. tb PM, (3) integration of the two into an HPM, and (4) analysis of the resulting HPM in different contexts. These analyses led to three general observations,
Integration of the PM into the HM both reduced the size and altered the solution space of both models. The changes in the PM were primarily because of the interfacial constraints set on the network.
Analysis of the PM portion of the HPM improved gene essentiality predictions.
The HPM enabled analysis of complex data sets from complex experimental conditions and was able to highlight differences in metabolism through the comparison of the three infectious states.
Constructing the HM: human alveolar macrophage, iAB‐AMØ‐1410
Reconstruction of the global alveolar macrophage network was a time intensive two‐step process beginning with algorithmic tailoring of Recon 1 followed by a set of manual curation steps (Figure 1). Using gene expression data from inactive macrophages, maximum enzyme fluxes, and exchange data for the macrophage (see Supplementary information), two preliminary models were built using two different previously published algorithms (Becker and Palsson, 2008; Shlomi et al, 2008). These algorithms produced context‐specific, not global, models that were the best fit of reaction flux and pathway length to the gene expression and exchange data. All other inactive and dead‐end reactions were removed by these algorithms. Unfortunately, macrophages are known to have varied gene expression states, such as during infection and inflammation; therefore, the algorithms do not take into account other potential states. In fact, both algorithms failed to include nitric oxide synthase, which is critical for nitric oxide production. This is due to the initial gene expression data used for inactive alveolar macrophages. Hence, a global, consensus model was constructed through manual curation and reconciliation of the two models. Manual curation was conducted using primary literature, enzyme, and immunohistological databases, and network structure and requirements. The global model included all inactive and dead‐end reactions to provide opportunities for further research with different exchange constraints and to allow for future gap filling, respectively. The resulting HM was named iAB‐AMØ‐1410 for i (in silico), AB (the primary author's initials), AMØ (alveolar macrophage), 1410 (number of open reading frames).
Characterizing the HM, iAB‐AMØ‐1410
The iAB‐AMØ‐1410 model has many of the capabilities of Recon 1 because of its high number of reactions (intracellular: 3012) and active genes (1410) (Figure 2A). In addition, the HM maintains a high reaction count from all major subsystems (Figure 2B). Macrophages are known to be highly metabolically active as well as have many different metabolic states because of varied patterns of gene expression (Gordon, 2007). A biomass maintenance function was formulated (see Materials and methods) and used to constrain the HM. The metabolic components of the biomass maintenance function are detailed in the Supplementary information. As macrophages do not readily multiply, the biomass maintenance function used here represents cellular maintenance requirements, such as lipid, protein, mRNA turnover, DNA repair, and ATP maintenance.
A series of tests were used to validate the macrophage model. We first tested ATP production of the macrophage model. We optimized the HM for ATP production yielding a flux of 0.6843 mmol/h/g cell DW (4C). In vitro experiments with macrophages yielded an ATP generation rate of 0.7121 mmol/h/g cell DW (Newsholme et al, 1999). This puts the HM at a production accuracy of 96.1%. Second, we simulated nitric oxide production, which is another important metabolic capacity of macrophages. Maximum nitric oxide production in silico was 0.0359 mmol/h/g cell DW, which is 1.7% lower than the maximum in vitro rate of 0.0365 mmol/h/g cell DW (Griscavage et al, 1993). The macrophage is also highly anaerobic in its growth despite its uptake of oxygen (Newsholme et al, 1999). Macrophages have high production rates of lactate (Burke and Lewis, 2002) because of high glucose oxidation very similar to the Warburg effect (Warburg, 1956). Consistent with this phenotype, the HM predicts a high glucose oxidation rate and produces 0.3093 mmol/h/g cell DW of lactate.
Recon 1 was characterized and validated by showing its potential in accomplishing 288 metabolic functions with the proper exchange constraints. A similar assessment of the HM was performed using the same 288 metabolic functions of Recon 1 (see Supplementary information). Unlike the reaction count (Figure 2B), the HM does not preserve the comprehensive metabolic functions of Recon 1 (Figure 2C), as would be expected. Important metabolic functions dealing with carbohydrate and central metabolism are well preserved. However, several peripheral metabolic functions such as amino‐acid, lipid, cofactor, and nucleotide metabolism were not preserved. There are two reasons for this. First, removing a few reactions from Recon 1 to build the HM does not significantly affect the reaction count, but it can have significant effects on functional pathways. Second, the 288 metabolic functions were originally characterized in Recon 1 using open exchange constraints. In this study, we ran all simulations using the macrophage exchange constraints, which can have substantive effects on output. For example, most metabolic functions in the miscellaneous category dealt with growth under different substrates. As there was no literature support for the macrophage to import some of these metabolites, the simulations were deemed as failures.
Infection in silico, the HPM: iAB‐AMØ‐1410‐Mt‐661
The alveolar macrophage network, iAB‐AMØ‐1410, was integrated with the M. tb network, iNJ661. The integration involved three steps: mathematical integration, setting interfacial constraints, and revision of the objective function (Figure 3). Integration had a few preparatory steps. We renamed reactions and metabolites for proper compartmentalization, added transporter reactions that linked the alveolar macrophage cytosol with a new phagosome compartment, and modified iNJ661 to uptake metabolites from the phagosome rather than the extracellular compartment. M. tb blocks certain signaling pathways to hinder maturation of the phagosome into the phagolysosome (Pieters, 2008). Without the formation of the phagolysosome, M. tb is able to survive within the phagosome. After these preparatory steps, the two stoichiometric matrices were mathematically combined. In doing so, we created a modularized network with M. tb residing in the macrophage phagosome, using resources from its host (Figure 4A). The final model was named iAB‐AMØ‐1410‐Mt‐661, which added the first letter of the genus and species name (Mt) of the pathogen and 661 for the number of open reading frames.
Upon integration, it was readily apparent that a proper infectious state could not be simulated unless interfacial constraints were set between the HM and PM, as the PM would take up any metabolite from the HM. Thus, we set interfacial exchange constraints based on literature. In the phagosome, M. tb is under hypoxic conditions (Honer zu Bentrup and Russell, 2001) and oxidative stress from metabolites such as nitric oxide (Schnappinger et al, 2003). In addition, glucose is depleted and the pathogen is dependent on glycerol and fatty acids (McKinney et al, 2000; Schnappinger et al, 2003) as a carbon source. Interestingly, the PM would not grow in silico unless there was a minimal oxygen uptake, consistent with functional hypoxia as opposed to anoxia. M. tb biomass production was compared with M. tb oxygen uptake showing a linear dependence on oxygen until an uptake of 0.129 mmol/h/g cell DW was reached, in which additional oxygen did not increase the biomass value (Supplementary Figure S1). The uptake represented 13.2% of the in vitro uptake rate (Jamshidi and Palsson, 2007), which was used for simulations. The oxygen M. tb biomass curve is consistent with existing knowledge that reduced oxygen is associated with slowed M. tb growth and disease progression (Rustad et al, 2009). In addition, in vitro studies lowering partial pressure of oxygen below 1% of normal conditions have shown to fully halt M. tb replication (Yuan et al, 1998). A detailed listing of the exchange constraints of in vitro (iNJ661) and in vivo (iAB‐AMØ‐1410‐Mt‐661) networks for M. tb is provided in the Supplementary information.
The third integration step involved revision of the iNJ661 objective function (biomass) to better represent an acute infection. The original biomass function for the PM was formulated to represent rapid growth in mid‐log phase. In mid‐log phase, the pathogen tries to accumulate biomass by dividing rapidly. However, under in vivo conditions, the pathogen survives in the phagosome's harsh environment, resulting in an altered metabolic state. In order to build a new objective function, we used M. tb gene expression data derived from in vivo mouse model studies (Talaat et al, 2004; Shi et al, 2005) as well as infection simulated in vitro studies (Bacon et al, 2004; Hampshire et al, 2004; Voskuil et al, 2004). Though the in vivo expression profiling studies came from mouse and in vitro simulated models, the studies represent the best available transcriptional profiling of M. tb during an infectious state. There have been some noted differences, especially in nitric oxide production mechanisms between murine and human macrophages (Schneemann and Schoeden, 2007), but the human alveolar macrophage is able to produce nitric oxide in a similar manner to murine macrophages (Rich et al, 1997).
The biomass reaction is a demand function that pulls resources from the metabolic model. Each metabolite in the biomass affects reaction activity and hence, the final solution space. We hypothesized that by using the known differential reaction activity determined by differentially expressed genes, an altered infection‐specific solution space could be defined by tailoring the PM's biomass objective function. Using randomized sampling, we characterized the solution space for each metabolic component in the biomass function, individually. Using linear regression in an iterative manner, we added and removed metabolic components of the biomass objective function to develop a new objective function (Table I) with a solution space that best fit the gene expression data. Particularly, we removed phenolic glycolipids from the M. tb infection objective function. These computationally driven revisions were confirmed with established experimental data that showed that M. tb H37Rv clinical strains lack biosynthesis of phenolic glycolipids (Constant et al, 2002).
Characterizing the metabolism of the HPM
The macrophage component of the integrated HPM is similar to iAB‐AMØ‐1410, with an M. tb component, which is starkly different from iNJ661. The integration radically affected the major objective function fluxes of the macrophage and M. tb components (Figure 4C). The maximum macrophage biomass was slightly reduced (∼12%), while maximum ATP, nitric oxide, and NADH production were significantly reduced (∼75, ∼55, and ∼70%, respectively). Such a large reduction in ATP, nitric oxide, and NADH production points toward a weakened or compromised alveolar macrophage. The M. tb biomass function was also greatly reduced. We then compared the size of the macrophage solution space in the HM and the HPM (Figure 4B). There was a significant reduction in the size of the solution space (51% reduction of the mean reaction flux span). By adding the M. tb model to the alveolar macrophage reconstruction, we were able to substantially decrease the solution space without adding any additional constraints on the internal reactions of the macrophage portion of the network. Only two reactions had a larger flux span: the exchanges of sulfate and phosphate. The reactions with the largest flux span reduction (>70%) dealt with fatty acid, lipid, and proteoglycan metabolism as well as transport reactions. The M. tb network tightened the flexibility of these reactions because of its need for fatty acids and glycerol.
Changes in reaction activity of iAB‐AMØ‐1410‐Mt‐661
Using randomized sampling, we identified the distribution of all feasible steady states for each reaction in iAB‐AMØ‐1410, iNJ661, and iAB‐AMØ‐1410‐Mt‐661. Statistically significant up‐ and down‐regulated reactions in the integrated HPM were identified after in silico phagocytosis of M. tb. The majority of differences are seen in the M. tb portion of the integrated model, which was expected because of the changes in exchange constraints and revised objective function. However, there are also some differences in the alveolar macrophage portion of the network.
The changes in flux states of the M. tb portion of the HPM show a shift in carbon uptake and overall usage. There are major shifts in flux states in central metabolism. Glycolysis is suppressed (Figure 5, ENO) with the production of acetyl‐CoA coming from fatty acids through the glyoxylate shunt (Figure 5, ICL). In addition, glucose is generated from gluconeogenesis (Figure 5, FBP). Beyond central metabolism, there are several changes. For example, there is an up‐regulation of fatty acid oxidation pathways, as fatty acids become a major carbon source. In addition, there is a shift toward mycobactin and mycolic acid synthesis. Production of nucleotides, peptidoglycans, and phenolic glycolipids is also reduced. Thus, the model more accurately described an infectious state for M. tb. Overall, 75 reactions were up‐regulated, while 128 reactions were suppressed. A full listing of reactions, associated genes, and their flux states are provided in the Supplementary information.
There are fewer changes in the alveolar macrophage portion of the integrated model. Only 8 reactions were up‐regulated and 46 reactions were down‐regulated. There is a higher reliance on glycolysis as reported by a much greater activity through phosphoglycerate mutase and enolase. In addition, reactions in the pathways for nitric oxide production are up‐regulated. However, the majority of the changes were down‐regulations of ATP production, nucleotide synthesis, and amino‐acid metabolism. The changes in flux states of the alveolar macrophage during infection point toward a suppressed metabolic state that has a high reliance on glucose oxidation.
Gene essentiality predictions
Interestingly, by modeling the host and pathogen together, we were able to increase the accuracy of gene deletion tests. Single gene essentiality predictions were completed for the M. tb portion of the HPM and compared with the results from iNJ661. Model predictions were improved for essential as well as non‐essential genes when compared with experimental data. In total, the PM predicts 188 essential metabolic genes, while the HPM predicts 162. Predictions from the two networks overlap with 153 similarly predicted genes (see Supplementary information). The non‐overlapping predictions are shown in Tables II and III.
The changes in gene essentiality were compared with published in vivo experimental data (Sassetti and Rubin, 2003). Using Transposon Site Hybridization (TraSH), Sassetti and Rubin identified M. tb genes that are essential for survival in a mouse lung model. A total of 374 of 2972 M. tb genes assessed using TraSH overlap with the metabolic models used here. For the nine genes that were predicted as essential solely in the HPM, two are supported by the TraSH data, while five contradict it. Four of the contradictions are components of nitrate reductase (narGHIJ). Though the TraSH data does not support it, Sohaskey (2008) has shown that nitrate reduction enhances M. tb survival under sudden stresses of hypoxia and nitric oxide, which M. tb is under during the infectious state. In addition, our gene essentiality results are further validated by in vitro M. tb tests (Pinto et al, 2007) of the essentiality of sulfite reduction (sirA).
Thirty‐five genes are predicted to be essential in the PM, but non‐essential in the HPM. Of the 35 genes, 12 are also non‐essential in TraSH. Only two genes (murI, embA) conflict with the TraSH data. The majority of the genes deal with peptidoglycan metabolism and membrane metabolism. Such a change from in vitro conditions shows a shift from division and production of biomass to a more specific infectious state dealing with a harsh phagosome environment. Hence, genes that deal with membrane and cell wall components are less likely to affect the survival of M. tb in vivo.
The changes in our new objective function (Table I) deduced by gene expression data is key for the non‐essential predictions. The 35 genes are predicted as non‐essential only in the HPM with the revised objective function for M. tb. Removal of peptidoglycan and phenolic glycolipid components of the objective function increased prediction accuracy with the HPM as evidenced with better correspondence with in vivo data (TraSH), which also further supported the condition‐specific objective function reformulation.
Host–pathogen interactions in different infectious states
M. tb is an insidious pathogen and can remain latent for many years before actively infecting a host. Furthermore, its clinical presentation can be extremely varying and it can infect every organ system. Antibiotic treatment regimens are often tailored for different types of tissue infections, based on their ability to penetrate different tissues. In an analogous manner, infections in different tissues may involve activation of different metabolic pathways in the host and pathogen. Hence, depending on the nature of the infection and the type of tissue affected, different treatment regiments may be appropriate. We sought to determine the activity of metabolic pathways in both the macrophage and M. tb in order to better understand M. tb infections and predict potential drug targets. In this study, we analyzed differences between active infections in two different tissues as well as latent tuberculosis.
An important application of genome‐scale metabolic network reconstructions is to provide a biologically meaningful context for the analysis of large data sets, that is a ‘context for content’ (Oberhardt et al, 2009). High‐throughput data such as gene expression profiles can be mapped to the reconstructions and subsequently analyzed to evaluate functional consequences of changes in expressions levels of genes. We mapped expression data (Thuong et al, 2008) of macrophages from different M. tb infectious states (latent, pulmonary, meningeal) onto the HPM (iAB‐AMØ‐1410‐Mt‐661) and built three infection‐specific models: iAB‐AMØ‐1410‐Mt‐661L (HPM‐L), iAB‐AMØ‐1410‐Mt‐661P (HPM‐P), and iAB‐AMØ‐1410‐Mt‐661M (HPM‐M), which can then be compared based on the differentially active and inactive reactions in the three states in both the macrophage and M. tb determined by flux variability analysis (Figure 6). A full detailing of the approach is provided in Materials and methods.
The majority of differences in active reactions were found in the macrophage (Figure 6B and C). This is expected because of the expression profile data coming from the macrophage. We saw two main discrepancies in active macrophage reactions between the three different types of infections. First, hyaluronan synthesis and exchange was active only in the HPM‐P. This suggests that hyaluronan is produced and secreted into the extracellular space. Hyaluronan is important for cell proliferation, especially for progression of some malignant tumors (Itano et al, 2008). It seems that hyaluronan could have an important function as a carbon source for pulmonary M. tb. In fact, it has been shown that hyaluronan synthase is active and important for extracellular replication of M. tb (Hirayama et al, 2009). As hyaluronan synthase is not critical for macrophage survival, inhibiting hyaluronan synthase (has123) could be a potential method for stopping activation of M. tb from a latent to active pulmonary state.
Second, we saw activation of vitamin D and folate metabolism in the active infectious states (HPM‐P and HPM‐M). Vitamin D injections were once used as a pre‐antibiotic era treatment for tuberculosis (Martineau et al, 2007). It is proposed that vitamin D is crucial for increasing nitric oxide production in the macrophage. In the active infectious states, the macrophage is readily fighting the infection and hence has active vitamin D metabolism. Folate is important as an enzyme cofactor for DNA synthesis and repair. It has been shown to be up‐regulated during activation of macrophages for arthrosclerosis (Antohe et al, 2005). The HPM‐P and HPM‐M predict that this also occurs for active M. tb infection.
Finally, we saw fewer but interesting changes in reaction activity of the M. tb portion of the three context‐specific HPMs (Figure 6B and C). Polyprenyl metabolism was only active in the HPM‐L and de novo synthesis of nicotinamide cofactors was active in latent and meningeal M. tb infections. Polyprenyl metabolism is important for the cell wall and it makes sense that it is active in the latent state when M. tb builds up its cell wall. Both of these pathways have been mentioned as potential drug targets for tuberculosis (Eoh et al, 2007; Boshoff et al, 2008). The lack of activity of the polyprenyl and nicotinamide pathways in the HPM‐P and the polyprenyl pathway in the HPM‐M suggests that such drug targets would not be as effective in eradicating the pathogen in a pulmonary or meningeal infection. However, the activity of these two pathways in the HPM‐L suggests that such drug targets could be potentially viable for latent infections. Using high‐throughput data from only one organism in the host–pathogen network, we were able to infer resulting changes in the other organism through simulations. This is critical for host–pathogen interactions in which high‐throughput data cannot be properly attained for one of the organisms. A listing of all the reactions in each of the three infection‐specific models is provided in the Supplementary information.
There has been growing interest in genome‐scale metabolic reconstructions and constraint‐based analysis over the past 10 years. Development of a global human metabolic reconstruction, Recon 1, as well as reconstructions of multiple human pathogens enabled us to simulate how metabolism in the host and pathogen changes as a result of infection.
In order to carry out this study, a manually curated human alveolar macrophage model of cellular metabolism, iAB‐AMØ‐1410, was created using Recon 1. Purely algorithmic approaches were inadequate to achieve a physiologically relevant model, hence extensive manual curation was used. Flux predictions for ATP and NO production for the macrophage were within 5–10% of experimental measurements. In addition, the in silico model exhibited high glucose oxidation and lactate production, also seen experimentally.
The HM was then integrated with an M. tb reconstruction, resulting in an HPM. We compared the two reconstructions to different experimental data sets for macrophages and M. tb to validate our models. Interestingly, integration of the pathogen with the host immediately shrunk the solution space without the imposition of additional internal constraints. As the objective function of M. tb in the infectious state is necessarily different than the biomass growth function, gene expression data was used to help define the metabolic objectives in the infectious state. These predictions were validated through comparison with in vivo gene essentiality studies (TraSH).
Among the observations of the HPM, it was noted that there was a higher reliance on glycerol and fatty acids in a hypoxic environment that relied on gluconeogenesis to produce sugars. The integrated model was reliant on the glyoxylate shunt. In addition, nucleotide synthesis was suppressed with increased mycolic acid synthesis. As with in vitro simulations, blocking nitrogen and sulfur reduction pathways were shown to significantly affect survival rates of the pathogen. Analysis of the HPM model also showed improved gene essentiality predictions compared with the PM alone.
We used gene expression data to specify the context for analysis of three different types of M. tb infection: latent, pulmonary, and meningeal. There were key differences in both the metabolism of the macrophage and M. tb in the different infection types. iAB‐AMØ‐1410‐Mt‐661 predicted hyaluronan synthase in the macrophage to be a potential drug target for inhibiting pulmonary M. tb infection. In addition, we showed that mapping high‐throughput data on one reconstruction of a host–pathogen network can elucidate physiologically relevant differences on the other reconstruction in the network. Polyprenyl metabolism and certain nicotinamide pathways were shown to be only active in latent and meningeal infections. Drug targets for these pathways have been discussed previously and should be only looked into for latent M. tb infections. A recent growing amount of evidence has shown M. tb to actively metabolize cholesterol (Brzostek et al, 2009; Hu et al, 2010). Interestingly, some reactions involved in cholesterol metabolism were highlighted by the model. As these pathways have not yet been fully characterized and were not in the current PM, the significance and interpretation of these findings were not clear at this point, but point to an area worth expanding the scope of the model and directing further future investigation.
Metabolic network reconstructions have many different uses and can be viewed as platforms for data interpretation, analysis, and prediction. Metabolism is being increasingly recognized as important in the role of pathogenicity (Kafsack and Llinas, 2010). Thus, metabolic HPMs are essential for detailing the molecular interaction between human cells and pathogens because of the importance of metabolic changes in both the human host cell and the pathogen. Host–pathogen reconstructions have several applications including providing a context for content analysis, biological discovery, and predictive power. Such reconstructions can serve as scaffolds for mapping high‐throughput data onto the host or pathogen and then determining the changes in the other organism. In addition, discrepancies in experimental data and in silico predictions can help elucidate unknown interfacial interactions between the host and pathogen. For example, we showed that some oxygen uptake was necessary for M. tb infection, consistent with functional hypoxia, as opposed to strictly anaerobic versus aerobic conditional dependence. Finally, host–pathogen reconstructions take a step toward being a more biologically faithful platform for predicting potential drug targets as the network better represents in vivo conditions. For example, the models can be used to serve as platforms for understanding the host–pathogen interactions through dual perturbation experiments in silico (Jamshidi and Palsson, 2009) for screening responses to different drugs simultaneously with different growth conditions, which could then be complemented with in vitro or even ex vivo experiments.
There are numerous opportunities for improvement and discovery in HPMs. The host and pathogen integration was performed in a modularized manner, thus allowing for replacement of the M. tb component with other bacterial pathogens that infect macrophages such as Salmonella typhimurium (AbuOun et al, 2009) or viral pathogens such as the human immunodeficiency virus. In addition, future revisions of the M. tb reconstruction that further expand the scope of the metabolic processes included in the current model, transcription/translation machinery, and/or Toll‐like receptors (Li et al, 2009) can be fitted within this framework. Another interesting application of the integration framework involves building bacterial community networks that show the interaction between more than two organisms in different environments, such as the gut. We hope that this successful demonstration of the integration of two genome‐scale models will spur future computational explorations into multi‐organism models to complement experiments and drive forward mechanistic understanding as well as generate new hypotheses.
Materials and methods
Construction of iAB‐AMØ‐1410
The reconstruction is derived from the global human metabolic network, Recon 1. Recon 1 consists of 1496 ORFs accounting for 3012 intracellular reactions and 2583 metabolites (Figure 2A). Recon 1 was used as described in Duarte et al (2007), with a few revisions. The originally described Complex IV reaction was replaced with three reactions. This was performed to ensure that the reaction was charge balanced, at the cost of decoupling superoxide production from the cytochrome c oxidase reaction (i.e. any coupling would have to be fixed by the modeler adjusting the relative constraints on the equations). In addition, the trans‐mitochondrial phosphate transporter was re‐written in the electroneutral form. The reaction formulas are explicitly shown in the Supplementary information.
An overview of the model building process is shown in Figure 1. Cell‐specific gene expression data for healthy, inactivated alveolar macrophage was obtained from Gene Expression Omnibus. The specific gene expression study (Kazeros et al, 2008) had a healthy patient sample size of n = 11. Presence/absence calls were made for the gene expression data using the PANP software package in the R computing platform. The PANP algorithm used signal intensities to determine the presence of transcripts. Exchange constraints as well as upper bounds for fluxes of major pathways of central metabolism were set from literature (Newsholme et al, 1986, 1999; Sato et al, 1987; Curi et al, 1988). The expression data was then mapped to the reactions and using two separate algorithms, GIMME (Becker and Palsson, 2008) and Shlomi‐NBT‐08 (Shlomi et al, 2008), we built two preliminary cell‐specific models. The two preliminary models were reconciled with primary literature databases (BRENDA (Chang et al, 2009) and HPRD (Keshava Prasad et al, 2009)), immunohistological staining data (Human Protein Atlas (Berglund et al, 2008)), transporter databases (HMTD (Yan and Sadee, 2000)), primary literature (see Supplementary information), and network stoichiometry to build a final, fully curated reconstruction.
The biomass was formulated similar to the methods described in previous eukaryotic reconstructions of Saccharomyces cerevisiae and Mus musculus (Famili et al, 2003; Sheikh et al, 2005). The biomass component percentages were determined from literature (protein/sugar (Iyengar and Vakil, 1985), lipid (Schmien et al, 1974), and DNA (Grutman and Orgel, 1970)) and are presented in the Supplementary information. The total composition of each biomass element with respect to the whole cell, except for RNA, was found in literature for macrophages. The total RNA content was adapted from the yeast and mouse reconstructions. All components were found in macrophages, but the RNA amount was adapted from the yeast and mouse reconstructions. DNA composition was set to 41% GC content. RNA and protein composition was determined by averaging the sequence composition of expressed transcripts from the gene expression data used to build the HM. Sugar was assumed to be stored solely as glycogen. Lipid components were adapted from studies on human alveolar macrophages (Sahu and Lynn, 1977). ATP maintenance requirements were calculated similar to the method for the mouse reconstruction. All calculations for biomass formulation are provided in the Supplementary information. The biomass reaction was not used as the objective function. Instead, a lower bound was set as a maintenance requirement for the HM at physiologically relevant levels (Mbawuike and Herscowitz, 1989) for all computational tests. Objective functions used for the HM include ATP and nitric oxide production.
Integration of iAB‐AMØ‐1410 and iNJ661
The integration workflow is shown in Figure 3. Before integrating iAB‐AMØ‐1410 and iNJ661, the reactions and metabolites were renamed to delineate between organisms and compartments. Once the reactions and metabolites were recompartmentalized, the stoichiometric matrix of the PM was added to the right‐hand side of the HM. The external environment of the PM was renamed as the phagosome environment. Transporter reactions were added that connected the cytoplasm of the HM with the phagosome environment for the exchange metabolites. The PM's exchange reactions were set to zero, except for ferric ion, copper ion, and hydrogen. These metabolites are not present in the HM and the exchanges were used as sinks because of the PM's need. The transporter reactions provided metabolites known to be in the phagosome to the PM, including ions (Wagner et al, 2005), nitric oxide, fatty acids, traces of glycerol, and oxygen (10% of iNJ661). A full list of metabolites available in the phagosome is provided in the Supplementary information.
The reconstructions were represented mathematically using the stoichiometric matrix (S). The matrix is m by n, containing m metabolites and n reactions. The relationship of reactions and the concentration time derivatives is as shown:
where x and v denote the concentrations of metabolites and fluxes of reactions, respectively. Steady‐state flux‐balance analysis is accomplished using the following linear optimization problem:
The optimization problem is constrained by thermodynamic and enzyme kinetics using upper and lower bounds (ub, lb) on the fluxes of reactions. The optimization vector (c) is a zero‐vector except for a one assigned to the reaction or function to be optimized for. A full description of flux‐balance analysis and its biological uses can be found in earlier publications (Varma and Palsson, 1994).
The flux span is calculated by the difference of the maximization and minimization of the linear optimization problem above for all reactions in the models. Reactions that carried no flux and were perceived to be involved in thermodynamically infeasible loops were not considered when calculating the flux span.
Characterizing the macrophage's metabolic functions
The global human metabolic network contained 288 metabolic functions. Sink reactions were added for the metabolite being produced in each metabolic function. Flux‐balance analysis was completed using the SimPheny software platform and the pathways were detected by manual inspection of flux maps. The original exchange constraints for iAB‐AMØ‐1410 were used for all simulations. A listing of all the simulations and results are provided in the Supplementary information.
Reaction activity tests
iAB‐AMØ‐1410, iNJ661, and iAB‐AMØ‐1410‐Mt‐661 were randomly sampled using a Monte Carlo method detailed below. The macrophage objective functions were maintained at the aforementioned physiologically relevant levels. The tuberculosis objectives were maintained at 75% of the maximum values.
As substrate uptake rates for in vitro versus in vivo M. tb are significantly different, the sampling results required normalization to properly predict differential expression of metabolic pathways. We assumed that the cells will use a slightly different set of reactions for the two growth conditions. In addition, the total flux through all pathways will be constrained by the amount of enzyme, assuming there is adequate substrate. To represent this, we normalized the sample points to have the same median total network flux level by summing the absolute flux through all tuberculosis reactions for each sample point. All points are scaled by the median network flux.
Differential reaction activity was determined under the assumption that if the distributions of candidate flux states, from sampling, for the same reaction under two different conditions do not significantly overlap, then the cell will have to adjust the gene expression to compensate for the complete shift in the solution space. For each tuberculosis reaction, a P‐value was computed, comparing the distributions of normalized sample points of the two models. Significance of change was determined at a false discovery rate of 0.05 (Storey and Tibshirani, 2003). The genes of the significantly changed reactions were determined from the gene–protein reaction association data and compared with published gene expression data.
Monte Carlo sampling
Monte Carlo sampling was used to generate a set of uniform, feasible flux distributions (points). The method is a modified version of the artificially centered hit and run (ACHR) algorithm (Kaufman and Smith, 1998). Initially, a set of non‐uniform pseudo‐random points, called warmup points, is generated. In a series of iterations, each point is moved in a random manner to a new point, always remaining within the solution space. This is performed by (1) choosing a random direction, (2) computing the limits of how far one can travel in that direction and (3) choosing a new point randomly along this line. After many iterations, the set of points will be mixed and approach a uniform random sample of the solution space.
Warmup points are generated by linear programming. For each point, the objective coefficients are set to a random vector with values in [−1,1]. This generates a point at random corners in the solution space. The direction of movement is chosen as described in Kaufman and Smith (1998). The center point of all points is computed and the direction is the difference of a randomly selected point and the center point. This has the effect of biasing the directions in the longer directions of the solution space and speeds up the rate of mixing while maintaining sample uniformity.
One of the problems with the ACHR is that the termination condition is not clearly defined. Here, we introduce the concept of the mixed fraction as a measure of how many iterations are required until proper mixing. A partition is created over the set of points by drawing a line at the median value with half the points on either side of the line. The mixed fraction is a count of how many points cross this line between the beginning of sampling and the end as a fraction of the total number of points. Initially, the mixed fraction is 1 as all points are on the same side of the partition. When perfect mixing is achieved, each point has a 50% chance of crossing the partition line so the mixed fraction will be close to 0.5. This value is approached asymptotically and a threshold of 0.53 was used for sampling.
Infectious state M. tb objective function formulation
The initial build of iAB‐AMØ‐1410‐Mt‐661 used the original tuberculosis biomass objective function of iNJ661. To build an objective function more specific to in vivo conditions of M. tb, the PM objective function was modified according to gene expression data. The original biomass was split into major metabolite components (Table I). Each component was randomly sampled to determine averages of reaction fluxes. In an iterative manner using linear regression, the components were either removed or weighted and sampling was redone. This process was continued till the sampling results from the specific objective function best fit the gene expression data.
Gene essentiality tests
Gene essentiality tests were completed by systematically setting all reaction fluxes associated with a gene to an upper and lower bound of zero. Flux‐balance analysis was then computed to determine the change, if any, of phenotypic state. Only tuberculosis genes were tested in the PM and the HPM. We set 0.2 of the maximal biomass as the threshold for essentiality.
Building M. tb infection specific models
We used gene expression profiling data from ex vivo macrophages infected with M. tb (Thuong et al, 2008). The study had three types of macrophages: derived from nurses with latent tuberculosis and treated patients of pulmonary and meningeal tuberculosis. The expression data was used to build three context‐specific models using the PANP software package and GIMME algorithm. All reactions were deemed active in the models that were able to carry a non‐zero steady‐state flux. This was performed using flux variability analysis, in which the minimum and maximum fluxes are independently calculated for each reaction in an iterative manner (Reed and Palsson, 2004).
We thank Richard Que for helping in putting together the Supplementary information. This work was supported in part by National Institutes of Health Grants GM068837 and Y1‐AI‐8401‐01, and the National Science Foundation IGERT Plant Training Grant DGE‐0504645.
Conflict of Interest
The authors declare that they have no conflict of interest.
iAB‐AMØ‐1410 SBML provided in SBML format
Spreadsheet detailing all 288 metabolic functions and which are functional in the macrophage.
Spreadsheet detailing the up‐ and down‐regulated reactions after integration.
Gene Essentiality tests for the M. tb portion of the networks.
Listing of the active and inactive reactions in the three infectious states of M. tb studied.
Spreadsheet containing reactions and metabolites of the macrophage model. Literature sources for manual curation and confidence levels of reactions are also provided.
Spreadsheet with listing of reactions and metabolites in the final integrated model.
Calculations for building biomass maintenance function for macrophage.
Table containing exchange metabolites of the macrophage model.
Table containing exchange metabolites of the M. tb and integrated host‐pathogen models.
Listing of all additional and updated reactions in Recon 1 as discussed in the methods section.
Figure S1: Oxygen ‐ Mycobacterium tuberculosis biomass curve
M. tb growth during infection in silico is dependent on oxygen uptake. There is a linear dependence on oxygen until about 13% of the measured in vitro oxygen uptake is reached. Any additional oxygen does not increase the in silico growth rate.
This is an open‐access article distributed under the terms of the Creative Commons Attribution License, which permits distribution, and reproduction in any medium, provided the original author and source are credited. This license does not permit commercial exploitation without specific permission.
- Copyright © 2010 EMBO and Macmillan Publishers Limited