Genome‐scale metabolic reconstructions can serve as important tools for hypothesis generation and high‐throughput data integration. Here, we present a metabolic network reconstruction and flux‐balance analysis (FBA) of Plasmodium falciparum, the primary agent of malaria. The compartmentalized metabolic network accounts for 1001 reactions and 616 metabolites. Enzyme–gene associations were established for 366 genes and 75% of all enzymatic reactions. Compared with other microbes, the P. falciparum metabolic network contains a relatively high number of essential genes, suggesting little redundancy of the parasite metabolism. The model was able to reproduce phenotypes of experimental gene knockout and drug inhibition assays with up to 90% accuracy. Moreover, using constraints based on gene‐expression data, the model was able to predict the direction of concentration changes for external metabolites with 70% accuracy. Using FBA of the reconstructed network, we identified 40 enzymatic drug targets (i.e. in silico essential genes), with no or very low sequence identity to human proteins. To demonstrate that the model can be used to make clinically relevant predictions, we experimentally tested one of the identified drug targets, nicotinate mononucleotide adenylyltransferase, using a recently discovered small‐molecule inhibitor.
Malaria remains one of the most severe public health challenges worldwide (WHO, 2008). Although several available drugs have been successful in controlling malaria in the past, most of them are rapidly losing efficacy due to acquired drug resistance in the most lethal causative agent, Plasmodium falciparum (Mackinnon and Marsh, 2010). This creates an urgent need for new drugs and treatments, as well as better understanding of the parasite physiology. With this in mind, we built a genome‐scale flux‐balance model of the P. falciparum metabolism. Given the complex life cycle of Plasmodium, the flux‐balanced model is of direct relevance to the ongoing search to identify new therapeutic drug targets. The model can be used to explore diverse metabolic states of the parasite and identify essential metabolic genes in the context of known alternative pathways (Oberhardt et al, 2009).
The reconstructed model, which is based on Plasmodium‐specific databases, genomic annotations, and literature reports, includes 366 genes, 1001 reactions, 616 metabolic species, and 4 cellular compartments. We applied flux‐balance analysis (FBA) (Orth et al, 2010) to identify the genes and reactions that are required to produce a set of necessary biomass components. Interestingly, compared with the yeast metabolic network (Duarte et al, 2004), a model eukaryote with a similar genome size, the Plasmodium network has a significantly higher proportion of essential genes; we confirmed this result using a comparative analysis of known gene knockouts in the two microbes. This low level of genetic robustness, which is likely due to the parasitic lifestyle, suggests that many metabolic genes of the parasite can be used as effective drug targets. Indeed, based on the in silico analysis we identified 40 essential P. falciparum genes with no or very little sequence identity to their human homologs.
We used a recently described small‐molecule inhibitor (compound 1_03; Sorci et al, 2009) to experimentally verify one of the enzymes identified as essential: nicotinate mononucleotide adenylyltransferase (NMNAT; Figure 2A). This enzyme, and the corresponding NAD synthesis and recycling pathway, have been recently used for anti‐microbial development (Magni et al, 2009). However, to the best of our knowledge, they have not been used against P. falciparum. The compound 1_03 was able to completely block host cell escape and reinvasion by arresting parasites in the trophozoite growth stage (Figure 2B). These results demonstrate that the inhibitory compound may be a good starting lead for new anti‐malarials.
Importantly, the metabolic model of the parasite can be also used to integrate various genomic data, such as gene expression (Oberhardt et al, 2009). To illustrate these possibilities, we applied gene‐expression data as constraints for the flux‐balance model (Colijn et al, 2009) in order to predict changes in metabolic exchange fluxes. We found that the model was able to correctly predict the changes in external metabolite concentrations (Olszewski et al, 2009) with about 70% accuracy (Figure 3). The availability of a human metabolic network reconstruction (Duarte et al, 2007) would allow, in the future, to analyze the combined parasite–host network, which would deepen understanding of the P. falciparum metabolic vulnerabilities.
Future improvements of the presented P. falciparum metabolic model, for example incorporation of missing activities and yet undiscovered pathways, will lead to a better understanding of parasite physiology. Ultimately, the improved understanding should significantly accelerate the identification and development of desperately needed new drugs against this devastating disease.
In the paper we present a metabolic reconstruction and flux‐balance analysis (FBA) of Plasmodium falciparum, the primary agent of malaria. The compartmentalized metabolic network of the parasite accounts for 1001 reactions and 616 metabolites. Enzyme–gene associations were established for 366 genes and 75% of all enzymatic reactions.
The model was able to reproduce phenotypes of experimental gene knockout and drug inhibition assays with up to 90% accuracy. The model also can be used to efficiently integrate mRNA‐expression data to improve the accuracy of metabolic predictions.
Using FBA of the reconstructed metabolic network, we identified 40 enzymatic drug targets (i.e. in silico essential genes) with no or very low sequence identity to human proteins.
We experimentally tested one of the identified drug targets, nicotinate mononucleotide adenylyltransferase, using a recently discovered small‐molecule inhibitor.
Malaria is an ancient disease, which can be dated back to 2800 bc (Nerlich et al, 2008), and remains one of the most severe public health challenges worldwide. Currently, about half of the Earth's population is at risk from this infectious disease according to the World Health Organization (WHO, 2008). Malaria inflicts acute illness on hundreds of millions of people worldwide and leads to at least one million deaths annually (Baird, 2005; WHO, 2008). It ranks as a leading cause of death and disease in many developing countries, where the most affected groups are young children and pregnant women (WHO, 2008). The disease is transmitted to humans by the female Anopheles mosquito and is caused by at least five species of Plasmodium parasites. The life cycle of the parasite is highly complex and includes various hosts and tissue types. During a blood meal, sporozoites are transmitted from the mosquito to humans and initiate infection in the liver where they reproduce prolifically but are asymptomatic. In the next stage of infection, the parasites are released from the liver cyst into the bloodstream in the form of merozoites, where they invade red blood cells (RBCs) and reproduce asexually (Aly et al, 2009). The destruction of RBCs coupled with the significant load imposed on the host metabolism is ultimately responsible for the major clinical symptoms of malaria, which are often fatal (Haldar and Mohandas, 2009).
Although several anti‐malarial drugs are currently available, most of them are losing efficacy due to acquired drug resistance in the most lethal causative agent, Plasmodium falciparum (Wongsrichanalai et al, 2002; Mackinnon and Marsh, 2010). The loss of drug efficiency in resistant strains poses a great threat to malaria control and has been linked to increases in worldwide malaria mortality (Hyde, 2007). There is an urgent need for new anti‐malarial drugs coupled with better administration strategies. Understanding the molecular mechanisms and interactions of the parasite's cellular components is essential for identification of new drug targets, especially given the difficulties associated with in vivo drug testing (Liu et al, 2008).
Various systems biology approaches have been applied to improve our understanding of P. falciparum physiology and to facilitate drug development (Dharia et al, 2010). The sequencing of the P. falciparum genome has provided researchers with a complete collection of parasite proteins and likely regulatory interactions (Gardner et al, 2002; Carvalho and Menard, 2005). Several large‐scale transcriptome (Bozdech et al, 2003; Le Roch et al, 2003; Llinas et al, 2006; Hu et al, 2010), proteome (Florens et al, 2002; Lasonder et al, 2002, 2008) and metabolome (Olszewski et al, 2009) analyses have been conducted in order to dissect functional interactions and define essential biological pathways. In addition to experimental studies, several databases have been developed to integrate functional knowledge of the parasite and its metabolism. For example, PlasmoCyc is an integrated database that links genomic data, protein annotation, enzymatic reactions, and pathway information (Yeh et al, 2004); the Malaria Parasite Metabolic Pathway Database (MPMP), on the other hand, is a manually curated resource that assembles annotated enzymes into likely metabolic pathways (Ginsburg, 2010).
A stoichiometric representation of metabolism can be effectively used to study functional properties of biochemical networks using a growing number of computational methods (Price et al, 2004). For example, flux‐balance analysis (FBA) considers steady‐state distributions of metabolic fluxes satisfying a set of biophysical constrains, such as bounds and mass balance of fluxes (Orth et al, 2010). Given the applied constraints, a likely distribution of fluxes in the network can be obtained by maximizing an appropriate objective function (e.g. biomass production) (Varma and Palsson, 1994) or applying minimal perturbation principles (Segre et al, 2002; Shlomi et al, 2005). The analysis of flux‐balanced genome‐scale metabolic networks is useful not only for the discovery of essential genes and potential drug targets, but also as a tool to better understand species‐specific biology (Breitling et al, 2008; Oberhardt et al, 2009). For example, among other applications, these models have been used to identify minimal media requirements for growth (Chavali et al, 2008), explore metabolic weaknesses in bacterial pathogens (Navid and Almaas, 2009), integrate gene expression and other types of data (Colijn et al, 2009), and investigate objective functions important under different growth conditions (Schuetz et al, 2007). Given the complex life cycle of Plasmodium, a flux‐balanced model of this organism is of direct relevance to the ongoing search to identify new therapeutic drug targets.
In this study, we reconstructed the genome‐scale flux‐balanced metabolic network of P. falciparum and used it to perform a systems‐level analysis of the parasite's metabolism. On the basis of in silico gene deletions, we identified potential new anti‐malarial drug targets with low sequence identity to human proteins. One of these targets, nicotinate mononucleotide adenylyltransferase (NMNAT), was experimentally tested in a growth inhibition assay using a recently discovered small‐molecule inhibitor. We also illustrate, using a previously published methodology, how the reconstructed metabolic model can be used to integrate flux analysis with expression data to more accurately simulate the physiology of this complex eukaryotic pathogen.
Scale of the reconstructed flux‐balanced metabolic network
The reconstructed flux‐balanced model is based on gene‐reaction associations reported in public domain databases as well as on a careful literature analysis. We used well‐curated microbial metabolic models and enzyme databases to determine the stoichiometry of most reactions. To produce a functional reconstruction, we also searched the literature for missing steps necessary for the model to produce a set of required biomass components (see Materials and methods). The model accounts for 366 genes, corresponding to 7% of all genes identified in P. falciparum. Compared with 61 metabolic models of various organisms compiled by Feist et al (2009), our model ranks 10th in terms of the smallest number of genes. Not surprisingly, the other metabolic models with small gene numbers include many parasitic/symbiotic species, such as Mycoplasma genitalium, Buchnera aphidicola, Haemophilus influenzae, and Helicobacter pylori. The P. falciparum network also includes 616 metabolites and 1001 reactions, 657 of which are metabolic transformations (Table I). In addition, there are 233 reactions corresponding to transport between different cellular compartments and 111 input–output exchange reactions that allow extracellular metabolites to enter and end products to be excreted from the network.
The metabolic reconstruction includes four distinct compartments: parasite cytosol, mitochondria, apicoplast (a non‐photosynthetic plastid), and the extracellular space (representing the host cell cytosol and host serum). The majority of all reactions (50%) occur in the cytosol. The apicoplast accounts for 10% of all reactions, such as the synthesis of isopentenyl diphosphate, fatty acids, and heme (Ralph et al, 2004). A special reaction is included in the model to account for the biomass components and essential metabolites needed for growth (Supplementary Table S1). Supplementary information provides a complete list of all network reactions and metabolite abbreviations.
Excluding metabolite‐exchange reactions, 74% of the reactions in the model are directly associated with P. falciparum genes, which compares well to other models of eukaryotes such as the iND750 yeast model (70%) (Duarte et al, 2004) and the iAC560 model for Leishmania major (63%) (Chavali et al, 2008). The remaining reactions include spontaneous transformations that can proceed without enzymatic catalysis and reactions required for the proper functioning of the metabolic model. Intracellular and inter‐compartmental transport reactions, most of which are not currently associated with any gene, account for about 6% and 15% of all reactions in the model, respectively (Figure 1A). Most of the transporter proteins in Plasmodium spp. are currently uncharacterized. However, it is well established that the parasite significantly modifies the permeability of the host cell membrane (Kirk et al, 1999; Martin et al, 2005) and several metabolic processes occur across different organelles. For instance, such metabolic pathways as heme biosynthesis and antioxidant defense have been shown to involve both host and parasite enzymes localized to multiple intracellular compartments (Bonday et al, 1997; Koncarevic et al, 2009). Given the importance of metabolite exchange, many transport reactions were included in the model, although the identities of the corresponding genes remain unknown (Figure 1A).
Transferases and hydrolases comprise the largest fraction of enzymatic reactions in the network (Figure 1B). In terms of specific metabolic processes (Figure 1C), most reactions in the network are related to lipid metabolism, followed by transport and exchange reactions. In comparison with Saccharomyces cerevisiae, a free‐living eukaryote of similar genome size, the most significant metabolic difference is the fraction of reactions involved in amino‐acid metabolism (Figure 1C). About 20% of the reactions in the iND750 yeast metabolic network (Duarte et al, 2004) are responsible for amino‐acid pathways; in contrast, this fraction is only 7% in P. falciparum. Amino‐acid biosynthesis pathways are absent in P. falciparum metabolism because of the unique ability of the parasite to catabolize the erythrocyte hemoglobin (Francis et al, 1997) and to scavenge free amino acids from the host serum (human stages) or hemolymph (mosquito stages).
Analysis of in silico single and double gene deletions
We simulated gene deletions using FBA of the reconstructed metabolic network. Even though sugars other than glucose do not support P. falciparum growth in culture (Saito et al, 2002), in performing the in silico deletions we initially allowed all exchange reactions to carry non‐zero metabolic fluxes, thereby permitting the potential import and utilization of other hexoses. Purines, such as inosine and adenosine, which are not normally included for in vitro culture but can be imported in vivo (LeRoux et al, 2009), were also made available to the network. We note that genes identified as essential when all exchange reactions are allowed will also be essential under more specific (constrained) conditions. The phenotypic effects of in silico deletions were classified into four groups: lethal, growth reducing (growth between 0 and 95% of the wild‐type network), slight growth reducing (growth between 95 and 100%), and with no effect. About 15% of all single gene deletions (Table II) were lethal, ∼1% were growth reducing, and 3.5% were slightly growth reducing.
Out of 366 genes in the P. falciparum metabolic network, 55 genes were predicted to be essential for growth (Supplementary Table S2). To assess the accuracy of these predictions, we compiled a list of experimentally validated gene knockouts and phenotypes resulting from targeted inhibitions of enzymatic activities with drugs (Table III). In the computational analysis, we assumed that the drug treatments resulted in a complete inhibition of targeted enzymatic activities. In this way, the available drug phenotypes were simulated with computational deletions of the corresponding genes. In total, 14 metabolic gene knockouts and 25 drug inhibition phenotypes for genes were retrieved from the literature for P. falciparum and Plasmodium berghei, a murine malaria parasite commonly used in experimental studies (Janse and Waters, 1995; Janse et al, 2006).
The FBA analysis was able to achieve 100% accuracy for predictions of both essential and non‐essential gene knockouts (14 cases in total). In contrast, about 70% accuracy was achieved for phenotypes resulting from drug inhibitions of metabolic enzymes. Interestingly, all mispredicted drug phenotypes (eight cases) involved genes for which the computational analysis predicted a non‐zero growth phenotype, whereas corresponding drugs were lethal to the parasite in experimental studies. In three cases, inconsistencies between the FBA predictions and experimental drug inhibitions can be explained by considering functions that are not explicitly represented in our model. These included the degradation of spontaneously forming toxic metabolites (e.g. methylglyoxal), and the synthesis of metabolites that are involved in the progression between the intraerythrocytic stages (e.g. sphingolipid, ceramide) (Hanada et al, 2002). The remaining discrepancies (five cases) can be resolved by taking into account specific literature‐based evidence (Table III, green rows), that is, by considering nutrient availabilities and directionality of exchange reactions. Interestingly, in one case, the source of discrepancy between the model and experiments was clearly related to off‐target drug effects. Specifically, the inhibitor of enoyl‐acyl carrier reductase (FabI), triclosan, has been shown to kill P. falciparum in vitro and in vivo (Surolia and Surolia, 2001) despite the fact that its presumed target, PfFabI, can be deleted with no apparent blood‐stage phenotype (Surolia and Surolia, 2001; Yu et al, 2008; Vaughan et al, 2009); this deletion phenotype is correctly predicted by our model.
For 15 metabolic genes identified in our analysis as essential, knockout experiments or drug inhibition assays in P. falciparum/P. berghei are already available in the literature (see Table III; Supplementary Figure S2). The remaining predictions include 24 genes coding for proteins with relatively low sequence identity (20–40%) to human transcripts (Supplementary Figure S1; see Materials and methods), and 16 genes with no significant sequence identity to any human protein (BLAST E‐value>10−2). This last group comprises six genes associated with isoprenoid metabolism, three genes involved in nucleotide metabolism, and genes related to CoA, shikimate, and folate biosynthesis (Table IV). Nine of the genes with no homology to human proteins are homologous to plant proteins; this is consistent with the essential functions of apicoplast‐associated genes in the Apicomplexa (Lim and McFadden, 2010). These 40 enzymes are of immediate interest as potential drug targets, as low homology to human proteins suggests that side effects for drugs targeting these enzymes may be minimized or avoided. Interestingly, 5 of the 16 enzymes with no detectable homology correspond to enzymatic activities (Enzyme Commission (EC) numbers) that are unlikely to be present in human metabolism according to the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al, 2008), HumanCyc (Romero et al, 2005), and UniProt (Uniprot.Consortium, 2010) databases. Among the predicted essential genes with low but significant sequence identity to human transcripts, only aminodeoxychorismate lyase/synthetase (22.214.171.124, 126.96.36.199, PFI1100w) is associated with EC numbers not reported in human metabolism. In addition to genes predicted as essential, nine internal metabolic reactions with no associated network genes (orphan reactions) were also predicted to be essential for growth. Four of these reactions are associated with the shikimate biosynthetic pathway; three with ubiquinone metabolism, one with nicotinamide, and one with porphyrin metabolism (Supplementary Table S3).
One metabolic pathway of significant interest in the parasite is the mitochondrial tricarboxylic acid (TCA) cycle. In most free‐living microbes, this pathway fully oxidizes available carbon sources to carbon dioxide, in the process generating high‐energy phosphate bonds (ATP/GTP). Within the Plasmodium species, however, the nature and function of the TCA cycle remains unclear (van Dooren et al, 2006). In the malaria parasites, the sole pyruvate dehydrogenase complex, which normally provides the key link between glycolysis and TCA metabolism, localizes not to the mitochondrion but to the apicoplast. In that compartment it is likely to be used primarily to generate acetyl‐CoA for lipogenesis (Foth et al, 2005). Incorporating this fact into our model, we find that the only TCA cycle enzyme predicted to be essential is the 2‐oxoglutarate dehydrogenase complex. This enzyme converts 2‐oxoglutarate into succinyl‐CoA, which is required for heme biosynthesis. Intriguingly, metabolic‐labeling experiments indicate that the TCA cycle of the parasite also reduces 2‐oxoglutarate to malate, generating acetyl‐CoA from citrate cleavage (Olszewski et al, 2010). A non‐zero flux through this reaction is observed in our model when additional constraints are applied to mitochondrial transport reactions.
We also extended the computational analysis of essential Plasmodium genes to pairs of genes that are not essential on their own, but are lethal if deleted simultaneously, that is synthetic lethal enzyme pairs with non‐trivial genetic interactions (Dixon et al, 2009) (see Table II). In total, deletion of 16 gene pairs gave rise to such synthetic lethality in the unconstrained model. The enzymes that were essential upon double deletions participate in glycolysis, metabolism of nucleotides, lipids, porphyrin, the pentose phosphate cycle, and transport of NO2 and phosphate (Supplementary Table S4).
The analysis of essential genes in the S. cerevisiae metabolic network iND750 (Duarte et al, 2004) can be used to put the results of the P. falciparum in silico deletions into an appropriate context. To make the proper comparison, we only considered deletions of genes carrying non‐zero metabolic fluxes in wild‐type models of both networks; the focus on non‐zero fluxes is necessary to prevent the difference in network sizes (S. cerevisiae 750 genes/1266 reactions, P. falciparum 366 genes/1001 reactions) from biasing the results. When both networks were allowed to simultaneously use all carbon sources, the fraction of essential genes associated with non‐zero fluxes in P. falciparum was 37%, whereas in S. cerevisiae it was 5% (Fisher's exact test, P‐value<10−10). The fraction of essential genes was also significantly higher in the parasite when single carbon sources were used in both models. For example, when glucose was used as the sole source of carbon, 50% of genes were essential in the parasite versus 26% in yeast (P‐value=10−6). To understand whether the significantly higher fraction of essential genes in P. falciparum arises exclusively from a smaller number of isoenzymes in that network (111 in P. falciparum versus 276 in S. cerevisiae) or if it is also related to the inherent differences in the networks’ architectures, we performed deletions of all isoenzymes associated with metabolic reactions, instead of individual gene deletions. As a result, when networks used all carbon sources, 39% of the reactions with non‐zero fluxes were essential in P. falciparum, compared with only 6% in S. cerevisiae (P‐value<10−10). In a glucose minimal media, 62% of the reactions are essential in the parasite and 47% in yeast (P‐value<10−10). These results demonstrate that a significantly smaller genetic robustness of the parasite's network arises, at least in part, because of a paucity of alternative metabolic pathways (Wagner, 2005). The lower redundancy of the P. falciparum network is likely to be a consequence of the adaptation to the relatively homeostatic and nutrient‐rich environments of the hosts in which it proliferates (Gardner et al, 2002).
Resolution of disagreements between in silico predictions and experimental data
Although the presented metabolic model achieves a high accuracy in predicting phenotypes of the experimental knockouts and drug inhibition assays, it is the disagreement between the model and experiments that often leads to model improvement (Thiele and Palsson, 2009). Thus, it is important to discuss the inconsistencies between modeling and experimental results, which were corrected in the process of model construction. In four cases, the disagreements between the predicted gene essentiality and experimental results reported in the literature were resolved through additional flux constraints. The adjustments included purine nucleoside phosphorylase (PFE0660c), adenosine deaminase (PF10_0289), ornithine decarboxylase (PF10_0322), and lactate dehydrogenase (PF13_0141/PF13_0144). The additional constraints applied to the network were based either on specific experimental conditions or known details of Plasmodium spp. physiology (Table III, green rows). For example, lactate is believed to be the main byproduct of glucose metabolism in P. falciparum (Vaidya and Mather, 2009), and lactate dehydrogenase is an essential enzyme that regenerates nicotinamide adenine dinucleotide (NAD+) from NADH (Berwal et al, 2008). In contrast, in our initial model, pyruvate was exported as the primary glycolysis byproduct and NAD+ was regenerated through the transformation of pyrroline‐5‐carboxylate to proline (EC: 188.8.131.52). As there is no evidence of extensive pyruvate export in Plasmodium, a constraint was added to the corresponding pyruvate exchange reaction. As a result, we observed reduction of pyruvate to lactate, followed by lactate export. In the adjusted model, lactate dehydrogenase carried a non‐zero flux and was correctly predicted to be essential for growth. In the other three cases, constraints on exchange fluxes were added to reproduce the composition of the media used in drug inhibition experiments. Specifically, purine nucleoside phosphorylase has been shown to be essential if hypoxanthine is not available in the environment (Kicska et al, 2002), whereas adenosine deaminase is essential in the absence of both inosine and hypoxanthine (Ho et al, 2009). When the fluxes through the inosine and hypoxantine exchange reactions were set to zero, the experimentally observed knockout phenotypes were reproduced. Similarly, ornithine decarboxylase was correctly predicted to be essential without putrescine in the growth media (Das Gupta et al, 2005).
In two cases, the inability of our initial model to reproduce experimental results was due to reactions involving metabolic dead ends; that is metabolites that are either only produced or only consumed in the network (Reed et al, 2003). In the first case, the model was not able to synthesize spermidine. The synthesis of spermidine from putrescine by spermidine synthase was accompanied by the production of 5‐methylthioadenosine (MTA) (Haider et al, 2005), which was a metabolic dead end in the initial model. Consequently, the spermidine synthesis caused MTA to accumulate, violating the steady‐state assumption of the constraint‐based approach. However, although it is known that in Plasmodium spp. MTA is first converted to 5‐methylthioinosine by adenosine deaminase and then recycled into methionine and hypoxanthine (Ting et al, 2005), not all enzymes involved in these reactions have been fully characterized. We addressed this problem by including in the model the PfADA and PfPNP activities, responsible for the hypoxanthine synthesis from MTA (Ting et al, 2005), and an additional hypothetical reaction that converts the resulting by‐product, 5‐methylthioribose‐1‐PO4, to methionine in order to represent the methionine salvage pathway. The second case was related to the folate biosynthesis pathway. In this pathway, the reaction catalyzed by 6‐pyruvoyltetrahydropterin synthase (184.108.40.206), in addition to the folate biosynthesis intermediate 6‐hydroxymethyl‐7,8‐dihydropterin, is known to produce a small amount of 6‐pyruvoyl‐5,6,7,8‐tetrahydropterin (6PTHP) (Hyde et al, 2008). Initially, both products were included in the same reaction and, because 6PTHP represented another metabolic dead end, the cell was not able to synthesize folate. We resolved this problem by including in the model separate reactions for each of the two alternative products.
The total number of remaining dead‐end metabolites in the final model (266) is comparable to that of other recently published genome‐scale metabolic networks; for example the iBsu1103 model for Bacillus subtilis (270 dead‐end metabolites; Henry et al, 2009), the iAC560 model of L. major (261; Chavali et al, 2008), or the iND750 S. cerevisiae model (194; Duarte et al, 2004). The remaining dead ends in the Plasmodium network include 109 metabolites that are consumed but are not currently produced or imported in the model, 80 metabolites that are produced but not consumed, and 79 metabolites, associated with reversible reactions, that can be either exclusively consumed or produced. As protein synthesis is not explicitly included in the metabolic model, a significant number (42) of the remaining dead‐end compounds correspond to tRNAs; this compares to 68 dead‐end tRNAs in the iND750 yeast network. Among the other functional categories associated with a large number of the metabolic dead ends are lipid metabolism (45%), transport reactions (15%), the metabolism of carbohydrates (5%), amino acids (8%), and nucleotides (5%) (Supplementary information). More precise measurements of the P. falciparum biomass composition, for example a detailed lipid composition of the parasite membrane (Hsiao et al, 1991; Mi‐Ichi et al, 2006), can be used in the future to significantly shrink the pool of the remaining dead‐end metabolites.
Validation of the predicted drug target nicotinate nucleotide adenylyltransferase
The most urgent motivation for flux‐balance reconstruction of the pathogen metabolism is to facilitate drug development. To illustrate the potential of the model to make clinically relevant predictions, we experimentally tested a predicted target for which candidate drugs have been reported in other microbial species. The ideal drug target will be an essential enzyme with no homolog in the human genome and in a pathway not currently targeted by any pharmaceutical. On the basis of these criteria, we selected for validation NMNAT (PlasmoDB ID PF13_0159, EC 220.127.116.11) (Table IV). This enzyme, a member of the plasmodial NAD synthesis and recycling pathway, catalyzes the conversion of nicotinic acid mononucleotide to nicotinic acid adenine dinucleotide (Figure 2A). NMNAT has recently been the focus of novel anti‐microbial agent development because of structural and metabolic differences between the enzyme in microbial and human cells (Magni et al, 2009). As NAD(P) is one of the most promiscuous redox molecules in metabolism, as well as a cofactor for important histone regulatory proteins such as sirtuins (Merrick and Duraisingh, 2007), inhibition of NAD(P) synthesis and recycling should have a profound impact on parasite metabolism. However, to the best of our knowledge, this pathway has not been previously targeted by pharmaceutical interventions in P. falciparum.
Recently, Sorci et al (2009) used a combination of in silico structure modeling and enzyme inhibition assays to identify several classes of small molecules that inhibit bacterial (Escherichia coli and B. subtilis) but not human NMNAT. Several of these drug candidates were able to inhibit bacterial growth in culture. We tested two of the designed candidate compounds (1_03 and 3_02), representing two different chemotypes, for their ability to inhibit P. falciparum growth using both the SYBR Green I fluorescence assay (Smilkstein et al, 2004) to measure DNA synthesis, and microscopic examination of morphological effects. The two compounds were tested at a range of concentrations for growth inhibitory effects in nicotinamide‐free culture medium, so as not to rescue any metabolic blocks induced by the drugs. Nicotinamide removal has been previously shown not to affect normal growth in culture (Divo et al, 1985), which we confirmed before running our growth assay experiments (data not shown). Although the compound 3_02 did not significantly affect parasite's growth at moderate concentrations (MIC50 >100 μM), the compound 1_03 exerted an inhibitory effect in a growth assay (MIC50 = 50 μM; Supplementary Figure S2) comparable to that previously observed for bacteria (MIC50 >80 μM for E. coli, MIC50 =10 μM for B. subtilis). At 100 μM, the compound 1_03 completely blocked host cell escape and reinvasion by arresting parasites in the trophozoite growth stage (Figure 2B). Importantly, the human NMNAT isoforms are insensitive to the compound at least up to the concentrations used in our assay (Sorci et al, 2009). This suggests that the parasite NMNAT enzyme and, more generally, the NAD(P) synthesis pathway are indeed potentially effective and druggable targets. The experimental results also highlight the ability of our model to identify promising candidates for pharmaceutical intervention.
Prediction of metabolite concentration changes based on expression data
Genome‐scale metabolic networks can be used not only to predict the effects of gene deletions, but also as a tool for the integration of diverse genomic and physiological data (Breitling et al, 2008; Oberhardt et al, 2009). For example, information on nutrient availability, uptake rates, and maximal rates of internal reactions can be used to further constrain the space of feasible metabolic fluxes. To illustrate the ability of the model to combine genomic data, we investigated whether available gene‐expression data sets can be used to predict shifts in concentrations of external metabolites caused by the P. falciparum exchange fluxes at different developmental stages. Investigation of the exchange fluxes is essential for understanding perturbations caused by parasitic infections in the metabolic state of their host tissues, and, consequently, main mechanisms of pathogenesis. As the FBA operates in the space of the fluxes and not in the space of metabolic concentrations, the model cannot be directly used to predict absolute concentration changes. Nevertheless, it is possible to use the model to investigate the direction of concentration changes for external metabolites, following the simple logic that an increase in the uptake rate or decrease in the excretion (output) rate should lead to a drop in the concentration of the corresponding external metabolite; similarly, an uptake rate decrease or excretion rate increase should increase the metabolite concentrations.
Although gene‐expression level does not perfectly correlate with the flux through the associated enzyme (Daran‐Lapujade et al, 2004; Shlomi et al, 2008), the recent study by Colijn et al (2009) demonstrated that mRNA abundance data, if used as additional constraints on maximal reaction fluxes, can significantly improve stoichiometric model predictions. For instance, if the expression level of a particular enzyme is low, it is unlikely that the enzyme will be used by the metabolic network to carry a large flux. Consequently, it should be possible to use gene‐expression data to obtain a more accurate view of the in vivo metabolic state. To test this, we used DNA microarray results collected from synchronized cultures of the P. falciparum 3D7 strain during the RBC phase of the parasite's life cycle (Llinas et al, 2006; Olszewski et al, 2009). The expression data were collected at the ring, trophozoite, and schizont developmental stages (see Materials and methods). Following Colijn et al (2009), the maximum flux allowed through enzymes was constrained proportionally to the relative expression level of the corresponding genes.
We compared the accuracy of our predictions to the experimentally measured metabolic changes in Plasmodium‐infected RBCs (Olszewski et al, 2009). In Figure 3, we show the predicted and experimentally measured changes, indicating either an increase or decrease in metabolic concentrations for the transition from the ring to trophozoite and from trophozoite to schizont stages. The predicted shifts in metabolic concentrations agree with the experimental results in 70% of the measurements (binomial, P‐value=9 × 10−4). In addition, we found a significant correlation between the magnitudes of the change in metabolite concentrations and the predicted flux values (Pearson's correlation: 0.34, P‐value=6 × 10−3, Spearman's correlation: 0.25, P‐value=0.04).
In order to further investigate the statistical significance of the results, we repeated flux predictions after randomly shuffling expression values between P. falciparum genes. In only 2% of these random trials, the accuracy of the predictions made with the shuffled data were higher than those obtained using the original expression values (Supplementary Figure S3). To explore the effects of multiple optimal FBA solutions (Mahadevan and Schilling, 2003) on the prediction accuracy, we used the centering hit‐and‐run algorithm (Kaufman and Smith, 1998), implemented in the COBRA toolbox (Becker et al, 2007), to randomly sample the solution space associated with the expression constraints. The 70% accuracy value, obtained for a single solution, is close to the mean of solutions sampled from alternative optima (mean 0.69, s.d. 0.05; see Supplementary Figure S3). Moreover, there is a significant difference (Mann–Whitney U, P‐value<10−10) between the results for randomized expression values and those based on the multiple alternative optima. These results illustrate the ability of the model, with appropriate constraints, to predict physiological changes unrelated to gene knockouts. It also suggests that expression and metabolomics measurements, which are being rapidly accumulated for various stages of parasite growth (Winzeler, 2008; Kafsack and Llinas, 2010), can be integrated with the model to gain a better understanding of the P. falciparum physiology.
The presented flux‐balanced model can serve as a valuable tool for quantitative predictions of P. falciparum metabolic states under various growth conditions and perturbations. The results of in silico gene deletions demonstrate that the model achieves high accuracy in reproducing available experimental measurements. In addition, our analysis suggests several dozen essential metabolic targets for therapeutic intervention. Although several studies that assemble and analyze plasmodial metabolic pathways have been performed previously (Yeh et al, 2004; Ginsburg, 2006, 2009, 2010), our contribution is important because the genome‐scale model can be used to investigate and predict genetic perturbations from a network‐level perspective.
Interestingly, our results suggest a limited degree of robustness in the P. falciparum network, which should lead to a relatively high‐success rate for inhibitors of metabolic genes. It is possible that the small robustness of the reconstructed model, for example in comparison with the yeast metabolic network, is due primarily to unannotated P. falciparum genes without significant homology to known enzymes in other organisms. To investigate this possibility further, we used the available collection of single metabolic gene knockouts/inhibitions in P. falciparum or P. berghei (Table III) and all metabolic gene knockouts in S. cerevisiae (Giaever et al, 2002). We calculated the fraction of orthologs for essential metabolic knockouts in the parasite, which are also essential in yeast, and, vice versa, the fraction of orthologs for essential metabolic knockouts in yeast, which are essential in the parasite (see Supplementary Table S5). Interestingly, while only about half of the orthologs for essential metabolic genes in P. falciparum are also essential in S. cerevisiae, all essential metabolic genes we analyzed in yeast are essential in the parasite (Supplementary Table S5; Fisher's exact test, P=0.04). This result independently supports the conclusion of the flux‐balance simulations about the relatively small robustness of the P. falciparum network.
We anticipate several immediate extensions of our work. First, the presented network can be used for effective integration of multiple genomic data types. For example, known regulatory interactions can be incorporated into the model (Covert and Palsson, 2002). Accurate measurements of gene expression (Hu et al, 2010), key protein–DNA regulatory interactions (De Silva et al, 2008), and post‐translational modifications (Chung et al, 2009) in the parasite will be especially important for modeling the dynamic behavior of the network under varying environmental conditions. Second, it will be important to model exchanges and interactions between the metabolic networks of the parasite and its hosts. The analysis of the combined parasite–host metabolic network should significantly improve understanding of the P. falciparum vulnerabilities. For example, several host cell enzymes are actively used by the parasite during its life cycle (Dhanasekaran et al, 2004; Ting et al, 2005). Although we did not consider these human enzymes in our analysis, they can be easily included in future applications of the model. The available global flux‐balanced metabolic human network (Duarte et al, 2007), metabolic network specifically active in the liver (Zhao et al, 2010), and well‐curated models of human RBC metabolism (Joshi and Palsson, 1989a, 1989b, 1990a, 1990b; Jamshidi et al, 2001) make such combined analyses possible. Third, it will be interesting to reconstruct stoichiometric metabolic networks for other clinically relevant Plasmodium species (e.g. P. vivax, P. malariae, P. ovale, and P. knowlesi) as well as the important model species P. berghei and Plasmodium yoelii. The comparative analysis of these networks may reveal important physiological and evolutionary differences between Plasmodium spp., and also help in the identification of common metabolic drug targets.
Taking into account the global health burden of malaria, it is essential to develop and implement new effective pharmaceuticals as quickly as possible. Systems biology approaches can be used to significantly facilitate drug identification and development (Yao and Rzhetsky, 2008; Kuhn et al, 2010). To date, we have only begun to see the application of such integrative methods in the context of malaria research (reviewed in Dharia et al, 2010). We believe that the presented network represents an important step in this direction. The experimental validation of a candidate drug, compound 1_03, targeting the parasite NMNAT illustrates the ability of the model to speed up development of novel anti‐malaria targets and pharmaceuticals. Importantly, although the identified compound is available, it has not been previously tested against Plasmodium spp. Even though 1_03 inhibited parasite growth only at relatively high concentrations (MIC50=50 μM), these were comparable to the inhibitory concentrations for the bacteria against which the drug was initially developed (Sorci et al, 2009). The incomplete growth inhibition at lower compound concentrations could be explained by incomplete drug inhibition. Our model predicts linear decrease in the P. falciparum biomass production as the level of NMNAT inhibition increases; for example, 90% inhibition results in 90% growth decrease. As Sorci et al initially screened for compounds that could selectively inhibit pure NMNAT enzyme, these compounds have not been optimized for cellular permeability, accumulation, or other pharmacokinetic parameters, and thus should primarily serve as a structural basis for further malaria drug development.
Future improvements to the reconstructed P. falciparum metabolic network, including adding experimental details for missing activities and precise metabolic measurements necessary to describe the growth‐related objective function, will lead to a better understanding of parasite physiology. Ultimately, such models should significantly accelerate the identification of desperately needed new drug leads against this devastating disease.
Materials and methods
Genome‐scale metabolic reconstruction
The reconstruction process started with the identification of enzyme‐coding genes in the P. falciparum genome. We considered a variety of resources, including PlasmoDB (Aurrecoechea et al, 2009), the MPMP (Ginsburg, 2006), PlasmoCyc (Yeh et al, 2004), and KEGG (Kanehisa et al, 2008). The identified enzymes were mapped to the corresponding metabolic reactions by consulting several well‐studied metabolic models, including the iAF1260 model for E. coli (Feist et al, 2007), the iND750 model for S. cerevisiae (Duarte et al, 2004), the genome‐scale human metabolic network by Duarte et al (2007), and the KEGG database (Kanehisa et al, 2008). On the basis of this set of enzymatic activities and their stoichiometry, we used FBA to see if the network was able to produce a set of basic biomass components; e.g. amino acids, lipids, nucleotides, and cofactors. For each metabolite that the network was unable to synthesize, we searched the literature for relevant publications concerning pathways and genes associated with the metabolite production or transport (relevant publications are included as notes in Supplementary information). Network enzymes were assigned to different cellular compartments based on experimental evidence, when available, and on computationally predicted localization information (Waller et al, 1998; Ralph et al, 2004; Chan et al, 2006; van Dooren et al, 2006). Transport and exchange reactions reported in the literature or in databases such as PlasmoDB and MPMP were initially included in the model. We added additional transport reactions required for production of the biomass components. All metabolic and transport reactions were used to formulate a stoichiometric flux‐balance model (Edwards et al, 1999). The model was improved following an iterative procedure as previously described (Feist et al, 2009; Thiele and Palsson, 2010).
The assembled network was manually inspected and compared with the MPMP (Ginsburg, 2010). Metabolic network gaps (Kharchenko et al, 2006) were identified and included in the assembled network model (Mullin et al, 2006; Quashie et al, 2008). The reactions for which no literature support is available and which are not essential for the biomass production were removed from the network. Additional adjustments related to reaction directionalities and metabolite availabilities were made following the computational analyses described in this work.
As the P. falciparum biomass objective function cannot be completely established based on the available literature, in our calculations we used a modified version of the yeast objective function reported in the iND750 model (Duarte et al, 2004). The objective function modifications included the lipid composition, which was adjusted as reported for Plasmodium (Hsiao et al, 1991), and amino acid and nucleotide compositions adjusted based on the proteome and genome sequences weighted by available expression data (Llinas et al, 2006). In particular, the percent prevalence of each ribonucleotide and amino acid across all open reading frames (ORFs) was calculated as the relative frequency of each monomer; the counts at each ORF were multiplied by the ORF's expression level (when available). The percent prevalence of the dNTPs was derived from the genome A+T content of 80.6%. These percentages were converted to mmol/gDW as described (Chavali et al, 2008). Systems Biology Research Tool (Wright and Wagner, 2008) was used to perform FBA (Edwards et al, 1999) of the network, including single and double in silico deletions of network enzymes. The reconstructed network was able to either synthesize or import all the biomass components presented in Supplementary Table S1. The assembled metabolic model is available as an Excel spreadsheet (Supplementary information) and in the Systems Biology Markup Language (SBML) format (Supplementary information). The SBML model was submitted to the BioModels database (Le Novere et al, 2006) with accession number MODEL1007060000.
Parasite culture, growth, and drug inhibition assays
The cultures of P. falciparum were maintained and synchronized by standard methods (Trager and Jensen, 1976; Lambros and Vanderberg, 1979). Briefly, RBCs, infected by P. falciparum (3D7 strain), were grown in the RPMI 1640 culture medium supplemented with sodium carbonate (2 mg/ml), hypoxanthine (100 μM), Albumax II (0.25%), and gentamycin (50 μg/ml) in a humidified incubator at 5% CO2, 6% O2, and 37°C. The growth synchronizations were carried out by incubating parasite‐infected RBCs in phosphate‐buffered saline (PBS) containing 5% w/v sorbitol for 5 minutes at room temperature, washing once with sorbitol‐free PBS, and resuspending in culture medium.
The compounds 1_03 and 3_02 were acquired from ChemDiv (http://chemdiv.emolecules.com; ChemDiv IDs 8003‐6329 and 5350‐0029, respectively) and resuspended at 100 mM in DMSO. Growth inhibition studies were carried out using the SYBR Green I fluorescence assay (Smilkstein et al, 2004). Briefly, synchronized parasite cultures (early ring stage, 1% parasitemia, 1% hematocrit, 100 μl total volume) were suspended in nicotinamide‐free RPMI 1640 containing 0.1% DMSO and differing concentrations of drug in 96‐well plates. After 72 h incubation, the plates were frozen at −80°C overnight, then thawed and mixed with 100 μl lysis buffer (20 mM Tris–HCl, pH 7.5; 5 mM EDTA; 0.008% w/v saponin; 0.08% v/v Triton X‐100; 1 × SYBR Green I) per well, incubated 1 h at room temperature and quantified using a BioTek SynergyMX plate reader (excitation 488 nm, emission 522 nm). The concentrations 0, 0.1, 1, 5, 10, 50, 100, and 250 μM were tested in triplicate in two independent growth assays.
Using an expression‐constrained network to predict shifts in external metabolite concentrations
In the expression‐constrained flux analysis, we used P. falciparum expression data (Olszewski et al, 2009) as described previously in Colijn et al (2009). Specifically, for genes with available mRNA‐expression data, the maximum flux through the associated metabolic reactions was constrained proportionally to their expression level; with the highest expression value corresponding to the maximum possible metabolic flux. In order to obtain absolute expression values, rather than the ratios between the microarray intensities at a given time point and those of a pooled sample, we multiplied each ratio with the average sum of the median intensity across the full intraerythrocytic developmental cycle (Llinas et al, 2006).
The ring, throphozoite, and schizont developmental stages in cultures were defined for hours 1–18, 19–30, and 31–48 after synchronization with D‐sorbitol, respectively. In the analysis, we used intracellular RBC metabolite concentration data obtained by Olszewski et al (2009). We compared the changes in metabolite abundances between infected and uninfected RBCs, from one development stage to the next. For each metabolite, the predicted concentration changes were considered to agree with experimental data, if the metabolite consumption in the network increased (or the metabolite production decreased) when the experimentally measured metabolite concentration decreased, or alternatively, if consumption of the metabolite decreased (or production increased) when the metabolite concentration increased.
We thank the members of DV and ML groups for stimulating discussions and comments on the manuscript. Research was funded in part by NIH grant GM079759 and BU GC207272NGC to DV, and a National Centers for Biomedical Computing (MAGNet) grant U54CA121852 to Columbia University. ML is funded in part by the Burroughs Welcome Fund, an NIH Director's New Innovators award (1DP2OD001315‐01) and receives support from the Center for Quantitative Biology (P50 GM071508). KO is supported by an NSF Graduate Research Fellowship.
Conflict of Interest
The authors declare that they have no conflict of interest.
Supplementary figures S1–3, Supplementary tables S1–5
This is an open access article under the terms of the Creative Commons Attribution‐NonCommercial‐NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non‐commercial and no modifications or adaptations are made.
- Copyright © 2010 EMBO and Macmillan Publishers Limited